A new snow module improves predictions of isotope-enabled MAIDENiso forest growth model

The representation of snow processes in forest growth models is necessary to accurately predict the hydrological cycle in boreal ecosystems and the isotopic signature of soil water extracted by trees, photosynthates and tree-ring cellulose. Yet, most process-based models do not include a snow module, consequently their simulations may be biased in cold environments. Here, we modified the MAIDENiso model to incorporate a new snow module that simulates snow accumulation, melting and sublimation, as well as thermal exchanges driving freezing and thawing of the snow and the soil. We tested these implemen5 tations in two sites in East and West Canada for black spruce (Picea mariana ) ::::: (Mill.) :::::: B.S.P.) and white spruce (Picea glauca ) :::::::: (Moench) ::::: Voss) : forests, respectively. The new snow module improves the skills of the model to predict components of the hydrological cycle. The ::::::::::: MAIDENiso model is now able to reproduce the spring discharge peak and to simulate stable oxygen isotopes in tree-ring cellulose more realistically than in the original, snow-free version of the model. The new implementation also results in simulations with a higher contribution from the source water on the oxygen isotopic composition of the simulated 10 cellulose, leading to more accurate estimates : of :::::::: cellulose ::::::: isotopic :::::::::: composition. Future work may include the development of inverse modelling with the ::: this : new version of MAIDENiso to produce robust reconstructions of the hydrological cycle and isotope processes in cold environments.


Introduction
In boreal regions of Canada and Alaska, snow represents about 30-50% of total precipitation (Mesinger et al., 2006). This 15 feature has notable influence on hydrological and ecological system functioning in these cold environments (Beria et al., 2018). From a hydrological perspective, snowpack dynamics greatly influence water infiltration in soils, groundwater and aquifer replenishment, runoff production and water supplies to both natural and artificial water bodies during spring flood (Li et al., 2017;Barnhart et al., 2016;Berghuijs et al., 2014). From an ecological perspective, snowpack accumulation protects exposed plant tissues and organs against cold winds (Boivin and Bégin, 1997). Snow melt contributes to mitigating ::::::: mitigate the 20 negative impacts of droughts on tree growth (St. George et al., 2009), while affecting photosynthesis (Perkins and Swetnam, 1996;Peterson and Peterson, 1994). Snowpack dynamics also have the potential to alter heat fluxes, temperature and depth of freezing in soils, all of which can impact the timing of critical ecophysiological processes , driving ::: that :::: drive : growth in high latitude forest stands.
molecules, undetermined but expected in a range of 24-30 ‰ (DeNiro and Epstein, 1979;Farquhar et al., 1998). * is the equilibrium fractionation due to the change of phase of water from liquid to vapor at leaf temperature (fixed at 21.4 • C, the temperature threshold for maximum carbon assimilation), with a value of 9.65 ‰ (Helliker and Richter, 2008). k is the kinetic fractionation due to the diffusion of vapor into unsaturated air through the stomata and the leaf boundary layer, set to 26.5 ‰ in Farquhar et al. (1989) but we consider it undetermined because it can vary over larger ranges (Buhay et al., 1996). 120 h air is the relative humidity, which is estimated in MAIDENiso from the daily air temperature and the dew point temperature (Running et al., 1987). δ 18 O V and δ 18 O XW are the δ 18 O of vapor and xylem (source) water, respectively. δ 18 O V is calculated from the δ 18 O of precipitation (δ 18 O P ) and the fractionation due to the phase change from liquid water to vapor at mean air temperature, * Tair (Horita and Wesolowski, 1994): The δ 18 O TRC time series produced through Eq. (B1) are daily, while the δ 18 O TRC measured from tree rings is commonly annually resolved, or occasionally with intra-annual resolution (e.g. Szejner et al. (2018)). To produce a yearly record comparable with observations, the daily series are weighted with GPP. This assumes that allocation of carbon to the trunk is proportional to daily GPP (GPP d ): The original version of MAIDENiso (Gennaretti et al., 2017a;Danis et al., 2012) includes one snow layer, where snow accumulates and melts following changes in atmospheric temperature. This module was implemented to simulate snow reflectivity and thus changes in albedo within :: as ::: part :: of : the calculation of the energy budget. However, this is :::: was a side-calculation that did not interact with any of the other subsystems in MAIDENiso, and thus the accumulated and melted snow was not taken into account in the water balance calculation. All :: In ::::::: addition, ::: all : water pools and fluxes in the model were liquid , regardless 135 of temperature. In boreal climate, this previous version of MAIDENiso was thus unable to predict snow accumulation during winter, and therefore did not include different source water signatures due to snowfall instead of rainfall, the fractionation of δ 18 O due to snow sublimation , or the rapid melting of snow in spring. Therefore, this ::::::: previous version of the model simulated unrealistic soil moisture and hydrological outflux (drainage and runoff) , and too depleted values of δ 18 O in source water in spring.

