Numerical Assessments of Excess Ice Impacts on Permafrost and Greenhouse Gases in a Siberian Tundra Site Under a Warming Climate

Excess ice that exists in forms such as ice lenses and wedges in permafrost soils is vulnerable to climate warming. Here, we incorporated a simple representation of excess ice in a coupled hydrological and biogeochemical model (CHANGE) to assess how excess ice affects permafrost thaw and associated hydrologic responses, and possible impacts on carbon dioxide and methane (CH4) fluxes. The model was used to simulate a moss-covered tundra site in northeastern Siberia with various vertical initializations of excess ice under a future warming climate scenario. Simulations revealed that the warming climate induced deepening of the active layer thickness (ALT) and higher vegetation productivity and heterotrophic respiration from permafrost soil. Meanwhile, excess ice temporarily constrained ALT deepening and thermally stabilized permafrost because of the highest latent heat effect obtained under these conditions. These effects were large under conditions of high excess ice content distributed in deeper soil layers, especially when covered by moss and thinner snow. Once ALT reached to the layer of excess ice, it was abruptly melted, leading to ground surface subsidence over 15–20 years. The excess ice meltwater caused deeper soil to wet and contributed to talik formation. The anaerobic wet condition was effective to high CH4 emissions. However, as the excess ice meltwater was connected to the subsurface flow, the resultant lower water table limited the CH4 efflux. These results provide insights for interactions between warming climate, permafrost excess ice, and carbon and CH4 fluxes in well-drained conditions.


