Exploring the interplay between soil thermal and hydrological changes and their impact on carbon fluxes in permafrost ecosystems

Accelerated warming of the Arctic can affect the global climate system by thawing permafrost and exposing organic carbon in soils to decompose and release greenhouse gases into the atmosphere. We used a process-based biosphere model (DVM-DOS-TEM) designed to simulate biophysical and biogeochemical interactions between the soil, vegetation, and atmosphere. We varied soil and environmental parameters to assess the impact on cryohydrological and biogeochemical outputs in the model. We analyzed the responses of ecosystem carbon balances to permafrost thaw by running site-level simulations at two long-term tundra ecological monitoring sites in Alaska: Eight Mile Lake (EML) and Imnavait Creek Watershed (IMN), which are characterized by similar tussock tundra vegetation but differing soil drainage conditions and climate. Model outputs showed agreement with field observations at both sites for soil physical properties and ecosystem CO2 fluxes. Model simulations of Net Ecosystem Exchange (NEE) showed an overestimation during the frozen season (higher CO2 emissions) at EML with a mean NEE of 26.98 ± 4.83 gC/m2/month compared to observational mean of 22.01 ± 5.67 gC/m2/month, and during the fall months at IMN, with a modeled mean of 19.21 ± 7.49 gC/m2/month compared to observation mean of 11.9 ± 4.45 gC/m2/month. Our results underscore the importance of representing the impact of soil drainage conditions on the thawing of permafrost soils, particularly poorly drained soils, which will drive the magnitude of carbon released at sites across the high-latitude tundra. These findings can help improve predictions of net carbon releases from thawing permafrost, ultimately contributing to a better understanding of the impact of Arctic warming on the global climate system.