New implementations in the model
The hydrology in the new version of MAIDENiso incorporates several pools of solid water: a canopy snow pool, a single snow layer on top of the soil and a pool of ice in each soil layer (Fig. 1). In addition, the snow layer is able to hold liquid water in its porous space, adding a new pool of liquid water. These water pools and the new water fluxes are shown in Fig. 2.
Input precipitation to the system is first partitioned into rainfall and snowfall driven by :::: based ::: on the average daily tempera-145 ture, following a linear partition between -2 • C and 4 • C (McCabe and Wolock, 2009). Following the same interception rule as liquid water, snow can be intercepted by the canopy and added to a ::: the canopy snow pool with a maximum capacity determined by LAI, from where it :::: snow : can sublimate. However, while the canopy water pool is always emptied :::::: always ::::::: becomes :::::: empty at the end of each day, the canopy snow pool does not. Canopy snow can still drip to the ground based on atmospheric temperature , following a drip model taken from the Community Land Model version 5 (CLM5) (Lawrence et al., 2019).

150
A single, uniform snow layer can cover the uppermost soil layer fully or partially, keeping track of snow thickness and the masses of snow and liquid water in the layer , ::: and : calculating the density of the layer dynamically. Freezing transfers water to snow mass without changing thickness, thus increasing density (with pure ice density as maximum), while melting and sublimation remove snow mass but keep density constant. The snow layer is forced to have a minimum thickness of 0.1 m (which is needed for numerical convergence), so the partial snow cover is decreased to avoid a thickness below this minimum 155 (i.e. partial cover of zero when no snow is present). While snowfall always accumulates ::::::: Snowfall :::::: always :::::::::: accumulate over the existing snow layer, increasing mass and thickness according to a temperature-variable density model of newly fallen snow (van Kampenhout et al., 2017). ::: In ::::::: contrast, the portion of rainfall that hits the snow layer is determined by the partial snow cover. Sublimation from the snow layer is calculated by modifying the version of the Penman-Monteith equation (Stigter et al., 2018) :: as :::::: follows: is the gradient of the saturation vapor pressure curve, γ (kPa • C −1 ) is the psychrometric constant (Loescher et al., 2009), R (MJ) is the net radiation over the snow surface, ρ air (kg m −1 ) is the air density, C P (MJ kg −1 ) is the specific heat of dry air, and r a (s m −1 ) is the aerodynamic resistance to water vapor transfer. r a typically depends on several factors, e.g. :::: such :: as wind (Blanken and Black, 2004). However, :::::: because : wind data is not :::::: usually available in tree-ring sites, which makes 165 formulas dependent on the wind unusable :::: some ::::::::::: assumptions :::: need ::: to :: be ::::: done :: to ::: use ::: the :::::::: equation ::::: above. Here, we assume a constant r a but optimize its value ::: that :: is ::::::::: optimized for each site using the available data of snowfall and snow pile's thickness and mass.
A pool of ice has been added to each soil layer. The pools of liquid water and solid water (ice) in each layer compete for the same porous space, so ::: thus : the ice content of a soil layer decreases its effective porosity. This decreases both the maximum 170 amount of liquid water the ::: that : a : layer can hold , and the hydraulic conductivity of the layer. Soil ice increases when soil temperature is below 0 • C , and decreases when soil temperature is above 0 • C.
To calculate the change in water phase from solid to liquid in :::: both the snow and soil layers, we have added a one dimensional ::::::::::::: one-dimensional : (vertical) thermal conduction model largely based on CLM5 (Lawrence et al., 2019). In this model, the system formed :::::::: composed : by the snow-soil layers is bounded at the top (as soil or snow) by the heat flux from the overlying atmosphere

175
, and at the bottom by a constant value representing the geothermal heat flux. The amount of water (or ice) that freezes (or melts) is calculated from the deficit (or excess) of energy to keep the temperature of the layer at 0 • C.
The new implementation of snow in MAIDENiso is now able to reproduce the dynamics of the snowpack, which influences :: is :::::::: connected :: to the rest of the components of the model (Fig. 1). The accumulation of the winter precipitation increases the amount of water available in the soil in spring, which in turn favors :::: may :::: favor : the onset of photosynthesis. A higher photosynthetic 180 activity results in more carbon assimilated by the canopy, potentially leading to a shorter budburst phase. A diagram showing the links between the different components in MAIDENiso is shown in Fig. 1.

Calibration of MAIDENiso
Different parameters that are species-dependent and site-dependent need to be defined before running MAIDENiso at a particular site. Most of these parameters can be obtained from direct observations at the studied site, such as the characteristics of the soil (composition and depth) , or the root-leaf proportions of the tree species. However, when ::::: When the values of the parameters are unknown, these are calibrated through a Bayesian optimization algorithm described in detail in Gennaretti et al. 200 (2017a). This optimization is based on Markov chain Monte Carlo (MCMC) sampling that retains combinations (blocks) of parameters that satisfy a condition, maximizing the coincidence between a series of observations and the equivalent products simulated by MAIDENiso. Here,we used 50 independent chains of parameter blocks and selected the most optimal block of parameters (called the "plausible block").
-GPP: 6 parameters (see appendix A). Calibrated by comparing observed and simulated daily GPP. Because the new and the original versions of MAIDENiso behave differently, an independent calibration of the parameters is needed to run each of them. Note that the original version of MAIDENiso does not need to be calibrated for the snow parameters because it does not include a snow module.
Some parameters can influence more than one process indirectly, e.g. the snowpile affects source water and therefore 215 δ 18 O TRC , or the GPP parameters control the amount of carbon produced, which in turn affects δ 18 O TRC . To avoid that the calibration of some processes affects parameters that are already calibrated, the parameter sets need to be calibrated in a specific order: Snow :::: snow : first, GPP second, and lastly δ 18 O TRC .