INTRODUCTION
The warming climate has resulted in changes in the Arctic system. A representative change is permafrost warming (Biskaborn et al., 2019), which is closely related to changes in ecological, hydrological, and biogeochemical functions. The changes are in turn fed back to global climate through altered water, energy, and carbon fluxes. For example, permafrost is estimated to store approximately 1,300 Pg of carbon, which is considerably larger than the current atmospheric carbon pool (Hugelius et al., 2014). If thawed and mobilized, this carbon could become a major source of greenhouse gas emissions (Schuur et al., 2008). Field observations have monitored increases in carbon dioxide (CO 2 ) release from permafrost soils to the atmosphere (Turetsky et al., 2019(Turetsky et al., , 2020, and numerical models also projected the deepening active layer thickness (ALT)-induced increase in organic carbon thawing under future warming scenarios (Koven et al., 2011(Koven et al., , 2015Lawrence et al., 2012;Nitzbon et al., 2020). These results imply that future climate conditions will accelerate the permafrost thaw through positive feedback.
The ice-rich permafrost, formed by ice wedge and Yedoma, is vulnerable to climate change (Jorgenson et al., 2006;Liljedahl et al., 2016). The warming temperature enhances permafrost degradation and melting of ground ice, inducing landscape changes such as subsidence and thermokarst. The thermokarst processes can yield spatially heterogeneous thaw lakes and thermo-erosional gullies, influenced by the microtopography of the terrain. The ground-ice melt supplies water to the thaw lakes in conjunction with snowmelt water and precipitation, so that the thaw lakes tend to continuously grow once initiated (Fedorov et al., 2014), causing severe permafrost degradation within a few years or decades (Ulrich et al., 2017). Thaw lake growth significantly alters landscape scales of water and energy fluxes. In situ measurements and remote sensing observations in several Arctic sites have documented the thermokarst-associated changes in local and regional hydrology (Liljedahl et al., 2016), with implications for a range of ecosystem functions such as, for example, the decomposition of organic carbon (Lara et al., 2015).
The contribution of thermokarst landscape evolution to hydrological processes in a warming climate is likely further amplified by ground subsidence inducing greater snow accumulation, reduced winter heat loss, and increased ground ice melt (Aas et al., 2019). A variety of numerical modeling studies have addressed the hydrologic processes changed by permafrost degradation, separately identifying lateral fluxes of subsurface water as well as redistribution of snow and surface water, with groundwater flow and related heat advection (Kurylyk et al., 2016;Sjöberg et al., 2016;Grant et al., 2017;Bisht et al., 2018). Some studies assessed the influence of microtopography on the biogeochemical processes in polygonal tundra areas (Cresto Aleina et al., 2013;Kumar et al., 2016;Abolt et al., 2018). These modeling studies have improved our understanding of thermokarst associated processes and feedbacks. The knowledge gained in earlier studies has become an important resource for recent model development that can simulate more realistic icerich permafrost thaw dynamics and successive fine-scale ground subsidence (Aas et al., 2019;Nitzbon et al., 2020). However, the improved models have limitations in terms of representations at regional and global scales. To overcome these limitations, the Community Land Model (CLM) included land subsidence related to permafrost thaw and ground ice melt in a global simulation (Lee et al., 2014). Using a sub-grid tiling approach allows land surface models to represent lateral fluxes of heat, water, and snow in the degrading permafrost terrain (Qiu et al., 2018;Aas et al., 2019;Cai et al., 2020;Nitzbon et al., 2020).
Previous observational and modeling studies have provided evidence of continuously degrading permafrost under climate warming (Liljedahl et al., 2016;Cai et al., 2020). The model projections likely have biases that arise from assumptions about model parameters and limitations in representing the inherent distribution of ground ice relevant to permafrost degradation.
Nonetheless, one of the advantages of incorporating excess ice in the model is that we can capture features of changes in permafrost landscapes and hydrology that result from thawing permafrost and melting excess ice. Therefore, we incorporated a simple framework of excess ice melt and subsidence (Bowling et al., 2008) into a coupled hydrological and biogeochemical model (CHANGE: Park et al., 2011Park et al., , 2018. The objective of this study was to examine the impact of excess ice melt on permafrost hydrology and associated changes in CO 2 and methane (CH 4 ) fluxes under future climate scenarios following the strongest representative concentration pathways (RCP8.5) in a tundra site in northeastern Siberia using the improved CHANGE model. We also quantified the contribution of excess ice meltwater to the surface-subsurface hydrologic system, which is likely achieved because of the inclusion of a water tracer scheme within CHANGE to track pathways of water originating from precipitation and ground ice meltwater (Park et al., 2021).

General Model Description
CHANGE is a process-based land surface model that calculates momentum, water, heat, and carbon fluxes in the atmosphere-vegetation-snow-soil system, soil hydrothermal states, snow hydrology, and stomatal physiology and photosynthesis, in combination with interactions and feedbacks between components and processes within the system (Park et al., 2011(Park et al., , 2018. The model simulates snow accumulation and melt, and the meltwater is partitioned into soil infiltration and surface runoff, even during precipitation events. The infiltrated water can move up or down within the soil, depending on the water potential gradient between soil layers, in which a portion is used for plant transpiration and the excess vertically infiltrated water at either the permafrost table or the bottom soil boundary layer is converted to lateral water flux that is assumed to be subsurface flow. CHANGE calculates the energy balance over the canopy, snow, and soil surface separately and also solves the water budget at both the canopy and soil based on the mass conservation principle. The soil evaporation, transpiration, and rainfall interception by the canopy, which consists of evapotranspiration, are calculated by the balance solution, and their summation is evapotranspiration. The model explicitly represents the dynamics of water and heat flows in permafrost soil considering the freezing/thawing phase transitions in a soil column of 70.4 m depth to explicitly capture thermal inertia in deep ground. The heat flux into the soil-snow is obtained by solving an energy balance equation and is used as the upper boundary condition for the solution of the heat conduction. Zero heat flux is prescribed as the lower boundary condition. The effects of soil organic carbon (SOC) on the soil thermal and hydraulic dynamics are represented by the parameterization of the vertical distribution of SOC within the soil depth of 1.2 m. The calculated ice content in frozen soil layers is used for the parameterization of effective porosity that represents the soil water stress conditions associated with drying and/or freezing, which is coupled to stomatal conductance and the maximum rate of carboxylation, consequently controlling plant productivity as well as phenology. CHANGE mechanistically calculates leaf photosynthesis via the Farquhar biochemical model (Farquhar et al., 1980) as a function of absorbed light, leaf temperature, CO 2 concentration within the leaf, and Rubisco enzyme capacity for photosynthesis. Leaf stomatal conductance, which is needed for the water flux, is coupled to the leaf photosynthesis and concentrations of CO 2 and water vapor (Collatz et al., 1991). The model calculates photosynthesis and stomatal conductance separately for sunlit and shaded leaves. CHANGE simulates biogeochemical processes across multiple biomes, including the carbon and nitrogen cycles in vegetation, litter, and soil. The carbon absorbed by photosynthesis is partitioned into vegetation components, and the coupled nitrogen cycle allows nitrogen limitation to influence plant productivity that is associated with the dynamics of SOC. The decomposition rate of SOC, controlled by soil temperature and moisture, in turn affects productivity through nitrogen availability.
Mosses covering the subarctic soil surface strongly affect water and heat fluxes because of their high water holding capacity and the provision of insulation. CHANGE was coupled with a moss process module, enabling an explicit representation of heat, water, and carbon exchanges in the ecosystem (Park et al., 2018). The moss cover that is assumed to be uniformly distributed over the soil surface is divided into three vertical layers for the calculations of heat and water flows. The temperature of each moss layer is calculated using the heat conduction equation, which is forced by the heat flux into the moss surface from the atmosphere and snow bottom layer, solving the energy budget equation (Launiainen et al., 2015;Park et al., 2018). The dynamics of water content between moss layers are calculated by the water budget equation, and then the fraction of liquid and ice water in individual moss layers is estimated by the simulated temperature. The estimation of thermal properties (i.e., heat conductivity and capacity) in the moss layers is dependent on the separated water and ice fractions. The moss net CO 2 exchange depends on climatic conditions and moss water content and temperature, based on the prognostic moss leaf area index over the annual time series. Therefore, the model assumes that the seasonal cycle of photosynthetic parameters of the moss layers was negligible, and all changes were considered to be instantaneous and reversible.
CHANGE includes a water tracer model that is designed to illustrate the spatiotemporal variability of water originating from precipitation and ground ice in the hydrologic system and quantify their contributions to hydrological processes (Park et al., 2021). It therefore has an ability to assess the contribution ratio of excess ice meltwater to the subsurface flow under future climate warming, which is important to understand the dynamics of hydrological processes.

Excess Ground Ice Model
CHANGE extended the soil column to a depth of 70.4 m, discretizing into 31 soil layers by default. The soil surface (e.g., ∼0.37 m) over which the soil water gradient is generally strong forms six soil layers based on an exponential equation, and then the soil layers have a 0.2-m node spacing to 3.2 m. The deepest remnant layers in turn form thicker node depths based on the exponential equation. The soil layers from the surface to 6 m are treated as hydrologically active with the remaining deeper layers representing bedrock. In the active soil layers, the appearance of permafrost operates the parameterization of ice blocking, thereby constraining the hydrological activity in the relevant soil layers regardless of the model configuration. The maximum ALT is updated during the model calculation period including spinup, and the frozen soil below the active layer which is defined as permafrost soil. Excess ice is initially distributed within the configured soil layers and below the ALT, with the same ice fraction (Table 1), which increases both the volume of the whole soil layer (including saturated soil moisture) and the soil layer thickness accordingly (Bowling et al., 2008;Lee et al., 2014). The increase in the soil layer thickness is linearly proportional to the volumetric content of excess ice. As excess ice melts in response to climatic warming, the soil layer thickness decreases, which corresponds to surface subsidence because of excess ice melt (Cai et al., 2020). The water holding capacity of the soil layer is simultaneously decreased by the melting. The subsidence induced changes in saturated soil moisture and soil layer thickness directly affect hydraulic characteristics such as hydraulic conductivity and saturated soil water suction. These changes influence the vertical fluxes of soil water and ultimately the moisture content of the soil layers, as described in the Supplementary Appendix. Because excess ice is located below the ALT, however, there is little to no change in excess ice and related variables during model spinup.
The excess ice substantially affects soil thermal conductivity and soil heat capacity, which is immediately linked to soil temperature through their influences on heat flux and phase change. Following Lee et al. (2014), CHANGE adopted the thermal parameterization that includes the influence of the excess ice, as described in the Supplementary Appendix. Because the excess ice is prescribed as an initial condition, it only melts and does not grow again. As thawing occurs at a layer, permafrost soil of the layer is preferentially thawed, and thereafter excess ice. These assumptions allow simplification of the model calculation. The excess ice meltwater is added to liquid water of the soil layer or an unsaturated layer above if the soil layer is saturated. The liquid water then participates in the hydrologic system of transpiration and subsurface flow generation in the way of soil water movement.

Methane Flux Model
The coupling of a CH 4 biogeochemistry model integrated in the community land model (CLM4Me, Riley et al., 2011) to CHANGE was made to examine interactions of climate-excess ice-CH 4 fluxes. CLM4Me has been improved by integrating several new process representations and parameterization, and the improved versions have been tested and evaluated at global and site scales (Meng et al., 2012;Melton et al., 2013). Because the detailed descriptions of CLM4Me have been provided in previous studies (Riley et al., 2011), we here describe the major characteristics of CLM4Me, including the major modifications.
In the model calculation, CLM4Me first controls the initialization of boundary conditions, inundation, and impact of redox conditions and then calculates CH 4 production, oxidation, transport through aerenchyma, ebullition, aqueous and gaseous diffusion, and the overall mass balance for unsaturated and saturated soils. Above all, mechanistically modeling net surface CH 4 emissions requires representing a time series of changing inundated fraction. CHANGE can optionally use a parameterization based on TOPMODEL (Beven and Kirkby, 1979) that includes a groundwater model calculating changes in water table controlled by water recharge and lateral flow. As the water table becomes higher, more of the surface area becomes saturated, thereby flooding the regions with higher topographic index first. Here, the saturated fraction (f sat ) calculated as in Gedney and Cox (2003) was alternated with the inundated fraction that is required in CLM4Me. The water table depth is calculated, and the probability distribution function of the topographic index within the site is then used to describe the relative frequency of occurrence of the topographic indices. f sat of the site can be found by integrating the probability distribution function of the topographic index. This is determined by numerical integration because a two-parameter gamma distribution is used for the probability distribution function.
An exponential distribution of f sat is fitted to the parameters of the gamma distribution. f sat is allowed to change at each time step, in which the calculation is made four times in total in the saturated and unsaturated fractions for CH 4 and O 2 . The total soil CH 4 quantity is conserved by evolving CH 4 to the atmosphere when f sat decreases and averaging a portion of the unsaturated concentration into the saturated concentration when f sat increases. In addition, our simulation neglected the solution of the CH 4 and O 2 mass balance for water bodies.

Model Experiments
Several simulations were conducted for the years 1980-2099, based on various modeling setups of excess ice, as summarized in Table 1. The case "CNTR" is the standard version of the model that excludes the process of excess ice. The cases that included excess ice had different settings for the surface depth and content of excess ice along the thickness of 3 m. We combined two types of daily forcing data for running the modelthe forcing data for the current climate from 1980 to 2013 were yielded by the combination of observations at the study site and ERA-Interim reanalysis data (Miyazaki et al., 2015) and the future forcing data from 2014 to 2099 were obtained from a single ensemble member of a projection simulation (RCP8.5) with Hadgem2-ES. Individual experiments were spun up for 1,200 years using the detrended forcing data of the initial 20 years and a CO 2 concentration of 350 ppm. Through the spinup runs, the model determines dynamic equilibriums for soil temperature and moisture and vegetation carbon and nitrogen. A static land cover type of arctic grass was defined for the simulation, while the vegetation phenology was prognostic based on the estimated carbon and nitrogen contents. The RCP8.5 scenario is a high greenhouse gas emissions scenario (Moss et al., 2010), which resulted in strong warming in the study site; for example, air temperature increased 3°C and there was three times higher precipitation in the 2090s when compared with the 1980s ( Figure 1A). It is noted that projections by climate models tend to contain biases in the Arctic, particularly excessive winter snowfall that likely results in unexpected thicker snow depth, as identified by the CNTR experiment ( Figure 1B). The snow depth leads to excessive soil insulation that likely results in rapid permafrost thawing in combination with the warmer air temperature than when the climate biases are removed. To explore the response of permafrost to the precipitation-induced snow depth change of strong emission scenarios, we conducted additional experiments that reduced precipitation by 50% when compared with the default condition (Table 1), which was called "the treated experiment." Although the reduced precipitation produced lower snow depth during the winter than the original experiment ( Figure 1B), the future snow depth of the treated experiment was considerably larger than the maximum snow depth (i.e., 40 cm) observed under the present climate.

Site Description and Observations
The Tiksi field site is located in northeastern Siberia (71.4°N, 128.5°E), near the mouth of the Lena River and 5 km west of the Laptev Sea coast (Figure 2). The study site receives about 173 mm of precipitation annually, with 60% in the form of rain during the summer (June to September). Annual mean air temperature is about −13.5°C with an annual amplitude from an average of −33.3°C in January to an average of 7.5°C in August. Strong wind speed results in usually heterogeneous and shallow snow cover, recording maximum snow depths of about 40 cm. The site represents a typical lowland tundra landscape, characterized by small lakes and low flat plains of wet tundra, and hill slopes and mountains of dry tundra. Various vegetation types, such as non-tussock sedge, dwarfshrubs, and moss, are distributed on flat plains and the lower parts of relatively gentle slopes. The top and steep-sloping part of the mountains is covered by gravel with lichen. Continuous and ice-rich permafrost underlies the ground. The thermoerosion induced by the climate warming increases the expansion of thermokarst lakes and small ponds (Liljedahl et al., 2016). The tundra soil consists of water-/ice-saturated sandy peat, with the water table close to the surface. The mineral soil is a clay loam consisting of 35% silt, 32% sand, and 33% clay. The quantity of organic matter around the site is significantly heterogeneous, with values ranging between 3 and 85 kg m −3 , and mean SOC values of 35 kg m −3 in the soil layer above 50 cm.
Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 704447 FIGURE 1 | Variations in annual precipitation (blue) and annual mean air temperature (red) (A) and snow depth (B) simulated by the original forcing climate (default experiment) and when precipitation was reduced by 50% when compared with the original data (treated experiment). Air temperature and precipitation until 2013 were yielded by the combination of in situ observations and ERA-Interim reanalysis data (Miyazaki et al., 2015), and thereafter those were obtained from a single ensemble member of a projection simulation (RCP8.5) with Hadgem2-ES (Collins et al., 2008), as described in Model Experiments section. Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 704447 5 Various meteorological variables have been recorded from a 10-m high observational tower since 1997. The meteorological variables (i.e., shortwave and longwave radiation, air temperature, precipitation, relative humidity, wind speed, and air pressure) were measured at different heights on the tower; however, heat, water vapor, and CO 2 and CH 4 fluxes have not been observed directly at this site.
The heterogeneous surface and subsurface properties of the study site likely produce large spatial variations in the soil thermal regime. ALT was monitored at a 1 km × 1 km CALM site established in 1998 (Iijima et al., 2016). The measurement was performed in grids with 100-m intervals, once a year in mid-September using a metal rod. Another site of mosscovered ground consisting of 12 points was established near the CALM site, to further monitor ALT in 2006. Two observation plots monitoring soil temperature were established in 2004the first plot was located in a flat and very wet area and the second was in a flattened wet area and a narrow valley about 500 m from the first site. The landscape of the plots was polygonal tundra, with low shrub, sedge, and moss cover. Moss cover was dominated by sphagnum. Each plot was equipped with two TR-52 single-channel temperature loggers (T&D Corp., Matsumoto, Japan). The temperature sensors were placed on the soil surface and in permafrost at depths of 0.5 and 0.6 m.

Model Performance
The CHANGE-simulated results for water and carbon fluxes and soil moisture and temperature profiles, which have previously been validated to observational records in this study site, showed statistically satisfactory performances for seasonal and interannual variability (Park et al., 2018). Although the CH 4 process was newly added to CHANGE in this study, deficiencies in the observational data constrain the validation of model performance. Here, the ALTs simulated under various model settings were compared against the observational records, which enabled us to address the influence of climatic and surface conditions on ALT. The comparison between the simulated ALT and the observation shows that CHANGE captured the interannual variability of the observed ALT with +6 cm deviation during the period of 1997-2013 simulated by the current forcing climate (Figure 3). The observational values were averaged by records of points with different microclimates, soil properties, and vegetation characteristics. The onedimensional CHANGE cannot accurately simulate the influences of these heterogeneities. This scale problem is one plausible reason for the difference in ALT between the simulated and observed values (Park et al., 2018). ALT was influenced more by snow cover than by moss cover (Figure 3). For example, the ALT simulated by EXCL MOS_BLS that excluded blowing snow was approximately 30 cm deeper than that by EXCL MOS that considered blowing snow, indicating the large cooling effect of thinner snow depth caused by blowing snow on ALT. The effect of moss cover on ALT was not larger than that of the blowing snow; CNTR, which considered moss cover, resulted in 5 cm lower ALT than EXCL MOS . Meanwhile, the comparison between CNTR and EXICE T1_C15 showed small differences in ALT, revealing the slight impact of excess ice on ALT. Because the excess ice was located 1 m deeper than the ALT, it was part of the permafrost. Therefore, the latent heat effect was dominant at the ALT bottom (i.e., the permafrost table). Instead, the excess ice has a large cooling effect on permafrost temperature. The differences in ALT between the experiments (Figure 3) suggest that the representation of moss cover and blowing snow is an inevitable consideration in models for the correct projection of permafrost dynamics in tundra regions influenced by climate change. Previous studies highlighted the critical impact of SOC on ALT. For example, a simulation for Samoylov Island of northern Siberia, close to this study site, addressed that the inclusion of SOC reduced ALT by 50 cm relative to the exclusion (Chadburn et al., 2015).

Active Layer Thickness and Subsidence
The simulated ALT exhibited evident deepening during 1980-2099, driven by the warming climate ( Figure 4C). The deepening was initiated from 2014 forced by the RCP8.5 climate conditions and became mostly significant from 2040 onwards when excess ice began to melt ( Figure 4B). The rate of deepening increased over time, particularly under conditions of moss-free cover and excess ice exclusion ( Figure 4C), where the resultant higher heat fluxes into the soil column were effective for the warming of permafrost. EXCL MOS_BLS of the default experiment resulted in 12 m of ALT during 2080-2099, which was considerably larger than that of the other simulations ( Table 2). Comparison of the experimental ALTs confirmed that the presence of excess ice limited the ALT ( Table 2). The limitation was largest in the cases with 30% ice contents of 3 m table; for example, ALT of EXICE T3_C30 was 4 m during 2080-2099, which was three times lower than that of EXCL MOS_BLS because of the larger latent heat requirements for the excess ice melt. The surface soil layers repeat seasonally thawing and freezing even under the strong future emission scenarios, and the warming climates likely develop annually unfrozen talik layers beneath the ALT. It is therefore defined that ALT in this study is the maximum depth of unfrozen soil including talik.
The treated experiments that reduced precipitation by 50% showed a time series of ALT variability similar to that of the default experiments. The absolute values of ALT under the same conditions were 45% larger on average in the default experiments than in the treated ones during 2080-2099 (Table 2), indicating the impact of snow on the ALT caused by the increased soil temperature. In the case of EXICE T1_C15 , for example, the snow depths in the two future periods (2030-2049 and 2080-2099) were considerably higher than in the present climate, despite later snow accumulation (Figures 5A,C). The thicker snow depth evidently has a larger insulation effect, preserving permafrost warming formed during summer and resulting in deeper ALT. A modeling study that treated snowfall magnitudes and snow seasons differently identified that changes in snow depth had larger impacts on thermal conditions of continuous rather than discontinuous permafrost and on colder Siberia than warmer North America (Park et al., 2015). The default experiment resulted in ALT of approximately 8 m on average during 2080-2099. The deeper ALT was also influenced by talik that was formed at the soil depth of 1 m or lower ( Figures 5B,D). A large volume of water stored in the talik soil layers increases heat capacity and thermal conductivity, effectively dampening the annual temperature cycle (Subin et al., 2013;Lee et al., 2014), as identified in the seasonal ALT of the default experiment ( Figure 5B). Interestingly, the talik was also found for the treated experiment with relatively thin snow depth ( Figure 5D), which highlights the strong influence of the warming temperature on the talik, with high insulation effect of the snow depth and moss cover. The different setups of simulations in both experiments showed that excess ice in this study site preserved the initial  Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 704447 state until 2040 because of the combined effects of less extreme emission scenario, the insulation of moss cover, and higher latent heat of excess ice, and thereafter it mostly melted by the end of the 21st century ( Figures 4A,B). However, the melting time series showed evident differences between the simulations ( Figure 4B). The timing and magnitude of excess ice melt was influenced by the distribution and amount of the initial excess ice; for example, the case EXICE T1_C15 with lower excess ice resulted in advanced excess ice melt by approximately 10 years relative to the case of EXICE T1_C30 ( Figure 4B). In the case of EXICE T3_C30 , interestingly, a portion of excess ice remained until 2099 ( Figure 4B) because of the higher latent heat effect associated with the higher excess ice content distributed at the deeper soil layer. The excess ice melt resulted in land surface subsidence. The estimated annual subsidence ranged from 2 to 40 cm in the default simulations ( Figure 4A). A simulation for Samoylov Island of northern Siberia, close to this study site, estimated surface subsidence of 35 cm by around 2030 and more than 1 m by the end of this century, based on the RCP4.5 warming scenario (Aas et al., 2019). Our results for subsidence show discontinuous occurrence in the time series of individual simulations. As permafrost thaws under the warming climate, excess ice simultaneously melts in nature. However, CHANGE used a physics model where excess ice starts to melt after permafrost of a soil layer is completely thawed. As the permafrost layer becomes active, the soil layer has already stored enough heat to rapidly melt the excess ice. Therefore, the physics model resulted in stepwise patterns of subsidence.

Water Flows
Permafrost and excess ice that are located at the base of the active layer greatly impede water flow below the layer. This leads to a relatively wet active layer and high water level. The simulations exhibited high water levels until 2040 ( Figure 6A) when the excess ice was still intact (Figure 4); the water level rapidly lowered after the excess ice melted ( Figure 6A). In the physics model, excess ice meltwater is added to soil moisture. The meltwater causes the base of the active layer to wet, so that the model groundwater becomes more active, leading to higher subsurface flow as observed in the simulations ( Figure 6B). The runoff abruptly increased after 2040, although there were differences in the amount of runoff between the simulations ( Figure 6B). In the treated experiments, the annual subsurface flow averaged 141.4 mm during 2080-2099, which corresponded to 31% of the annual precipitation averaged during the same period.
As the ALT thickens, the meltwater from both excess ice and permafrost increases. Substantial amounts of unfrozen liquid water, defined as older soil moisture in the tracer scheme, is stored in permafrost soil. As the permafrost thaws, the liquid water becomes active and is involved in talik formation and subsurface flow. The tracer model quantified the largest contribution of the old soil moisture, which was the initial unfrozen water in the model setting, to the subsurface flow during 2080-2099 ( Figure 7D). The excess ice meltwater accounted for 26% of the subsurface flow during the same period, which was not identified in the present climate. Both the old soil moisture and excess ice meltwater maintained high contribution rates to the subsurface flow with changing seasons, strongly implicating the talik in the runoff processes. The influence of the snow meltwater and rainwater on subsurface flow considerably decreased in the future when compared with the high fractions in the current climate ( Figures 7B,D). The simulated results also displayed low fractions of the snow meltwater and rainwater constituting deeper soil moisture in the future, particularly at 1.6 m or deeper, while the excess ice meltwater increased at the same soil depths relative to the current climate ( Figure 8).
The future warming climate resulted in a drier environment in the upper soil layers relative to the present (Figure 8), despite the larger precipitation ( Figure 1A). As ALT thickens, the higher conversion of excess ice meltwater to subsurface flow, as described above, was related to soil drying. The dried soil layers need more water for soil wetting. A portion of the available precipitation is used for soil wetting, which ultimately reduces water flowing to the deeper soil layers. In terms of soil moisture distribution, evapotranspiration increases moisture loss from the upper soil layers, thereby resulting in increased soil drying. The dry soil resulted in the decline of the water table ( Figure 6A). The water table exhibited the declining trend toward the end of the 21st century, in which the interannual variability was also similar to the patterns of the subsurface flow and ALT, representing their interaction.

Carbon Dioxide and Methane Fluxes
Soil moisture in the upper soil layers (i.e., 0-1 m) during 2080-2099 mainly distributed within the range of 0.2-0.3 m 3 m −3 , although it was drier than during 1980-1999 ( Figure 8A). The future soil moisture was never low enough to constrain vegetation productivity or SOC decomposition. The simulated net primary productivity (NPP) of the study site continuously increased during the study period ( Figure 9A), with slight differences in values between the model experiments during 2080-2099 ( Table 2). The increase in NPP was large from 2014 when it was forced by the warming climate of RCP8.5. In the warming climate, there was no large difference in the annual mean NPPs of the two experiments. Interestingly, however, NPP of the default experiment showed slightly lower values than the treated one from 2060 onwards, despite the larger precipitation in the default experiment. The deeper ALT in the default relative to the treated experiment was linked to the decrease in soil moisture, as demonstrated by the lower water table (Table 2), which likely somewhat limits ecohydrological functions. The high precipitation during the summer season wets a large fraction of the vegetation, limiting temporary photosynthesis by plants  (1980-1999, left) and future (2080-2099, right) climates in the case of EXICE T1_C15 of the treated experiment. The contributions of source waters to runoffs were separated by the tracer scheme in CHANGE. The old water in (D) indicates the water that was existent as unfrozen type in the soil column by the initial model setup.
Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 704447 (Park and Hattori, 2004). The simulated heterotrophic respiration (HR) exhibited similar variations to those identified for NPP ( Figure 9B). The simulated CH 4 efflux displayed complex interannual variability and decadal trends ( Figure 9C), considerably different to those identified for the simulated NPP and HR ( Figures 9A,B). The warming climate forced increasing CH 4 emission and HR from the permafrost soil to the atmosphere until 2040 when permafrost was still stable (Figure 4C), and the water level was relatively high ( Figure 6A). Then, CH 4 effluxes generally tended to decrease, although the timing was quite different according to the simulations; for example, CNTR of the default experiment represented the earliest decrease of CH 4 efflux by 2040, because of the decreased soil moisture induced by the early permafrost thawing, while the experiments with excess ice simulated higher CH 4 effluxes until the mid-2070s. The excess ice meltwater temporarily creates a wetted soil that efficiently generates CH 4 effluxes. The model effectively simulated higher CH 4 effluxes during excess ice melt. In contrast, the high conversion of excess ice meltwater to subsurface flow induced soil drying, as explained above, resulting in an abrupt decrease in CH 4 effluxes because of the aerobic conditions. Despite of the soil drying, both NPP and HR continuously increased responding to the warming climate, inconsistent with the decreasing CH 4 efflux. The precipitation in the future climate ( Figure 1A) was as large as to alleviate the drying of the surface soil layers ( Figure 8B), which contributed to higher vegetation productivity and respiration. The same variations in CH 4 effluxes were identified in the treated experiments as those in the default experiments, although the timing of the projected higher CH 4 effluxes was later relative to the default experiments. These model experimental results suggest that excess ice is an important component in the complex interannual variability of CH 4 efflux.

Sensitivity Analysis
The excess ice constrained the deepening of ALT; the magnitude of this response fluctuated depending on the ice content ( Figure 4B). For example, at the cases that the surface of Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 704447 excess ice was 1 m, EXICE T1_C30 with excess ice content of 30% lowered ALT by 50 cm relative to EXICE T1_C15 with 15% during 2080-2099 (Table 2). At the cases of 3 m, the difference in ALT between EXICE T3_C15 and EXICE T3_C30 increased ALT by 2.7 m during the same period ( Table 2), indicating larger latent heat effect by higher excess ice content. The permafrost including excess ice impedes water flowing downward, wetting soil layers above the permafrost table and thereby resulting in a higher water level ( Figure 6A). Under those conditions, the warming temperature warms the active soil layers, consequently increasing SOC decomposition. Excess ice melt contributed to increased CH 4 efflux, with a larger response obtained for high ice content than for low ice content ( Table 2). The depth of the excess ice table was also closely associated with the CH 4 efflux, which was lower for an ice table depth of 3 m than of 1 m ( Table 2). As the table of excess ice deepens, the thicker soil column intuitively slows the soil wetting and warming. These conditions are effective for limiting the CH 4 efflux. Even though soil wetting and warming are extended to the deeper layers, CH 4 efflux under the changed conditions might not be large because of smaller SOC storage in the deeper soil layers. The maximum estimated ALT was 60 cm during the model spinup of 1,200 years. SOC is mainly stored in the active layers. As ALT is thickened by the warming climate, SOC also flows down. However, the amount of SOC flowing to 3 m was low, hence the low CH 4 efflux.

DISCUSSION
Because of the lack of excess ice datasets and observational evidence, our model projections of excess ice melt and associated processes likely have biases that arise from model settings and simplifications of the excess ice initialization. However, the sensitivity experiments can potentially improve our understanding of the impacts of excess ice melt on hydrological and biogeochemical processes. The water tracer Frontiers in Earth Science | www.frontiersin.org September 2021 | Volume 9 | Article 704447 model in CHANGE tracks pathways of excess ice meltwater in the soil system and assesses their contribution to hydrological processes. The excess ice and permafrost act as a barrier to water infiltration because of the low conductivity of icy soils, storing the inflowing precipitation water within the active layers. The large volume of stored water induces higher heat capacity and conductivity, dampening the annual soil temperature cycles (Subin et al., 2013), in which warming climate forces permafrost thawing and excess ice melt. As soon as the excess ice melts, the ALT bottom layers are wetted by the melt water. The wetted soil moisture was connected to the subsurface flow (Liljedahl et al., 2016). Once excess ice melts, there is no additional water production from the soil layer. During these processes, subsidence occurs, and the excess ice meltwater is linked to the upper soil moisture that is mainly the result of precipitation (i.e., snowmelt and rainwater). However, evapotranspiration and the recharge of soil moisture preferentially use precipitation water inflowing from the surface, and thus there is limited connectivity of the precipitation to the deeper soil moisture. As a result, the soil column tends to be dried as ALT deepens, and the water level is lowered. CHANGE captured both the weakening signal of the connectivity of surface soil moisture to the deeper soil, and the hydrological regime of thawed permafrost tipping as a consequence of melted excess ice to a dominance of unsaturated soil condition (Aas et al., 2019;Nitzbon et al., 2020). The excess ice has a contrary effect of temporarily limiting rapid soil drying in the melting stage. At the bottom of the active layer, the excess ice melt water-induced wetting reduces water input from the upper soil layers because of the reversed large moisture gradient. The moisture gradient gradually decreases with the increasing subsurface flow, so that soil is dried. The excess ice melting lasted for 15-20 years dependent on the excess ice settings. In practice, the model setups that excluded excess ice simulated earlier soil drying than experiments that included the excess ice ( Figure 6A).
The simulated regime of soil moisture is strongly dependent on the assumptions for the modeling. Under water-logged conditions, the melting of excess ice would develop deep water-filled troughs, and potentially form thermokarst lakes and ponds (Nitzbon et al., 2020). Furthermore, excess ice melt and associated ground subsidence increases the lateral hydrological connectivity of the natural landscape, causing higher water levels in the depressed lowlands. However, models lacking in the representation of thermokarst-inducing processes likely underestimate permafrost thaw dynamics and resultant carbon fluxes. A modeling study highlighted differences of 2-12 times in SOC decomposition between two experiments that included and excluded thermokarst-inducing processes in northeastern Siberia (Nitzbon et al., 2020). As mentioned before, CHANGE employed one-dimensional physics of permafrost thaw dynamics considering only the vertical heterogeneity of excess ice. It was assumed that the excess ice meltwater was connected to soil moisture dynamics and subsurface flow generation. The simplistic representation of CHANGE likely yields biased projections for the subsurface hydrology crucially influenced by heterogeneous microtopography, lateral hydrological connectivity and ground subsidence. In reality, CHANGE simulated wetter soil and higher subsurface flow after the melt of excess ice (Figure 7) and subsequent higher carbon releases (Figure 9). However, our simulations provide important insights for excess ice melt and associated changing hydrology and carbon fluxes under well-drained conditions. Following the simulated results, the default experiments with large precipitation produced large HR and CH 4 (Figure 9) until 2040 when permafrost was stable, while larger production of HR and CH 4 was found in the treated experiments after this time. However, the differences between the two experiments were not large. This is probably because of the relatively similar moisture environments in their upper soil layers that were relatively sensitive to precipitation. Although the default experiments simulated dry soil and lower water table level (Figure 6), the high precipitation was efficient in alleviating drying of the upper soil layers.
The CH 4 fluxes showed different variations than that of HR. After 2060 when excess ice melted, the interannual variability of CH 4 effluxes was large and exhibited considerable differences between simulations (Figure 9). Earlier occurrence of high CH 4 efflux was found in the simulations with the surface of excess ice of 1 m than in ones with 3 m. Moreover, CH 4 efflux was further increased in the simulations with higher volumetric contents of excess ice. These results indicate that the initialization of excess ice greatly affects CH 4 efflux. Therefore, the model projection of CH 4 efflux utterly relies on the initial settings of excess ice. Because of the scarcity of observational data, however, this study employed the active layer-dependent excess ice initialization for modeling, together with statistical and empirical knowledge of the vertical extent and depth of excess ice as used in previous modeling studies (Lee et al., 2014;Aas et al., 2019;Cai et al., 2020). These configurations of excess ice still likely deviated from the actual ice features. We do not therefore expect the modeled excess ice melt and CH 4 efflux in this study to be an adequate representation of reality. However, our simulations could be used to infer the timing and magnitude of excess ice melt and associated carbon and CH 4 fluxes, and possible changes in soil hydrology under the well-drained condition.
Locally developing thermokarst lakes have recently been observed in this study site (Liljedahl et al., 2016), which could have been induced by natural or anthropogenic disturbances like the accumulation of snow in topographic depressions, the removal of vegetation or moss cover, and forest and tundra fires (Jones et al., 2015;Nauta et al., 2015;Abolt et al., 2018). Here, our model results provide evidence for the early excess ice melt as a result of higher insulation of thicker snow and thus enhanced soil drying. The simulation that excluded moss cover also estimated considerably earlier excess ice melt than the undisturbed case. In contrast, the high cooling by moss cover limited ALT deepening, hence stable permafrost. A simulation including excess ice initialization for the Lena River delta close to our site estimated the initiation of excess ice melt before 2020 (Cai et al., 2020). The melt timing is quite fast relative to our result showing stable permafrost until 2040 (Figure 4). The different model structure and settings make it impossible to directly compare the two results. Considering that the model used by Cai et al. (2020) did not consider the influence of moss cover, here we could expect further easy warming of permafrost under the moss-free condition. Park et al. (2018) found that 4-cm-thick moss cover reduced the ALT by 14 cm relative to the moss-free case. These results suggest that model improvement including moss cover as a land component is required for realistic representation of permafrost features under warming climate conditions.
Recent observations reported complex landscape development in the terrestrial Arctic, such as thermokarst ponds and lakes of meter-scale distribution (Liljedahl et al., 2016;Ulrich et al., 2017).
To represent the small-scale permafrost features, models employed sub-tile interactive setups that simulated microtopographic changes associated with degradation of icerich permafrost, considering the lateral water flux (Aas et al., 2019;Nitzbon et al., 2020). The sub-tile representation of excess ground ice within permafrost soil was used for the global simulation without the lateral flux (Cai et al., 2020). The subtiling approach preferentially needs available observational data for the model simulation, which is substantially consistent with the requirements of our model. Therefore, the introduction of the tiling system could be a direction for the future development of CHANGE, particularly by configuring an optional function whereby lateral heat and water fluxes could be switched on and off dependent on topography. Laterally coupled tiling likely builds more potential for better modeling subscale variability of soil and snow, consequently improving permafrost-water-carbon feedbacks.

CONCLUSION
This study coupled an excess ice scheme to CHANGE and examined the impacts of excess ice on permafrost, soil water dynamics, and carbon and methane fluxes in a Siberian tundra site under the strong future emission scenarios. The model results indicated that the warming air temperature and higher snow depth were major factors deriving ALT increasing and the resultant permafrost degradation. Under the warming permafrost, in contrast, excess ice thermally stabilized permafrost until it melts, which was further significant under the conditions of larger excess ice content and moss cover. As excess ice melts, however, the meltwater was connected to hydrologic and thermal regimes, and the resultant wetted and warm soil was positively fed back to higher subsurface flow, permafrost degradation, and temporarily higher CH 4 efflux. These results provide insights for interactions and feedbacks between climate change, permafrost excess ice, and carbon fluxes in a well-drained condition.
The insights obtained by this study represented that the magnitude and timing of excess ice melt are largely dependent on the model settings. The model parameterization and initialization of excess ice remain necessarily simplistic, which highlights the importance of the connectivity of the excess ice meltwater to the subsurface flow and subsequent implications for carbon fluxes. Therefore, future model improvement emphasizes the necessity of coupling tiling hierarchy to CHANGE, which would enable a more comprehensive examination of the possible trajectory of the permafrost carbon-climate feedback under changing climate.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
HP developed the ideas, wrote most of this paper and drew most of the figures. All authors participated in data processing and preliminary analysis; HP and TH coordinated the model experiments and analyzed the simulation results; AF and PK collected and analyzed in situ observational data.

FUNDING
This study was partly supported by the Japan Society for the Promotion of Science KAKENHI Grant Number 17H01870, 19H05668, and 21H04934. It was also supported by JST Belmont Forum Grant Number JPMJBF2003, Japan, and the Arctic Challenge for Sustainability II (ArCS II; JPMXD1420318865) project funded by the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan.