Permafrost thermal conditions are sensitive to shifts in snow timing

Changes in snow precipitation at high latitudes can significantly affect permafrost thermal conditions and thaw depth, potentially exposing more carbon-laden soil to microbial decomposition. A fully coupled process-based surface/subsurface thermal hydrology model with surface energy balance is used to analyze the impact of intra-annual variability in snow on permafrost thermal regime and the active layer thickness. In the four numerical scenarios considered, simulations were forced by the same meteorological data, except the snow precipitation, which was systematically altered to change timing of snowfall. The scenarios represent subtle shifts in snow timing, but the snow onset/melt days, the end of winter snowpack depth, and total annual snow precipitation are unchanged among scenarios. The simulations show that small shifts in the timing of snow accumulation can have significant effects on subsurface thermal conditions leading to active layer deepening and even talik formation when snowfall arrives earlier in the winter. The shifts in snow timing have a stronger impact on wetter regions, especially soil underneath small ponds, as compared to drained regions. This study highlights the importance of understanding potential changes in winter precipitation patterns for reliable projections of active-layer thickness in a changing Arctic climate.


Introduction
A huge amount of frozen organic carbon is stored in Arctic permafrost regions (Schuur et al 2015, Zimov et al 2006, most of which lies in the upper 3-4 m of soil (Tarnocai et al 2009, Hugelius et al 2014. Widespread permafrost degradation in the high-latitude regions has already been recognized in multiple field studies, e.g. Lachenbruch and Marshall (1986), Romanovsky et al (2002), Osterkamp (2003), Hinzman et al (2005), Jorgenson et al 2006), Osterkamp (2007b), Wu and Zhang (2008), Batir et al (2017). We refer to permafrost degradation as thawing of permafrost or increase in soil temperature within permafrost (Dobinski 2011). Degradation of near-surface permafrost can result in decomposition of soil organic matter and then possibly release to the atmosphere, which could cause further increase in the temperature (Schaefer et al 2014, Koven et al 2011, Schuur et al 2008. A number of factors influencing permafrost degradation have been identified such as air temperature, seasonal snow cover, soil stratigraphy, vegetation, soil moisture, and terrain (Zhang and Stamnes 1998, Osterkamp 2007b, Atchley et al 2016, Nicolsky et al 2017. However, due to the complex Arctic environment, it is often unclear which factor or combination of factors contributed most to permafrost degradation.
Snow is one of the most important physical variables of the high-latitude regions. Snow has low thermal conductivity and high albedo and thus acts as an effective insulator between the atmosphere and the ground (Callaghan et al 2011b). Theoretical understanding of the role of snow cover in soil warming is well established (Goodrich 1982, Callaghan et al 2011b, Westermann et al 2013, Langer et al 2013, Ekici et al 2015, Domine 2011, and several field and modeling studies have highlighted snow as an important climate factor in permafrost regions (Zhang 2005, Zhang et al 2018, Park et al 2015, Park et al 2014, Morse et al 2012, Ling and Zhang 2007, Zhang et al 1996, Lawrence and Slater 2010, Gouttevin et al 2012, Jafarov et al 2014, Domine et al 2016, Gisnås et al 2016, Calonne et al 2011. Snow cover dynamics such as variations in its intensity, frequency, duration, characteristics, and the timing of onset, melt, and maximum snowpack depth can significantly affect permafrost thermal regime (Goodrich 1982, Zhang and Osterkamp 1993, Osterkamp and Romanovsky 1996, Stieglitz et al 2003.
Climate models have projected considerable changes in precipitation during the 21st century in the high latitudes of the Arctic. These changes, likely increases but highly uncertain, in precipitation are expected to be more profound in winter than in summer (Kattsov et al 2005, Kattsov et al 2007, and may include increases in fall and early winter snow precipitation (Kattsov et al 2007, Collins et al 2013, Bintanja and Selten 2014, AMAP 2017. Given the strong control of snow on permafrost thermal conditions, changes in snow precipitation have potential to affect permafrost degradation and active layer thickness. In recent decades, a number of field studies have indicated that snow cover has played an important role in permafrost soil warming. For example, increased permafrost temperatures and thermokarst (soil subsidence due to the melting of ice) formation have been observed in field studies (Osterkamp 2007a, Osterkamp et al 2009, Sannel et al 2016, and in snow manipulation experiments (Johansson et al 2013, Lafrenière et al 2013, O'Neill and Burn 2017a. Those studies show a significant positive correlation between winter precipitation and soil warming. Complementing those observational studies, snow manipulation experiments in the Arctic regions have a relatively long history. These experiments have been used as a tool to understand potential effects of snow on soil warming, subsidence, and carbon dynamics under a changing environment; see, for example, Natali et al (2011) and references therein.
Although observational and snow-manipulation studies have clearly established the importance of snow cover on permafrost thermal conditions, understanding of the relative importance of various snow-related processes is also needed. In particular, the effects of potential shifts in intraseasonal patterns of snow precipitation on permafrost is difficult to untangle from other influences in observational and snow manipulation studies. Models capable of simulating the current state of permafrost in conjunction with well-designed numerical experiments allow us to study the effects of important processes influencing permafrost degradation in different climate scenarios. Here we take advantage of new modeling capability that couples key thermal hydrological processes across the surface and subsurface to investigate how intra-annual variability in snow precipitation (enhanced snow accumulation during fall and early winter) would affect permafrost degradation and active-layer thickness. We use the recently developed Advanced Terrestrial Simulator (ATS), which provides integrated surface and subsurface thermal hydrology models including surface energy balance, snow processes, overland flow, and advanced representations of interacting flow and heat transport in variably saturated, variably frozen soil (Painter et al 2016). We isolate the influence of snow accumulation timing by repeating a 365-day time series of other surface forcing variables such as air temperature, rain precipitation, relative humidity, wind speed, and incoming shortwave and longwave radiation. Four simulation scenarios have been examined on a two-dimensional transect with a small shallow pond (1 m deep) at the valley bottom. We force the model with observed climate data from a continuous permafrost site near Utqiaġvik, Alaska. The model domain is considered to be representative of cold continuous permafrost tundra hillslope environments. It is important to note that the total annual average snow precipitation, snow onset/melt day, and the end of winter snowpack depth are unchanged in all scenarios.
Models have previously been used in simulating the effects of snow cover on permafrost soils and active-layer thickness (e.g. Stieglitz et al 2003, Ling and Zhang 2007, Ge and Gong 2010, Atchley et al 2016, O'Neill and Burn 2017b. Using simplified one-dimensional models, those studies mostly focused on snowpack depth, as opposed to this study where we focus on the effects of enhanced snow accumulation in early winter with other snow-related factors held constant. Here we keep the winter season length, number of snowy days, maximum snow depth, and total annual snow precipitation the same and simulate solely the impacts of temporal shifts in snow accumulation, herein referred to as 'snow timing' . Of the previous works just mentioned, only O'Neill and Burn (2017b) examined the effect of changing the duration of maximum snow accumulation. They simulated dry (no surface inundation) tundra with different snow accumulation periods using a one-dimensional subsurface heat transport model. Here we simulate in two dimensions at the hillslope scale and focus on the interaction between surface inundation and snow timing, which requires a more complete cryohydrology model that includes surface and subsurface flow, lateral heat transport, surface energy balance, and important snow processes such as dynamic snow density, aging, and compaction.

Domain description
The simulation domain represents a two-dimensional cross-section of tundra hillslopes in cold continuous permafrost regions (figure 1). The domain is 100 m long. The surface topography is a hillslope with 10% slope and a parabolic-type valley bottom on the left (x = 0−10 m). For discharge, a spill point is set 1 m above the datum (Z 0 in figure 1). The bottom boundary (50 m below the datum Z 0 ) of the domain is flat. The model domain is resolved by 100 × 80 grid cells with horizontal size of 1 m and vertical discretization that increases with depth, starting at 2-5 cm in the top 3 m. The domain consists of three vertical layers: moss (2 cm), peat (8 cm), and mineral soil (rest of the domain). The surface and subsurface thermal and hydraulic boundary conditions are discussed later. It should be noted that the thermal signal from the surface does not penetrate deeper than 15-20 m, which is well short of the bottom boundary, thus avoiding boundary artifacts.
For observations only, the top 25 m model domain is split into six vertical zones. The label, thickness, and depth (from the surface) to base and top of each zones are provided in table 1. We further subdivide into two horizontal zones: soil beneath the inundated region at the valley bottom (wetland, length 0-5 m), and soil beneath the drained region (upland, length 5-100 m).

Cryohydrology simulator
In recent years, cryogeohydrology models have become available ( . However, those models only represent subsurface processes and must be driven by surface boundary conditions obtained from empirical correlations to represent the effects of meteorological forcings and surface thermal hydrology processes. More recently, integrated surface/subsurface permafrost thermal hydrology models started to appear, for example, the Advanced Terrestrial Simulator (ATS) (Coon et al 2019). This emerging class of models couple the subsurface processes represented in cryogeohydrology models to overland flow, surface energy balance, and snow processes.
Our numerical experiments used the Advanced Terrestrial Simulator (ATS, version 0.88), which provides a comprehensive process-based permafrost thermal hydrology modeling capability (Coon et al 2019). The permafrost configuration of ATS solves fully coupled conservation equations of mass (a modified form of Richards equation) and energy for variably saturated and variably frozen flow in the subsurface. That subsurface representation is coupled to 1) a surface flow system represented by a modified diffusion wave equation for flow and energy, 2) surface energy balance model driven by meteorological forcing data, and 3) a phenomenological snow distribution model (not exercised in this study). ATS represents important physical processes such as snow processes (thermal conduction, aging, compaction, effects of depth hoar, and spatial distribution), hydrological processes (lateral surface and subsurface flows), soil/surface thermal processes (advective heat transport), and cryosuction. More details about implementation of mathematical models in ATS and its capabilities can be found in Painter et al (2016) In previous studies, ATS has been evaluated against borehole temperature data for Utqiaġvik, Alaska in one-dimensional (Atchley et al 2015) and two-dimensional (Abolt et al 2018) simulations. A very recent study evaluated ATS against multiple field observations (such as snow depth, water table, soil temperatures, and evaporation) in multiyear threedimensional simulations for polygonal tundra in Utqiaġvik, Alaska (Jan et al 2019). These studies have shown ATS to provide observations-consistent integrated thermal hydrology modeling capability in continuous permafrost regions.

Model initialization
The model is initialized using a multi-step procedure described as follows: 1) Run a one-dimensional column simulation to establish water table, 2) Freeze the water table from below by setting a constant boundary temperature of −6 • C (1000-year simulation), 3) Map the one-dimensional ice-table to twodimensional model domain, 4) Run the model using observed forcing data for 1985 to obtain approximate annually periodic steady state (300-year loop), and 5) Drive the model with observed meteorological data from 1986 to 2015. This completes the model initialization process. In all simulations, the last day of the previous step serves as the initial condition for the next step. Steps 4 and 5 perform fully integrated surface/subsurface thermal hydrology with surface energy balance simulation using Daymet data (Thornton et al 2018) for Utqiaġvik, Alaska, adjusted for snow undercatch (30% increase) as described in Atchley et al (2015). This spinup process provides initial conditions for four numerical experiments considered in this work.

Climate scenarios
For the basecase climate scenario, denoted by P 0 , we constructed a one-year time series by averaging each day of the observed meteorological data in the years 2010 to 2015 from the Barrow Environmental Observatory (BEO) located near Utqiaġvik, Alaska. This 2010-2015 meteorological forcing data were compiled from several sources. More details are provided in the supplemental material (section S1: Climate scenarios).
In order to generate three additional scenarios with earlier snow accumulation compared with the basecase, we first split the basecase winter season into two time intervals, before and after the time of maximum snowpack depth (≈37 cm). We refer to these as the snow-accumulation period and the latewinter period. The new scenarios, denoted P 30 , P 60 , and P 90 , are then generated by shortening the snowaccumulation period and inserting before the latewinter period a plateau period during which snow depth is approximately constant. These scenarios are shown schematically in figure 2. Daily snow precipitation for the scenario denoted by α = 0, 30, 60, 90 can be expressed mathematically as Here, P [m,n] means basecase precipitation in the index interval between m and n and zero otherwise, where the index corresponds to a winter day (air temperature <0 • C). The notations l and N denote the lengths of the basecase accumulation period and the entire winter season, respectively. It is worth pointing out that the late-winter period (last term in the equation) remains unchanged.
In this weighted average technique, the basecase precipitation is reproduced by setting α = 0, β 1 = 1, and β 2 = 0. In snow precipitation scenarios P 30 , P 60 , and P 90 , the parameters β 1 and β 2 were chosen so that the accumulation phase is shortened by roughly 30, 60, and 90 days relative to basecase snow accumulation period. It is important to note that we have systematically chosen β 1 and β 2 such the maximum snowpack depth (37 ± 2 cm) and the end of winter snowpack depth (35 ± 0.25 cm) do not deviate significantly from those of P 0 (the basecase). Additionally, the annual average snow precipitation (220 mm per year) is unchanged among scenarios as shown in figure 3(bottom left). The time series of the snow precipitation in our four scenarios is presented in figure 3(top). The scalars β 1 and β 2 in scenarios P 30 , P 60 , and P 90 are (1.03, 0.735), (1.08, 0.635), and (1.22, 0.512), respectively. Specifically, the maximum snowpack depths (37 ± 2 cm) were obtained on 214th, 188th, 161th, and 131th day in scenarios P 0 , P 30 , P 60 , and P 90 , respectively, indicated by vertical lines in figure 3(bottom right). It is important to note that the snow onset/melt days, end of winter  . The cumulative snow precipitation and winter snowpack depth are shown in the bottom left and right, respectively. Snow onset date, end of winter snowpack depth, and the total annual snow precipitation are unchanged. The horizontal line in the bottom right corresponds to 37 cm snowpack depth, approximately the maximum snowpack depth in the basecase. The gray horizontal band represents ±2 cm variation. The vertical dashed lines show the day when snowpack depth is 37 ± 0.5 cm. The notation S0, S30, S60, and S90 represents simulations forced with P0, P30, P60, and P90, respectively. snowpack depth, and the annual time series of other climate factor (air temperature, rain precipitation, wind speed, relative humidity, and incoming longwave and shortwave radiation) are the same in all scenarios.

Simulations description and model inputs
Four numerical experiments denoted by S 0 , S 30 , S 60 , and S 90 were performed corresponding to four climate scenarios P 0 , P 30 , P 60 , and P 90 , respectively. In all simulated scenarios, models were initialized with the same permafrost conditions that were obtained by the end of spinup, which is the last day of the year 2015. To remove the initial influence of climate data, we repeatedly forced the model with the four climate scenarios to obtain an equilibrium state, which was achieved after approximately 10year loop. Our focus in this work is to study the influence of intra-annual variations in snow precipitation on permafrost thermal regime, therefore, simulations ignored the spatial variability of snow. Inputs to our model include meteorological forcing data, thermal and hydraulic boundary conditions for both surface and subsurface, soil stratigraphy, threedimensional mesh, intrinsic permeabilities, porosities, wet/dry soil thermal conductivities, and water retention curves parameterized by van Genuchten model. Soil thermal and hydraulic properties used in our study are provided in table S1 in the supplemental document (stacks.iop.org/ERL/15/084026/mmedia), and are consistent with those of Jan et al (2019).

Boundary conditions
A constant temperature of −6 • C is applied at the bottom boundary (Romanovsky et al 2010). The simulations use no-flow boundary conditions on the vertical sides of the subsurface domain (see figure S1 in the supplemental document). The top-surface is subject to transient (daily) observed meteorological data. The meteorological forcing data used in our simulations include air temperature, snow precipitation, rain precipitation, relative humidity, wind speed, and incoming longwave and shortwave radiations. The surface hydraulic boundary condition is a spill point at 1 m above the datum (Z 0 in figure 1), which allows outflow when the simulated ponded depth is above that value and no flow otherwise. It should be noted that the pond water level can fluctuate due to evaporation or infiltration or exfiltration, but always stays below or at the level of the spill point.

Results
To evaluate the effects of intra-annual variations in snow precipitation on permafrost temperatures, active-layer thickness, and discharge, results from S 0 (the basecase) are compared to those obtained in scenarios S 30 , S 60 , and S 90 . Figure 4 shows soil temperature profiles in the top 4 m of the model domain for the four scenarios. The snapshots were taken on the last day of winter (day before the air temperature goes above freezing) from the last simulated year. The rows (top to bottom) correspond to scenarios S 0 , S 30 , S 60 , and S 90 , respectively. Results indicate that intraannual variations in snow caused significant warming in soils in wetter regions (area underneath small pond). In the S 90 scenario, a small region of the soil beneath the pond remains unfrozen at the end of the winter, potentially indicating the formation of a talik.
To explicitly demonstrate the role of intra-annual variations in snow on permafrost warming, we divided the top 25 m of model domain into six vertical and two horizontal zones as described in section 2.1. The annual average ground temperature for the active layer, shallower permafrost, and shallow permafrost zones are displayed in figure 5. The columns correspond to two horizontal zones: subsurface beneath inundated surface (left) and drained upland region (right); and rows correspond to the top three vertical zones. Results for the remaining permafrost zones (deep, deeper, and deepest) are provided in the supplemental document (figure S2). These results indicate that ground temperatures increase as early winter snow precipitation increases. For instance, in scenario S 90 as compared to S 0 , the annual average ground temperature increased by 1.43, 2.03, and 1.35 • C in the active layer, shallower, and shallow permafrost zones underneath the pond center (x = 0 m), respectively, (figure 5(left)). It can be seen that shallower permafrost (permafrost right beneath the active layer) was more affected by the snow timing as compared to other permafrost zones. An increasing trend in the soil temperatures in the upland region is illustrated in figure 5(right). Soil temperature in the permafrost deeper than 8 m, in the uplands, is not significantly affected by the intra-annual variability of snow in scenarios S 30 and S 60 . In scenario S 90 , the thermal signal of the snow timing penetrates to depth 15 m in the upland region. However, in the wetland regions, the top 15 m permafrost is affected by the snow timing in all scenarios. Moreover, irrespective of the numerical scenario, simulations show that the soil beneath the pond is 3-4 • C warmer than the upland region, highlighting the effects of soil moisture and small water bodies (<1 m) on the permafrost soil thermal regime. The colors black, green, blue, and red correspond to scenarios S 0 , S 30 , S 60 , and S 90 , respectively, and are consistent among figures. Figure 6 shows time series of the ground temperature at depth 56 cm (location of permafrost table in the reference case) at x = 0 m (center of the pond) and at x = 50 m (center of the model surface domain). In scenario S 90 , ground temperatures beneath the center of the pond are warmer than the reference case by 5 • C in the winter and 0.5 • C warmer in the summer. In the upland, winter temperatures show an increase by 4 • C, but summer temperatures are not affected by intra-annual variations in snow precipitation. This suggests that projected increase in early winter snow accumulation will have a more pronounced effect on permafrost temperatures beneath inundated regions (wetter soil) than drained regions (drier soil).
Time series of depth-to-base and depth-to-top of the thaw layer at locations x = 0 m and x = 50 m are shown for all scenarios in figure 7. Soil underneath the pond mostly freezes from the bottom and thaws from the top, unlike the drained areas that freeze up from both the top and the bottom. The first day of the time axis on figure 7 corresponds to the fall onset date, the first day when the air temperature falls below the freezing point, 0 • C. There is a significant time lag between fall onset and freezing beneath the pond. In the scenario S 0 , the soil column beneath the pond completely freezes, but not until mid-March. In scenarios S 30 and S 60 the ground stays unfrozen until April and May, respectively. This significant delay in complete freeze-up and increased soil temperatures in several permafrost layers (as shown in figure 5) are indications of permafrost warming.
Significantly, enhanced snow accumulation in early winter resulted in the formation of a closed talik in scenario S 90 . That is, soil beneath the pond from 5 cm to 55 cm depth remains unfrozen throughout the year. The development of a talik is also evident Finally, results from our simulations show that more snow accumulation in early winter leads to an increase in cumulative discharge, the result of a decrease in cumulative evaporation (figures S3 and S4 in the electronic supplement). However, the effect is small, less than 5% of the discharge in the S 90 . The water content in the active layer also decreased slightly in the non-basecase scenarios ( figure S4(right)). As expected, the water content in the remaining vertical zones (e.g, shallower permafrost etc) was not affected by the snow timing. It is important to note that the total water mass in the system is conserved in all cases.

Summary and discussion
The role of snow in permafrost warming has been reported in many field observations and in modeling studies (Stieglitz et al 2003, Osterkamp 2007a, Lafrenière et al 2013. However, the thermal response of permafrost is governed by tightly coupled thermal and hydrological processes on the surface and in the subsurface, which makes it difficult to untangle sensitivities and creates uncertainties in projections of permafrost degradation. The emerging class of  integrated surface/subsurface thermal hydrology models like the ATS provide capability for developing process-based understanding of those sensitivities and dependencies. Here we performed numerical snow manipulation experiments using ATS to explore sensitivity of permafrost thermal regime and active layer dynamics to changes in the intra-annual variability in snow precipitation patterns. We find that the depth of the active layer, annual average temperature, and winter temperatures are all sensitive to intra-annual patterns in snow precipitation, with more snowfall in early winter resulting in warmer subsurface conditions for the same air temperature and end-of-winter snowpack. This is consistent with field studies that have shown that at some locations permafrost has warmed without any warming trends in air temperature, for example, Sannel et al (2016). Winter temperatures are more sensitive than summer conditions. For example, in scenario S 30 , simulated ALT is similar to the basecase ALT but permafrost winter temperature are 2 degrees warmer. That result is consistent with the findings of Lafrenière et al (2013) that changes in the snow depth has stronger effects on winter soil temperatures than on the ALT.
Soil beneath small shallow ponds is warmer than the drained (upland) regions in our numerical experiments. For the conditions of our simulations, annual average soil temperature beneath the shallow pond is about 3.5 • C warmer than the drained regions. Winter temperatures beneath ponds are more sensitive than annual average temperatures to changes in snow timing. In the S 90 scenario, a talik formed beneath the pond despite no imposed change in air temperature or maximum snowpack depth. Talik formation has the potential to further accelerate permafrost thaw because taliks can retain and transmit more energy (Connon et al 2018) and expand both vertically and laterally (Devoie et al 2019).
Greater sensitivity under small ponds can be understood by considering the competition between heat conduction and latent heat effects. Water bodies conduct more heat than dry soil due to the higher thermal conductivity of water, leading to colder soil temperatures in winter. However the energy required to melt ice (the latent heat of fusion) works in the opposite direction of that effect. The strong insulating effects of snow and enhanced snow accumulation in fall can considerably reduce winter heat conduction and allow the latent heat of fusion effect to become dominant, thus keeping soil under ponds warmer. Permafrost degradation underneath small ponds can cause early thermokarst formation as compared to regions without preexisting ponds (Langer et al 2016).
Our results show that snow timing has a modest effect on evaporation, discharge, and total water content in the system. Increased subsurface lateral flow caused by the prolonged unfrozen period and active layer deepening resulted in decreased active layer water content. Although the decrease in the active layer water content over a single year is not significant as compared to the bulk water, it may have a larger impact on the total water budget in the long-term.
An earlier study by O'Neill and Burn (2017b) investigated the effects of prolonging the duration of snow accumulation using one-dimensional heat transport simulations. Their results show considerable influence of snow timing on thaw depth (decrease by 1-3 m) for dry lands (no water inundation). While we find a similar relationship between annual average subsurface temperature and the length of the snow accumulation period (figure 5), in the dry region that effect comes mostly from changes in winter temperatures with little effect on thaw depths in the summer (figure 7). Active layer thickness in our simulations is more sensitive in inundated areas but still considerably less sensitive than the simulations of O'Neill and Burn (2017b), who only considered dry regions. Those differences in results are caused by different assumptions about total snow precipitation and different scenario definitions, and likely by model differences.
Climate models have projected changes in winter precipitation during the 21st century in the highlatitude regions (Kattsov et al 2007, Collins et al 2013, Bintanja and Selten 2014, AMAP 2017, with early winter season projected to experience enhanced snow. However, these projected changes in snow precipitation are uncertain and likely regionally nonuniform (Callaghan et al 2011a). Our finding of high sensitivity of the permafrost thermal regime to temporal shifts in snowfall implies those previously identified uncertainties in projected atmospheric forcing will propagate to greater uncertainty in permafrost degradation projections.
Our numerical experiments suggest that reliable projections of air temperature and total snowfall may not be sufficient to understand permafrost degradation and associated thermal and hydrological changes. Intra-annual variability in snow precipitation is also important and subtle shifts in those patterns may bring substantial hydrological change to Arctic ecosystems. Thus, knowledge of the changes in snowfall timing are also needed for reliable projections of permafrost degradation. These results also underscore the importance of snow process representations and reliable field observations of snow parameters (such as density, thermal conductivity, and compaction) throughout the winter period.

Acknowledgments
This work was supported by the Next-Generation Ecosystem Experiment-Arctic (NGEE Arctic) project.

Code/Data availability
The Advanced Terrestrial Simulator (ATS) Coon et al (2019) is open source under the BSD 3-clause license, and is publicly available at https://github.com/amanzi/ats. Simulations were conducted using version 0.88. The data that support the findings of this study are available from the corresponding author upon reasonable request. In addition, the forcing data, model input files, and jupyter notebooks for pre-and post-processing will also be made publicly available through the NGEE Arctic long-term data archive.