Study sites and input meteorological data
The studied :::: study : tree-ring sites are located in Tungsten, Northwest Territories, Yukon border ( MAIDENiso requires daily meteorological inputs for a continuous period of time overlapping with the period of available observations. Daily CO 2 data were obtained from the Mauna Loa Observatory observations (Keeling et al., 1976) corrected with the CarbonTracker measurement and modeling system (Peters et al., 2007).
In-situ snowpile data are needed to test the predictive skills of MAIDENiso to simulate the snowpile. ::: The : Snow Water Equivalent (SWE) is the ideal snowpile data to use because addition (from precipitation) and removal (from sublimation 270 and melting) of snow to or from the snowpile is calculated in units of mass. Alternatively, snow depth (SNDP) data, most commonly available, can be used as well to compare with observations but requires knowledge about snow density. SWE field measurements were only available for the Caniapiscau site at discrete (biweekly) intervals during Winter and early Spring ::::: winter :::: and :::: early :::::: spring between 1971-1993 (data provided by Hydro-Québec, personal communication). Therefore, in order to make the results from both sites comparable, we used SNDP data from NARR to calibrate the snowpile at Tungsten and

305
(2017), we compared the predicted δ 18 O TRC from the reference simulations with those obtained from two experiments. First, to isolate the contribution of the source water on δ 18 O TRC , we set the relative humidity (h air ) and δ 18 O V constant using the average values of h air and δ 18 O V obtained from the reference simulations. Second, to isolate the contribution of the isotopic enrichment of the leaf water during transpiration on δ 18 O TRC , we set :::: δ 18 O :: in :::::: xylem ::::: water : (δ 18 O XW : ) : constant using the average value of the reference simulation. We then compared the reference and experimental simulations using the coefficient 310 of determination (R 2 ).

345
The calibration process converged and constrained the values of r a well (resistance to water vapor transfer) ::: well : at both sites, as shown in Table 4. The values obtained for both sites were quite different, equal to 47.88 s m −1 in Caniapiscau but to 87.35 s m −1 in Tungsten. This potentially indicates that the physical factors that are neglected when assuming a constant ::: We ::::::: obtained : a ::::: value :: of : r a (such as wind speed) are quite different between the two sites. The higher values at Tungsten indicates that snow sublimates at a slower rate, likely because the average wind speed during winter at this site is smaller ::::: almost ::::: twice 350 :: as :::: large :: at :::::::: Tungsten :::::::::::: (87.35 s m −1 ) than at Caniapiscau ::::::::::: (47.88 s m −1 ).

GPP calibration
GPP was calibrated twice at each station: first for the original version of MAIDENiso and second for the new version with 370 snow. Both versions were able to predict GPP observations , and in a similar way for both the Uaf and the EOBS stations (Fig. 5). The observed and simulated GPP were in good agreement regarding the timing (onset and offset) of the yearly peaks.
The maximum of these peaks was higher for the observations, but this is due to exceptional days of very high observed GPP.
Overall, the average GPP during the growing season was well reproduced at both GPP stations.
The distribution of the optimized parameters controlling δ 18 O TRC can help to understand ::::: assess how the model compensated for the presence or absence of snow (Table 4 and Fig. 7). For instance, a lower value for the dampening factor than the expected range of f 0 = 0.4 − 0.5 obtained from previous studies (Roden and Ehleringer, 2000;Saurer et al., 1997) would indicate that the model does not simulate the δ 18 O XW signal correctly. Higher values of the biochemical fractionation 0 , or the kinetic fractionation k , in one version of MAIDENiso than the other, would indicate that the model is compensating for comparatively lower values of δ 18 O XW and/or δ 18 O V in Eqs. (B1,B2). In our simulations, adding snow induced an increase in the dampening factor f 0 at both sites, suggesting that the signal of the source water on δ 18 O TRC was stronger than without considering snow.
The existence of the snow layer induced a higher δ 18 O at the leaf level, especially in spring, due to the higher δ 18 O XW level.
In this study, we implemented a new snow module in MAIDENiso to simulate snowpack dynamics and improve the model representation of the soil hydrology and the isotopic fractionation of oxygen in water and tree-ring cellulose. In the following 435 paragraphs we discuss the impacts of the snow module addition on the different components of the model (i.e. the hydrological, photosynthetic and δ 18 O modules), addressing the skills and limitations of our approach. We also discuss the implications of our new snow module implementation in MAIDENiso for future studies.

No effects of snow on photosynthesis
The calibration process yielded two different sets of parameters when using the two versions of the :::::::::: MAIDENiso : model but resulted in similar predicted GPP values. However, for both versions of MAIDENiso, we obtained :: Th ::::::::: parameters :::::::::: controlling :::: GPP ::: wee :::::::: obtained :::::: showed : very similar posterior distributions and values in the plausible block for the parameters controlling 460 GPP (Table A1, Figs A1, A2). These similarities indicate that, at our study sites, photosynthesis is ::: was : not very sensitive to the availability of additional water from snow melt, suggesting ::: that : radial tree growth is ::: was : not limited by water availability. This was expected, because in :: is :: in ::::::::: agreement :::: with ::::::: previous :::::: studies ::::::: showing :::: that :: in high latitudes soil humidity is not :::: often : a major constraint on tree growth and trees are usually mostly sensitive to temperature (Boisvenue and Running, 2006; D'Orangeville from snowmelt (Du et al., 2014). Because GPP is not affected by the snow module, our study sites are ideal to investigate the effects of snow on δ 18 O TRC variations because it allows us to discard GPP as a possible cause for the differences observed between the two model versions.
(2) could be calibrated to compensate for the absence of snow. These findings indicate that the addition of snow allows the model to increase the contribution of the source water to δ 18 O TRC .
The addition of the snow module to MAIDENiso therefore frees the calibration process from having to overcompensate for the artificially depleted values of δ 18 O XW ::::: values : during the growing season (Fig. 8b). As our results have shown, this significantly increased the correlation between the observed and simulated δ 18 O TRC compared to the version without snow 515 (Fig 6; r = 0.52 for Caniapiscau and r = 0.57 for Tungsten, p < 0.01 versus non-significant, respectively). The improvement of the predictive skill of the model with the snow module reflects the influence of winter precipitation on physiological processes.
Without snow, all winter precipitation passes through and out of the hydrological system without affecting the trees. In contrast, including snow allows winter precipitation to affect δ 18 O TRC indirectly through the source water.
Overall, the improvements found in the δ 18 O TRC simulations at both sites indicate that snow plays a critical role in δ 18 O of 520 the source water and thus on the final signature of δ 18 O TRC . Even if the addition of snow would not had resulted in a significant improvement of the correlation between the simulated and observed δ 18 O TRC , accounting for snow-related processes along the mechanistic chain is necessary for the application of a process-based model in an environment where snow is present. Processbased models are useful to understand complex processes, and while they may not necessarily produce better simulations (closer to observations) than response functions, they can be used under different conditions from those of the data used for 525 their calibration :::::::: calibrated ::::: under :::::::: favorable :::::::: conditions :::: and :::: then :::: used :: for :::::::: different :::::: datasets : . The incorporation of the snow module in MAIDENiso is therefore required for predicting tree-ring isotopic composition in forests located in cold environments where snow is present.

Implications for future studies
Based on our results and comparison with other studies, we can conclude that the snow module predicts ::::::: predicted : more realistic 530 and robust fluxes of water within the soil-plant-atmosphere continuum. The improvement of MAIDENiso to disentangle the contribution from the source (xylem) water and the δ 18 O leaf enrichment signal on δ 18 O TRC can help to track the origin of the isotopic signal and eventually improve the interpretations of the climate signal recorded in the tree rings. This is important because tree-ring isotopes are important climate proxies (Cernusak and English, 2015). The inclusion of the new snow module in the model can provide a more accurate representation of the physical and physiological processes taking place than in earlier 535 studies which ignored ::: that ::: did ::: not :::: take :::: into ::::::: account the additional effects of snowpack dynamics on δ 18 O TRC , e.g. Lavergne et al. (2017). Now, MAIDENiso can simulate more reliably the interaction :::::: reliable :::::::::: interactions between the coupled water and carbon cycles and the tree physiology ::: tree ::::::::::: physiological ::::::::::: mechanisms : in cold environments. The ::: Our ::::::: findings :::: will ::::::::: contribute :: to ::::: reduce : uncertainties in the predictions of the response of forest productivity to hydrological changesmay also be reduced, leading to better forward predictions that can eventually be used to reconstruct seasonal and long-term hydroclimatic variations.

Conclusions
In this paper we presented the new snow module incorporated into MAIDENiso, which consists of new hydrological calcula-555 tions of snow dynamics and a thermal module. Our findings :::::: results show how this snow module improves the simulation of outputs associated with the hydrological cycle at cold and high latitudes :::::: latitude :::: sites : without affecting simulations from the carbon cycle component. These findings were expected because GPP and tree-ring growth at the studied boreal high-latitude sites are not constrained by soil moisture availability but by surface air temperature and light (Jarvis and Linder, 2000). The MAIDENiso. Based on the development presented here, the potential for the application of MAIDENiso is notably increased.
Code and data availability. Appendix A: Photosynthesis model GPP (g C m −2 day −1 ) in MAIDENiso derives from a coupled photosynthesis-stomatal conductance system. The leaf photosynthesis is modelled following Farquhar et al. (1980), scaled to the canopy followingDe Pury and Farquhar (1997) as explained in Misson (2004). Daily Vcmax (Vcmax i ) is modelled as: The parameter Vmax determines how daytime temperature Tday controls the maximum carboxylation rate at day i. Because there was no explicitly known mechanistic formula relating Vcmax and Tday, three parameters were introduced to control this relationship in a non-linear way, i.e. Vmax, Vb and Vip. These parameters control the asymptote, the slope and the inflection point of Vcmax, respectively, and have to be calibrated.

575
The stomatal conductance for carbon (µmol m −2 s −1 ) is calculated using the Leuning et al. (1995) model, modified by Gea-Izquierdo et al. (2015) to incorporate soil water stress: where g 0 = 0µmol m −2 s −1 and g 1 = 10µmol m −2 s −1 are fitted parameters representing the residual conductance as the net assimilation rate (A n ) approaches zero and the slope of the function, respectively. P atm is the atmospheric pressure (Pa).

580
C a is the atmospheric CO 2 pressure (Pa). Γ * is the CO 2 compensation point in the absence of dark respiration (Pa), which is calculated following Bernacchi et al. (2001). VPD is the vapor pressure deficit (kPa), and VPD 0 is an empirically fitted parameter representing the sensitivity of stomata to changes in VPD (usually around 15 kPa; Knauer et al. (2015)). θ g is the empirical soil water stress factor, a non-linear function ranging between 0 when the soil is too dry for the roots and 1 in absence of water stress: The water stress level depends on the soil water content (SWC, mm), but the current version of MAIDENiso lacks a mechanistic model to explain the relationship between soil water content and water stress. For this reason, this relation is modelled as a logistic function, introducing the calibration parameters soilb and soilip as the slope and the inflexion point of θg.
Finally, there is a time lag between the recovery of photosynthesis and the temperature increase in spring that is taken into 590 account by the model. This is done by replacing Tday in Eq. (A1) by the temperature transformation S, defined as: where τ is a parameter representing the number of days needed by the tree to adapt the photosynthesis to changing temperatures.

595
These parameters were calibrated at the two eddy covariance flux stations described in section 2.4, for both versions of the model. These parameters, their prior distributions and their posterior distributions are shown in Table A1. For better visualization, the probability distribution function (pdf) of the posterior distributions of the GPP parameters are also shown in Figs. A1 and A2.