Introduction
The stability of permafrost is the result of complex interactions between global and local environmental factors that influence soil thermal and hydrological regimes.The Arctic amplification of climate warming is recognized as one of the main drivers of permafrost thaw observed across the Arctic [1][2][3].As a result, roughly 7% of permafrost has thawed over the last 40 years [4], and models project widespread continued permafrost thaw during the 21st century (e.g.[5,6]).Despite efforts to incorporate in situ observations to assess permafrost carbon dynamics across the Arctic, harsh environmental conditions and the remoteness of many sites have made it difficult to collect high-quality data.This has resulted in limited site-level data and has made upscaling large ecosystem models challenging.Therefore, there is a pressing need to incorporate more field data to improve the performance of these models and reduce their uncertainties through improved parameterization.By doing so, we can gain a more accurate understanding of permafrost carbon dynamics across the Arctic [7][8][9].Process-based biosphere models show that carbon loss from permafrost thaw represents a significant positive feedback to global climate change (e.g.[10][11][12][13][14][15] However, large spatial heterogeneity in the rates and patterns of permafrost thaw and the lack of uniform responses reveal the importance of other environmental drivers such as snow cover, soil composition, organic-layer thickness, vegetation structure, and ground ice content on permafrost stability [2,16,17], in addition to increases in air temperature.The interplay between soil texture, drainage, and moisture, as well as the impact of moisture on thermal conductivity, collectively have significant implications for permafrost dynamics and stability [18,19]. Additionally, the linkages between soil thermal and hydrological regimes exert critical control over the response of permafrost to a changing climate.Hydrology is a fundamental determinant of permafrost stability which modulates the effect of geomorphology on permafrost cryostratigraphy [20][21][22].On the other hand, permafrost influences soil permeability and water flow path, confining it mainly to the seasonally unfrozen parts of the soil column (i.e.active layer [19,23,24].As permafrost thaws, subsurface hydraulic properties are altered, allowing water to infiltrate and move to deeper soil horizons, accelerating permafrost thaw and talik development [25][26][27][28].The fate of soil carbon as permafrost thaws depends on thaw-induced transitions in soil hydrologic conditions and their implications on decomposition rates, the proportion of soil aerobic to anaerobic respiration, and plant, water, and nutrient acquisition.Permafrost thaw in uplands can result in drainage and drying that enhances aerobic respiration (e.g.[29][30][31]) and carbon loss from the decomposition of deep/old organic matter (e.g.[32][33][34]).The increased hydrological connectivity resulting from permafrost thaw may enhance water flow and runoff and increase the mobilization and lateral transport of dissolved carbon [35][36][37].In contrast, water saturation resulting from ground subsidence following thaw of ice-rich permafrost in lowlands enhances anaerobic processes and methane production and emissions [37,38], which may not be fully offset by an increase in post-thaw vegetation productivity (e.g.[39][40][41][42][43]).
Substantial uncertainties remain associated with model projections of the permafrost carbon-climate feedback related to soil thermal and hydrological regimes on the ecosystem carbon balance.The uncertainty associated with ecosystem carbon balance in high-latitude ecosystems is the result of the direct effects of climate and vegetation characteristics on ecosystem productivity and respiration [3,6,12,42,44,45].The main objective of this study is to assess the capacity of a terrestrial ecosystem model to represent the response of soil thermal and hydrological regimes to environmental changes and their consequences on ecosystem carbon balance at two sites differing in soil drainage conditions related to contrasting stages of permafrost degradation.We synthesized observations from two Alaskan tundra sites to inform a model sensitivity analysis and model-data comparison and evaluate the capacity of a terrestrial ecosystem model to represent (1) how changes in soil thermal and hydrological regimes affect both belowground and aboveground carbon dynamics and (2) how the linkages between permafrost, hydrology, and ecosystem carbon dynamics are modulated by climate, drainage conditions, and soil properties.

Sites
We used data from two long-term ecological monitoring tundra sites in Alaska: Eight Mile Lake (EML) and Imnavait Creek Watershed (IMN) shown in figure 1. EML is located on the northern foothills of the Alaska Range (∼670 m elevation) and the home of the Carbon in Permafrost Experimental Heating Research (CiPEHR, 63 • 52 ′ 59 ′′ N, 149 • 13 ′ 32 ′′ W) project [32], which was initiated in 2008.A series of sites were established in tussock tundra, combining controls and treatments that warm air and soil inducing substantial permafrost thaw [46].Landscapescale CO 2 fluxes have been monitored across a permafrost thaw gradient since June 2004 [31] using a combination of clear chambers and the eddy covariance approach (EC).Ice-rich permafrost thaw resulting from the warming treatment resulted in increased waterlogging [41].
Imnavait Creek site is located in the northern foothills of the Brooks Range, in the Arctic region of Alaska (68 • 37 ′ N, 149 • 18 ′ W, ∼760 m elevation).IMN has three eddy covariance towers positioned about 500 m apart, spanning across a topographic gradient with three distinct tundra types, for which landscape-scale EC CO 2 fluxes have been monitored since 2008 [48].The first tower's footprint lies within heath tundra, the second tower characterizes moist acidic tussock tundra, and, the third tower is characterized by wet sedge tundra, where the soil is poorly drained.Additional detailed descriptions of each study site can be found in the supplementary materials.
The EML and IMN sites differ in several aspects.EML is classified as an alpine tundra and is located within the discontinuous permafrost zone, while IMN is an Arctic tundra site located in the continuous permafrost zone.For this study, we set initial soil drainage conditions, to accurately represent current soil and hydrological conditions at each site.We initialized our model runs at IMN specifying a welldrained condition type, while EML was set to poorlydrained [49], to represent soil drainage conditions from subsidence at the site, limiting lateral drainage.The Imnavait site is located within a region of continuous permafrost and is characterized as a well-drained site, while (c) Eight Mile Lake is located within discontinuous permafrost and is characterized as a moderately-drained site.Reported Net Ecosystem Exchange totals (b), (c) are from site observation flux data where EML is shown to be a net source, while IMN is relatively neutral, neither sink nor source.Permafrost distribution source: [47], Imnavait Site Photo source: AmeriFlux.Adapted with permission from Greg Friske.(Main map layer) Adapted with permission from [47].Adapted with permission from Sebastian Biraud.

Model description
The Terrestrial Ecosystem Model with Dynamic Vegetation and Dynamic Organic Soil Layers (DVM-DOS-TEM) is a process-based biosphere model designed to simulate biophysical and biogeochemical processes between the soil, vegetation, and the atmosphere [50].DVM-DOS-TEM v0.6.1 is the latest version of the Terrestrial Ecosystem Model [51,52] that focuses on representing carbon and nitrogen cycles in high latitude ecosystems and how they are affected at monthly to centennial scales by climate, disturbances and biophysical processes such as soil thermal and hydrological dynamics [49,53], snow cover [54] and vegetation dynamic [55].DVM-DOS-TEM has been used extensively in arctic and boreal permafrost regions to evaluate soil moisture, permafrost distribution dynamics [49,53,56], tundra vegetation dynamics [55,57], regional carbon storage and fluxes [58,59], in response to wildfire disturbances [60][61][62] and climate change [63].Additional description of DVM-DOS-TEM can be found in the supplementary materials in addition to equations representing the linkages between soil thermal and hydrological processes equations (S1)-(S4).

Input and study site observations
Historical climate data  for use in DVM-DOS-TEM were extracted from the Climatic Research Unit time series version 4.0 (CRU TS4.0; [64]) downscaled at 1 km × 1 km spatial resolution on a monthly time step, prepared by the Scenario Network for Alaska Planning (https://uaf-snap.org/,[65]).Site observations from each study site were aggregated to a monthly resolution to match the temporal resolution of model outputs.
Observations for EML were compiled mainly from the database of the Bonanza Creek Long Term Ecological Research (LTER) Program, while data for IMN, is archived and can be downloaded from the AmeriFlux database and the Arctic LTER [66].Additional flux and ancillary observations for each site were collected from the Arctic-Boreal CO 2 Flux Database, which has compiled observations provided by the site principal investigators (ABCflux; [67]).A full description of the variables included can be found in (table S1).The impact of permafrost hydrology on methane emissions was not included in the present analysis, as DVM-DOS-TEM does not explicitly represent methanogenesis yet.Interannual variability of model and observational data are presented with standard deviations.

Parameter sensitivity and data analysis
We conducted a parameter sensitivity analysis (PSA) on 14 parameters to identify the most influential parameters affecting subsurface thermal and hydrological conditions and ecosystem carbon cycling.These parameters included soil porosity, bulk density, soil thermal and hydrological conductivities of the fibric and humic organic horizons, and the mineral layers.For detailed insights on snow and n-factor influence on soil conditions, refer to the supplementary materials (equations (S5)-(S8)).Thawed and frozen n-factor (n thawed , n frozen ), which is the ratio of degreeday sum temperature at the soil surface to that in the air [67,68], is used in DVM-DOS-TEM to estimate ground surface temperature from air temperature [49].For the full list of parameters used within the study and their descriptions refer to table 1.
The initial ranges of parameter value were set to ±15% based on the initial PSA ranges, in-situ observations, and literature review, and sampled using a uniform distribution, to better assess the model response to a multitude of possible combinations of parameters within its given range.The results of the PSA were used to optimize parameter values and minimize the RMSE between model outputs and observations [69][70][71].After running the PSA, we used a random forest machine learning method, the RandomForestRegressor in Python's sklearn package [72], to identify the most influential ecosystem parameters on model outputs, which were compared to site observations at a monthly time step for all periods of overlap.

Modeled ecological relationships
We analyzed first-order correlations between climate input, modeled subsurface physical outputs, and biogeochemical outputs for thawed and frozen seasons to better understand the interconnections between physical and biogeochemical processes (see supplementary materials).For model correlation analysis, we defined the thawed season as the period of year when monthly mean air temperatures were ≥0 • C, and the frozen season when monthly mean air temperatures were <0 • C. For both sites at 10 cm, model outputs consistently underestimated temperatures during the annual thawed season; however, the insulating effect of the soil was dampened at the deeper depth, 40 cm which led to a better match between the model and observation (figures 2(g) and (h)).Our model simulations of soil temperature at 10 cm depth at IMN show a slight warming trend, especially over the last 3 years of observations (2013-2016), which is consistent with borehole measurements that observed a warming trend, particularly during the frozen season.These measurements show a significant increase in the 50 cm borehole's minimum winter temperatures, from −9 • C in 2012 to −3 • C in 2015 [48].

3
While DVM-DOS-TEM was able to match the overall seasonal pattern of soil moisture at 5 cm depth at IMN (figure 2(i)), precipitation inputs were significantly lower than estimates from the weather station between 2005 and 2011 (figure S1).These discrepancies in precipitation between the gridded inputs and in situ data also contributed to discrepancies in modeled snow depth.At IMN, the model tended to overestimate snow depth, especially during the fall season for the years 2014 and 2015 by approximately 20%.Similarly, the model overestimated snow depth for the year 2014 by approximately 20% and underestimated for the year 2015 by less than 5% at EML, with the mean RMS equal to 0.3108 ± 0.010 m.For additional modeled and observation monthly means for each site, refer to table S2.On average the model overestimated Reco at EML over the frozen season by 4.51 gC/m 2 /month and 1.12 at IMN gC/m 2 /month.The best-performing simulation for NEE resulted in an RMSE of 13.56 ± 1.61 and 16.10 ± 1.89 gC/m2/month at EML and IMN, respectively.For total annual carbon budgets for each study site, refer to figure S2.

Sensitivity analysis
The distributions of coefficient influence between the 14 parameters included in the sensitivity analysis varied at each site (figure 4).ALT was primarily influenced by the porosity of the fibric layer, a component of the organic soil that has undergone relatively little decomposition, at both sites (figure 4(a)).High n thawed resulted in a deeper active layer at EML. Maximum snow density and the n frozen were the most influential parameters on soil temperature at 10 cm at IMN (figure 4(b)), while soil temperature at 10 cm at EML was most influenced by the porosity of the fibric layer and the q10 coefficient for heterotrophic respiration (rhq10).GPP, Reco, and NEE were primarily influenced by the maximum rate of photosynthesis for ericoid shrubs and sedges, the two most dominant PFTs at both sites.The increase in the maximum rate of photosynthesis (cmax) increased ecosystem carbon sequestration through an increase in GPP that was partially offset by an increase in Reco.The increase in Reco was driven by (1) an increase in autotrophic respiration resulting from increased vegetation productivity and (2) an increase in heterotrophic respiration resulting from increased litterfall and associated availability of active soil carbon pool.Reco was also secondarily influenced by the q10 coefficient for heterotrophic respiration.

Soil moisture
In this study, we analyzed how the linkages between permafrost, hydrology, and ecosystem CO 2 dynamics are modulated by drainage conditions and soil properties at two tundra study sites in Alaska.While DVM-DOS-TEM showed agreement with field observations at both study sites, discrepancies between observations and model simulations remained for snow depth and soil moisture at 5 cm depth.Modeled soil moisture at 5 cm depth was at both sites, which may be influenced by the uncertainty in measurements of soil moisture, as sensor calibration is known to be challenging in porous organic horizons [73].Variations in soil moisture at EML from 2005 to 2009 could be partly attributed to inconsistencies in the climatic drivers, specifically precipitation [74,75].During these years, observed precipitation was substantially higher than CRU precipitation, which could explain the underestimation of soil moisture by DVM-DOS-TEM for these years (figure S1).Observations at EML showed that factors related to permafrost thaw, such as subsidence could lead to changes in drainage conditions, soil porosity, and soil moisture and could contribute to uncertainties within land models to accurately represent soil hydraulic properties [30,31,34,42,43,[76][77][78].While the underestimation of thawed season soil temperature at 10 and 40 cm depths at both sites could not be associated with bias in the air temperature, it may be due to a potential underestimation of soil moisture at these depths.

Porosity of fibric layer
Our analysis revealed the parameters which significantly influence the dynamics of soil cryohydrology and carbon shown in figure 4. The porosity of the organic layer, particularly the fibric layer, emerged as the critical parameter influencing soil moisture, soil temperature, and ALT across both study sites [79,80], the exception of soil temperature at the IMN site (figures 4(a)-(c)).Interestingly subsurface parameters such as bulk density and thermal conductivity showed a relatively minor influence on soil temperature at 10 cm and soil moisture at 5 cm (figures 4(b) and (c)) [81][82][83].The ability of a soil layer to store water is determined by its porosity.In peatland ecosystems, high porosity results in a more saturated organic layer at both poorly drained (EML) and welldrained (IMN) sites.Porosity influences soil temperature through its impact on soil thermal conductivity [84].The empirical studies have demonstrated that while soil moisture positively affects soil thermal conductivity, this effect is more constrained in organic soils than in mineral soils [85].The soil thermal regime was affected at both sites by n thawed and n frozen [86,87].At IMN, soil temperature at 10 cm is more sensitive to changes in n frozen than at EML, suggesting a strong influence on subsurface thermal conditions [88].Furthermore, the sensitivity analysis revealed that parameters associated with snow density and n frozen were the most influential on soil temperature at multiple depths at both sites.This emphasizes the insulating effect of snow on soil (table S3), as well as the effect of vegetation on surface temperature via nfactor [86,89].These results underscore the intricate interplay between seasonal n-factors, snow dynamic and soil structure, and their ultimate impacts on soil temperatures in the model.This suggests a continuing shift at the EML site from climate-driven ecosystemprotected to ecosystem-protected [18].The IMN site is located in the climate-driven permafrost zone and also exhibits the impact of fibric layer, where soil temperatures are mainly influenced by snow and n frozen , indicating an ongoing transition towards climatedriven and ecosystem-protected permafrost.

Ecosystem carbon balance
The model was able to represent the spatio-temporal patterns of carbon fluxes observed at both studied sites.The seasonal pattern of monthly GPP, Reco, and NEE showed general agreement between measured and model outputs across all years at both sites (figure 3).Both GPP and Reco were largely influenced by parameters driving vegetation productivity (i.e.cmax, rate limiting parameter of GPP, figure 4).The influence of vegetation productivity on Reco resulted from (1) the direct relationship between vegetation productivity and autotrophic respiration [90] and (2) the positive relationship between vegetation productivity and litterfall and between litterfall and soil organic matter decomposability [91].Annual means in NEE observations from both sites showed differences in magnitude and trends (figures S2(c) and (f)).Observations of yearly NEE fluxes indicate a net carbon source with an annual mean of 50 ± 93 g CO 2 C m −2 y −1 at EML [31,92].In contrast, our model results show EML to be a net carbon sink, with a mean of −85 ± 35 g CO 2 C m −2 y −1 .The overestimation of NEE uptake in the model can be attributed to the overestimation of GPP during fall shoulder season months (figure 3(b)).Based on both observations and model simulations, IMN could not be distinguished from carbon neutral across measurement years (−21 ± 33 g CO 2 C m −2 y −1 ), switching from source to sink across the study period years (2009-2016), which was consistent with other studies at IMN [48].

Effect of climate on GPP
While the direct effect of rising air temperatures has been linked to an increase in vegetation productivity in warming experiments across northern ecosystems [93,94], the increased risks of extreme warming events, droughts, disturbances, and heightened competition among plant communities that have been associated with warming can also result in reduced vegetation productivity [95,96].The correlation analysis showed that GPP and air temperature have a significant negative correlation at IMN (table S4).While air temperature is positively linked to vegetation productivity in DVM-DOS-TEM [57], this effect is offset at IMN by the fact that during the period of analysis, warm growing summers are associated with low precipitation (figure S3(a)).Increased water demand resulting from the effect of warm temperatures on photosynthesis, and low water availability resulting from low precipitation, led to high water stress (estimated by the ratio of actual and potential evapotranspiration, EET:PET, figure S3(b)) and a decrease in GPP (figure S3(c)) during warm summers.This effect is not as significant at EML because the level of precipitation is generally higher than at IMN (figure S3(a)) and the level of water stress is lower (i.e. higher EET:PET ratio, figure S3(b)).Our results suggest that the response of vegetation productivity to warming will strongly be modulated by the level and trends in water availability, which is controlled by precipitation and permafrost hydrology in high latitudes.

Effect of climate and drainage on Reco
In the permafrost region, soil respiration typically increases as a result of rising air temperature through two main processes: (1) the increase in soil temperature resulting in faster decomposition rate [38] and (2) the increase in the amount of unfrozen organic matter available for decomposition due to permafrost thaw [97].Both processes are represented in DVM-DOS-TEM [58,61] and are the main drivers of the differences in Reco between sites (figure S4(a)), with higher overall Reco at EML located in a warmer boreal climate, than at IMN located in a colder arctic climate.
At a site level, soil temperature was not the main driver of interannual variation in Reco due to its limited range of variation within the 30 year period of analysis within each site (figure S4(a)).this scale, water supply was found to be greatly modulating the influence of soil temperature on ecosystem respiration, as shown by the following two pieces of evidence.Firstly, we found that at IMN, the inverse relationship between ecosystem respiration and air temperature (table S4).We investigated how this relationship translates to the soil temperatures at 10 cm depth.We found that an outlier corresponding to year 2007 emerged as the driest year (figures S4(a) and (b)).The condition that led to unusually high soil temperatures-via the influence of soil moisture on thermal conductivity [49]-and low respirationvia the effect of soil moisture on heterotrophic and autotrophic respiration (figure S4(c)) [57,98].Upon exclusion of the year 2007 outlier, the correlation between Reco and soil temperatures transitioned to a positive relationship.This phenomenon underscores the pivotal influence of precipitation on soil thermal dynamics and the overall ecosystem carbon balance.Secondly, we restricted the drainage conditions at EML by reducing baseline flow and runoff, which resulted in a shallower water table in comparison to the well-drained IMN site (figure S4(d)).Within each site, our model simulations showed that shallow water table limited Reco (figure S4(d)), due to the associated reduction of aerobic conditions in the soil column [38].However, recent synthesis of anaerobic incubation experiments showed that anaerobic production of CO 2 and CH 4 can offset the relative decrease in aerobic CO 2 production when soil saturation increases [38].It is thus critical for large-scale ecosystem models like DVM-DOS-TEM to include anaerobic processes when evaluating the impact of permafrost thaw in lowlands [99].
Lastly, winter respiration was strongly correlated to air temperature at both sites (table S4), following findings from the recent flux meta-analysis [100], and emphasizing the importance of ecosystem models parameterized for high-latitude ecosystems to account for winter respiration in regional carbon assessments.

Implications
The changes in soil thermal and hydrological regimes play a pivotal role in shaping both belowground and aboveground carbon dynamics.Soil properties and thermal regime affect soil respiration, and vegetation productivity, impacting NEE.The variance in NEE at the two sites further underscores the modulation of ecosystem carbon dynamics by climatic conditions, permafrost, and hydrological states.Our findings suggest that while climate warming typically enhances vegetation productivity in northern ecosystems, this trend is not always positive, highlighting the importance of precipitation and the complex interactions between permafrost, hydrology, and ecosystem carbon dynamics.Additional research and data are required to determine whether this is a transient occurrence or a longer-term pattern.
Finally, it is critical to consider the broader implications of our findings for predicting ecosystem responses to climate change.Climate projections carry over uncertainty in the projections of ecological processes, not only over long-term periods (e.g.[101]), but also across areas with scarce weather monitoring networks [102].Our analysis has shown that an alteration of the linkage between air temperature and precipitation in climate forcing can have a substantial influence on carbon dynamics and the response to warming.As air temperatures continue to rise, we anticipate shifts in the relationship between ALT and ecosystem carbon dynamics, contingent on permafrost thaw stages and hydrological changes.This underscores the delicate balance between environmental drivers and ecosystem responses, which can be further complicated by local and regional climate variability and the capacity of models to capture complex ecological processes.

Conclusions
This analysis revealed the importance and feasibility of using a parameter sensitivity approach to diagnose model simulations of Arctic tundra ecosystems.Both site observations and model simulations indicated the well-drained Arctic tundra site (i.e.IMN) to be carbon neutral over the period from 2008 to 2016.In contrast, the poorly drained boreal tundra site (i.e.EML) is identified as a net emitter of carbon according to observations, while model projections depicted it as a carbon sink.However, DVM-DOS-TEM does not yet represent carbon loss related to anaerobic processes, thus underestimating carbon loss in poorly drained sites.
Vegetation productivity drives changes in GPP and Reco and is influenced by both temperature and moisture levels.We highlight the complex interplay between warming air temperatures, moisture availability, and permafrost conditions, revealing that while warming generally increases vegetation productivity, this response is highly sensitive to precipitation and, therefore, water stress.
Soil aerobic respiration typically decreases with increasing soil water saturation and increases with warmer soil temperatures due to faster decomposition and increased availability of organic matter from thawing permafrost, which underscores the significant influence of climate and drainage on Reco.The inverse relationship between Reco and air temperature at the IMN site was traced back to an anomalously dry year, indicating the substantial role of precipitation in soil thermal regulation and overall ecosystem carbon balance.
While warming is expected to raise soil temperatures, the effects on carbon dynamics are modulated by factors such as soil saturation and permafrost thaw stages.These insights highlight the necessity for models (1) represent and benchmark interactions between soil thermal and hydrological regimes and their influence on the ecosystem carbon cycle, and (2) integrate anaerobic processes and winter respiration to predict ecosystem responses to climate change.Further improvements to better capture changes to the subsurface soil due to thaw, such as subsidence, are still needed, highlighting the complexity of capturing complex permafrost, hydrology, and ecosystem carbon dynamics.However, our study also suggests that by forcing drainage conditions to reflect the hydrological consequences of subsidence could be used to investigate its ecological impacts.It should be noted that the influence of permafrost on methane release was not considered in this analysis, as the DVM-DOS-TEM model does not yet incorporate methanogenesis.Nevertheless, it is expected that methane release due to the high moisture levels at EML would counterbalance the modeled carbon absorption [9].
These findings build on the understanding of the complex network of interactions between permafrost and ecosystem carbon dynamics across highlatitude sites and emphasize the importance of model validation against site-level observations, specifically through an in-depth multi-layer analysis.Finally, this analysis emphasized the need to further finetune driving parameters related to soil structure and soil characteristics and to better represent the linkage between permafrost thaw and ecosystem carbon dynamic in poorly drained conditions to reduce model uncertainty and allow for more accurate projections of landscape change in response to ongoing and future climate change in the permafrost region.

Figure 1 .
Figure 1.Study sites, (a) Eight Mile Lake and Imnavait Creek, used within this study with permafrost extent across Alaska.(b).The Imnavait site is located within a region of continuous permafrost and is characterized as a well-drained site, while (c) Eight Mile Lake is located within discontinuous permafrost and is characterized as a moderately-drained site.Reported Net Ecosystem Exchange totals (b), (c) are from site observation flux data where EML is shown to be a net source, while IMN is relatively neutral, neither sink nor source.Permafrost distribution source:[47], Imnavait Site Photo source: AmeriFlux.Adapted with permission from Greg Friske.(Main map layer) Adapted with permission from[47].Adapted with permission from Sebastian Biraud.

Figure 2 .
Figure 2. Model-data comparison of physical variables at Imnavait (IMN) and Eight Mile Lake (EML) sites for active layer depth (a), (b), soil temperature at 5 cm depth (c), (d), soil temperature at 10 cm depth in (e), (f), soil temperature at 40 cm depth in (g), (h), soil moisture at 5 cm depth in (i), (j), and snow depth (m) in (k), (l).Lighter gray shaded areas show the spread of model output simulations, with model means noted in a solid blackline, and ±1 standard deviation shown in light gray.Site observations are shown in red for each plot.

Figure 3 .
Figure 3. Model-data comparison of the ecosystem carbon balance at IMN and Eight Mile Lake sites for Gross Primary Productivity (GPP) (a), (b), Ecosystem respiration (Reco) (c), (d), and Net Ecosystem Exchange (NEE) (e), (f).Lighter gray shaded areas show the spread of model output simulations, with model means noted in a solid blackline, and ±1 standard deviation shown in light gray.Site observations are shown in red for each plot.Ecosystem respiration was computed as the sum of heterotrophic and autotrophic respiration, while NEE is the difference between Reco and GPP.

Figure 4 .
Figure 4. Model parameter importance rank for all model output variables: (a) active layer thickness, (b) Soil temperature at 10 cm, (c) Soil moisture at 5 cm, (d) Gross Primary Productivity (GPP), (e) Ecosystem Respiration (Reco) for each study site.Parameter importance was calculated from model simulations (n = 1000).

3.2. Data model comparison: CO 2 fluxes
Modeled seasonal maximum GPP matched observations at IMN well for all years (2008-2015), with an RMSE of 21.39 ± 2.14 gC/m 2 /month (figure 3(a)), however, the model overestimated GPP during the thawed season of 2008 by an average of 50.7 gC/m 2 /month.EML had a slightly greater total GPP, both modeled and observed, than IMN, with a mean RMSE of 26.91 ± 4.2.Modeled monthly GPP at IMN averaged 94.77 ± 6.26 gC/m 2 /month, while the observed mean GPP at IMN was 67.99 gC/m 2 /month.At EML, monthly modeled GPP averaged 115.46 ± 8.13 gC/m 2 /month, while observations averaged 108.86 ± 23.02 gC/m 2 /month.There was a general agreement between modeled and observed Reco at both IMN and EML (figures 3(c) and (d)) with a mean RMSE of 20.0 ± 5.0 gC/m 2 /month at IMN and 13.4 ± 1.84 gC/m 2 /month at EML. Model simulations of NEE showed a minor overestimation during the frozen season winter months (higher CO 2 emissions) at EML with a mean NEE of 26.98 + 4.83 gC/m 2 /month compared to the observational mean of 22.01 + 5.67 gC/m 2 /month, and during the fall months at IMN, with a modeled mean of 19.21 + 7.49 gC/m 2 /month compared to observation mean of 11.90 + 4.45 gC/m 2 /month, due to overestimations of computed Reco over the corresponding months.