Articles | Volume 15, issue 5
Model evaluation paper
09 Mar 2022
Model evaluation paper |  | 09 Mar 2022

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

Ignacio Hermoso de Mendoza, Etienne Boucher, Fabio Gennaretti, Aliénor Lavergne, Robert Field, and Laia Andreu-Hayles

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 implementations in two sites in eastern and western 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 cellulose, leading to more accurate estimates of cellulose isotopic composition. Future work may include the development of inverse modelling with this new version of MAIDENiso to produce robust reconstructions of the hydrological cycle and isotope processes in cold environments.

1 Introduction

In boreal regions of Canada and Alaska, snow represents about 30 %–50 % of total precipitation (Mesinger et al.2006). This feature has a 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égin1997). Snowmelt contributes to mitigate the negative impacts of droughts on tree growth (St. George et al.2009), while affecting photosynthesis (Perkins and Swetnam1996; Peterson and Peterson1994). 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 that drive growth in high-latitude forest stands.

For decades, tree-ring proxies such as ring widths (Nicault et al.2015; Ols et al.2018), wood density (Boucher et al.2017) or stable isotope ratios of tree-ring cellulose (Naulier et al.2015b, 2014; Porter et al.2014) have been used to track inter-annual changes in forest response to climate variability. Most studies emphasized the dominant role of summer temperatures on key ecophysiological processes controlling proxy formation. This has helped to clarify the response mechanisms of the boreal forest to growing-season temperatures (Gennaretti et al.2017a) and enabled long, millennial summer temperature reconstructions to be produced in this region (Gennaretti et al.2017a, b; Naulier et al.2015a). However, despite their ecological and hydrological significance, snow-related processes were rarely taken into account in these tree-ring studies (Coulthard et al.2021; Woodhouse2003; Huang et al.2019; Zhu et al.2021). Consequently, the impacts of these changes in snow cover properties (Meredith et al.2019) on vegetation growth and ecophysiological response remain highly uncertain.

Predicting the effect of snow dynamics on tree growth is a complex task as both phenomena occur in distinct seasons (Coulthard et al.2021). Inter-seasonal heat and moisture fluxes attributable to snow need to be accounted for in order to accurately model the impact of snow on tree-ring formation. The timing and magnitude of these transfers, however, result from a complex interplay between snowpack properties (snow depth, density and water content) and processes that control snow accumulation and melt (precipitation, sublimation, redistribution by wind, rain-on-snow events, among others) (Rutter et al.2009). These transfers also modify the isotopic signature of the water used by trees. Indeed, snow is more depleted in the lighter isotope 18O than rainfall (Kurita et al.2004), but sublimation-driven enrichment of snow may also change the isotopic composition of the source water used by trees. Ultimately, this should be recorded in the δ18O of tree-ring cellulose (Beria et al.2018). Correlation-based tree-ring analyses based on statistical relationships cannot take into account this mechanistic level of complexity; thus, there is a need to explicitly integrate snow dynamics in forest growth models.

Process-based models developed for simulating tree growth are important tools to study the relationship between climate and tree-ring proxies (Guiot et al.2014). These models are driven by meteorological and environmental variables and integrate a wide number of equations that represent state-of-the-art knowledge on how physical and ecophysiological processes determine tree response to climate variability. A number of process-based models have been developed over the years, such as the Vaganov–Shashkin (VS) model (Fritts et al.1991; Vaganov et al.2006; Shishov et al.2016), MAIDEN (Misson2004), StandLeap (Girardin et al.2008), CAMBIUM (Drew et al.2010), ECOPHYS (Hölttä et al.2010), Biome3 (Rathgeber et al.2003) or the T model (Li et al.2014). Despite the importance of snow for tree growth, most process-based models do not include a snow module, mostly because they were not designed to be used in boreal and alpine environments or even in mid-latitude temperate forests where snow accumulates during winter. Among the previously mentioned models, exceptions are the Vaganov–Shashkin (VS) model (Shishov et al.2016) and the Biome3 model (Rathgeber et al.2003), which incorporate basic models of snow accumulation and melt driven by air temperature but do not consider processes such as sublimation, energy balance or stable isotope fractionation of water isotopes during the cold season. Among the available models, MAIDEN (Misson2004) was specifically designed to improve the interpretation of tree-ring proxies based on our knowledge about ecophysiological processes and relationships between climate and tree growth. MAIDEN simulates the water and carbon fluxes exchanged between forests and the atmosphere, including the influence of phenology on the production and allocation of carbon to different parts of the tree. Because it requires a very limited number of meteorological inputs, the application of the model is possible in regions where data are scarce. The isotope-enabled version, MAIDENiso (Danis et al.2012), incorporates calculations of the stable isotopic composition of oxygen (δ18O) and carbon (δ13C) in the different components of the tree. MAIDEN was originally created for tree species in Mediterranean climates, and it has been optimized for Quercus petraea (Matt.) Liebl. and 12 Mediterranean species (Misson2004; Gaucherel et al.2008; Boucher et al.2014; Gea-Izquierdo et al.2015). Since then, the phenology and physiological processes have been adapted to simulate tree radial growth in boreal northeastern American forests (Gennaretti et al.2017a) and used to simulate tree-ring cellulose δ18O in boreal and temperate forests of eastern Canada and southern South America (Lavergne et al.2017). MAIDENiso provides two main advantages over other process-based models. (1) The outputs are directly comparable to tree-ring proxies. (2) It is an isotope-enabled model, allowing users to track down the origin of the climate signal recorded therein. However, the use of MAIDENiso in high-latitude forests has been limited by the fact that its hydrological cycle was never adapted to boreal conditions and the lack of an adequate representation of snow dynamics.

Here, we incorporate a new snow module in MAIDENiso and test the simulated data against real observations. This module is driven by a new thermal conduction model to improve the simulations when the model is used in cold environments where snow is present. This snow module allows MAIDENiso to reproduce the basic dynamics of the snowpack, targeting a more realistic water balance and water isotope fractionation sequence by representing the δ18O signal of snowfall, the sublimative fractionation at the snow surface and its final imprint in tree-ring cellulose (TRC). Despite this added complexity on processes, the snow model can work with the same small number of environmental variables that MAIDENiso currently requires. In this study, we evaluate the impact of the new snow module on the simulation of soil moisture, water outflux and the δ18O signal in soil and TRC in two forest sites in Canada: a black spruce (Picea mariana (Mill.) B.S.P.) forest in the Caniapiscau basin (Quebec) and a white spruce (Picea glauca (Moench) Voss) forest in Tungsten (Northwest Territories).

2 Materials and methods

2.1 MAIDENiso model

2.1.1 Original model

MAIDENiso (Misson2004; Danis et al.2012; Gea-Izquierdo et al.2015; Gennaretti et al.2017a) simulates the mechanical and physiological processes of a tree and its immediate environment. The model requires daily meteorological inputs of maximum and minimum temperature, precipitation, and atmospheric CO2 concentration (optional inputs are relative humidity, radiation, wind speed and atmospheric δ13C). MAIDENiso simulates gross primary production (GPP) and carbon allocation on a daily basis based on inputs of meteorological and tree phenological data. Carbon is allocated explicitly to several pools (leaves, roots, stem and a carbon reservoir) using mechanistic rules dependent on phenology. A diagram of the model is shown in Fig. 1 with the original components of the model depicted in black. These original components include the photosynthesis module and the isotopic module, which are described in Appendix A and Appendix B, respectively.

Figure 1Diagram showing the main features in the new version of MAIDENiso, with old components and fluxes in black and new ones in blue for snow/ice and in red for the thermal module. Processes are in italics, boxes are carbon and water pools, broken lines are links between processes, and solid lines are carbon and water fluxes. Figure was modified from Misson (2004).

MAIDENiso simulates the hydrological processes in the immediate environment around the tree: at canopy (interception and canopy evaporation), ground surface (infiltration, evaporation and runoff) and underground (hydraulic transfers and root absorption) levels. These processes are modelled through a series of water pools and fluxes (Fig. 2). For instance, the canopy can intercept a portion of the precipitation water up to a maximum determined by the leaf area index (LAI), which can be evaporated or dripped to the ground overnight. The surface of the soil cannot hold any stagnant water, so daily incoming water from throughfall infiltrates the soil (up to a maximum determined by soil properties) or exits the system as runoff. The soil consists of four layers of distinct thickness, with porosity and hydraulic conductivity determined by the composition of the soil, and water moves between these layers following Darcy's law. Soil water is replenished through infiltration and depleted by root absorption for transpiration (at all layers), soil evaporation (at the upper layer) and drainage (at the bottom layer).

Figure 2The hydrological system in the new version of MAIDENiso. Pools (flasks) and fluxes (arrows) are shown for liquid water in dark blue and for snow/ice in light blue.


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; thus, changes in albedo as part of the calculation of the energy budget. However, this 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. In addition, all water pools and fluxes in the model were liquid regardless of temperature. In boreal climate, this previous version of MAIDENiso was thus unable to predict snow accumulation during winter; therefore, it did not include different source water signatures due to snowfall instead of rainfall, the fractionation of δ18O 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 values of δ18O in source water in spring that are too depleted.

2.1.2 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 based on the average daily temperature, following a linear partition between 2 and 4 C (McCabe and Wolock2009). Following the same interception rule as liquid water, snow can be intercepted by the canopy and added to the canopy snow pool with a maximum capacity determined by LAI, from where snow can sublimate. However, while the canopy water pool 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).

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 (i.e. partial cover of zero when no snow is present). Snowfall always accumulates 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:

(1) λ E pot,snow = Δ R + ρ air C P δ atm / r a Δ + γ ,

where Δ (kPaC-1) is the gradient of the saturation vapour pressure curve, γ (kPaC-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, CP (MJ kg−1) is the specific heat of dry air and ra (s m−1) is the aerodynamic resistance to water vapour transfer. ra typically depends on several factors, such as wind (Blanken and Black2004). However, because wind data are not usually available in tree-ring sites, some assumptions need to be made to use the equation above. Here, we assume a constant ra 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; thus, the ice content of a soil layer decreases its effective porosity. This decreases both the maximum amount of liquid water 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 (vertical) thermal conduction model largely based on CLM5 (Lawrence et al.2019). In this model, the system composed of the snow–soil layers is bounded at the top (as soil or snow) by the heat flux from the overlying atmosphere 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 is connected to the rest of the components of the model (Fig. 1). The accumulation of winter precipitation increases the amount of water available in the soil in spring, which in turn may favour the onset of photosynthesis. A higher photosynthetic 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.

The new MAIDENiso version also includes new isotopic fractionation processes for the sublimative fluxes and for the phase changes between liquid water and ice. In cold regions where snowfall is a considerable portion of the yearly precipitation, fractionation from snow sublimation is expected to produce a significant enrichment of the δ18O isotopes in the snow layer, which after melting may be incorporated into the soil water, and ultimately reflected in TRC.

Given the already high number of parameters in MAIDENiso (121 in the new version with snow, 117 in the previous version), one of our goals during the development of the new snow module has been to keep the number of new free parameters to the minimum possible. Despite the complexity and the new processes incorporated into the model, the new snow module only added four new parameters to MAIDENiso, which are listed in Table 1. Three of them are site parameters (determined externally to MAIDENiso) that correspond to a linear regression model of precipitation δ18O (more information in Sect. 2.3). Therefore, we only added a single free parameter that requires calibration: the resistance to vapour transfer ra. This parameter controls snow sublimation, which fundamentally depends on wind speed and therefore varies considerably between sites, requiring independent calibration at each site that, as explained above, we computed using observations of snowfall and snow pile's thickness and mass, because wind data are not available.

Table 1New parameters introduced to MAIDENiso in the new version.

Download Print Version | Download XLSX

2.2 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. When the values of the parameters are unknown, these are calibrated through a Bayesian optimization algorithm described in detail in Gennaretti et al. (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”).

The series of observations used in the MCMC consist of observed time series of snow (depth or mass of the snow pile), GPP and δ18O in tree-ring cellulose (δ18OTRC). The parameters to be determined via MCMC for each component of the model are the following:

  • Snow pile: one parameter, ra from Eq. (1). Calibrated by comparing observed and simulated daily snow depth (SNDP).

  • GPP: six parameters (see Appendix A). Calibrated by comparing observed and simulated daily GPP.

  • δ18OTRC: three parameters (see Appendix B); that is, f0, ϵ0 and ϵk in Eqs. (B1) and (B2). Calibrated by comparing observed and simulated yearly δ18OTRC.

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 as it does not include a snow module.

Some parameters can influence more than one process indirectly; for example, the snow pile affects source water and therefore δ18OTRC, or the GPP parameters control the amount of carbon produced, which in turn affects δ18OTRC. 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 first, GPP second and lastly δ18OTRC.

2.3 Study sites and input meteorological data

The tree-ring study sites are located in Tungsten, Northwest Territories, Yukon border (61.98 N, 128.25 W; 1145 m a.s.l.), and in the Caniapiscau basin, Quebec (54.86 N, 69.72 W; 530 m a.s.l.).

MAIDENiso requires daily meteorological inputs for a continuous period of time overlapping with the period of available observations. Daily CO2 data were obtained from the Mauna Loa Observatory observations (Keeling et al.1976) corrected with the CarbonTracker measurement and modelling system (Peters et al.2007).

The closest meteorological stations to the study sites were located 100 km away from Tungsten and 186 km from Caniapiscau. Therefore, temperature and precipitation data were taken from the NARR (North American Reanalysis) dataset (Mesinger et al.2006) at the coordinates of the studied sites. The NARR has a spatial resolution of 32.5 km×32.5 km and spans the period 1979–2013. Meteorological inputs were also needed at two additional sites to calibrate GPP, which are described in detail in Sect. 2.4. During the period 1979–2013 based on NARR, average summer temperatures (June–July–August) ranged 5.4–12.6 C in Tungsten and 8.2–16.2 C in Caniapiscau. Yearly temperatures at Tungsten are stable during the whole period, while at Caniapiscau average temperatures steadily rise after 1990 by 0.1 Cyr-1. Average yearly precipitation values during this period were 584 and 796 mm in Tungsten and Caniapiscau, respectively. In Tungsten, 56 % of the yearly precipitation was snowfall, while in the warmer site of Caniapiscau snow was only 45 %.

MAIDENiso also needs information about δ18OP. Two different approaches can be used to infer δ18OP. The first and most direct way is to use daily values of δ18OP as another meteorological input. However, these values are often not available. The second approach, which we used here, is to use precipitation (P, mm), air temperature (Tair, C) and δ18OP from an observed dataset to obtain a linear regression model for daily values of δ18OP based on temperature and precipitation:

(2) δ 18 O P = a T air + b P + c .

This approach has the advantage that, once the model is obtained, the parameters can be used with a different dataset of air temperature and precipitation to obtain the corresponding δ18OP values. In this paper, we calibrated this regression model using meteorological data from the gridded dataset IsoGSM (Yoshimura et al.2008) from the grid points that contain the coordinates of the Tungsten and Caniapiscau sites. We discarded the direct use of the IsoGSM meteorological and δ18OP data for MAIDENiso as the precipitation amounts derived from IsoGSM were too low compared to the amounts from observations recorded in meteorological stations nearby the study sites. However, the IsoGSM meteorological data were still useful to obtain the parameters for our regression model. We obtain different equations for liquid (rainfall: arain, brain, crain) and solid precipitation (snowfall: asnow, bsnow, csnow) using separately data corresponding to temperatures below 4 C for snowfall and higher than 2 C for rainfall (see Table 2).

Table 2Parameters obtained for the regression models for snow and rainfall using the IsoGSM dataset at the Tungsten and Caniapiscau sites. The thresholds for the null hypothesis are NNSE=0.5 and KGE=0.41.

Download Print Version | Download XLSX

2.4 Tree-ring δ18O, GPP and snow data

We used published δ18OTRC chronologies for Tungsten (Field et al.2021) and Caniapiscau (Nicault et al.2014). These chronologies span between 1900–2003 for Tungsten and 1948–2013 for Caniapiscau; however, for this study we used the isotopic records for periods that overlap with the NARR meteorology: 1979–2003 for Tungsten and 1979–2013 for Caniapiscau.

We used GPP data available from the closest eddy covariance flux stations to estimate the parameters controlling GPP, assuming that the obtained parameters were similar at the studied sites. To calibrate GPP in Tungsten, we used the University of Alaska Fairbanks (Uaf) station from the Ameriflux network (64.87 N, 147.85 W; data period 2003–2018; Ueyama et al.2021) at 1023 km from our study site. For Caniapiscau, we obtained daily GPP data from an eddy covariance station located in a mature black spruce forest in northern Quebec (“Quebec Eastern Old Black Spruce station” – EOBS; 49.69 N, 74.34 W; (last access: 26 January 2016); data period 2003–2010; Bergeron et al.2007) at 650 km from our study site. Although these eddy covariance flux stations are geographically distant from our study sites, they provide GPP data for the same tree species in our sites. Because the parameters used to calibrate GPP are more related to species-specific traits (Gennaretti et al.2017a) than environmental conditions at a given site, it is a reasonable assumption to calibrate GPP at these stations and use the obtained GPP parameters in our study sites. MAIDENiso simulated GPP at both stations, using the following meteorological inputs. For the Uaf station, we used the meteorological inputs available at the station. For the EOBS site, the meteorological inputs were taken from the gridded interpolated Canadian database of daily minimum–maximum temperature and precipitation for 1950–2015 (Hutchinson et al.2009), used in Gennaretti et al. (2017a).

In situ snow-pile data are needed to test the predictive skills of MAIDENiso to simulate the snow pile. The snow water equivalent (SWE) data are the ideal snow-pile data to use, because addition (from precipitation) and removal (from sublimation and melting) of snow to or from the snow pile 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 between 1971–1993 (data provided by Hydro-Québec, personal communication, 2019). Therefore, in order to make the results from both sites comparable, we used SNDP data to calibrate the snow pile at Tungsten and Caniapiscau and only used the SWE measurements at Caniapiscau to validate the simulations. The SNDP data (1979–2013) were extracted from NARR for Caniapiscau and from observations of a meteorological station for Tungsten.

2.5 Model evaluation and experiments

To evaluate the agreement between observed and simulated δ18OTRC for the two versions of MAIDENiso, we calculated the Pearson correlation coefficient and associated p values (p<0.05 were considered significant). To determine that the simulated δ18O at the leaf and cellulose level were statistically different, we used the Welch t test (which tests the null hypothesis that the difference between the means of two curves is zero).

While MAIDENiso does not calculate river discharge as an output (which would be possible through the implementation of a routing model), water discharge (water leaving the system in liquid form) can be calculated as the addition of runoff (water overflowing the infiltration capacity of the soil) and drainage (downwards water flux from the lowest soil layer) and be compared to measurements of river discharge. For this study, these comparisons were only done for Caniapiscau due to the availability of river discharge observations. To evaluate the coincidence between the observed and simulated water discharge, we used the Nash–Sutcliffe model efficiency (NSE) coefficient (Nash and Sutcliffe1970), which is equivalent to a coefficient of determination:

(3) NSE = 1 - t ( Q sim t - Q obs t ) 2 t ( Q obs t - Q obs ) 2 ,

where Qsimt and Qobst are the simulated and observed discharge at time t, respectively. The NSE ranges between −∞ and +1. To facilitate the interpretation of NSE, we rescaled the NSE within the range of (0,1) with the normalized Nash–Sutcliffe efficiency (NNSE) coefficient (Nossent and Bauwens2011):

(4) NNSE = 1 2 - NSE .

Another useful metric is the Kling–Gupta efficiency (KGE) (Gupta et al.2009), which addresses several shortcomings of the NSE and is increasingly used for model calibration and evaluation:

(5) KGE = 1 - ( r - 1 ) 2 - σ sim σ obs - 1 2 + μ sim μ obs - 1 2 ,

where r is the linear correlation between observations and simulations, σ is the standard deviation and μ is the mean.

Values of NSE>0 (NNSE>0.5) are typically used as the benchmark to establish a model as a “good” model (NSE=0 indicates that the model is a predictor as good as the mean of the observations). The KGE equivalent values to consider that a model is skilful are KGE>1-2-0.41 (Knoben et al.2019).

To estimate the effect of the new snow module on predictions of δ18OTRC, we compared the parameters influencing δ18OTRC obtained by independent calibrations.

We also investigated the relative contributions to δ18OTRC of the source (xylem) water and of the fractionation processes during transpiration in the leaf (see the Eqs. B1 and B2 in Appendix B). Using the same approach as in Lavergne et al. (2017), we compared the predicted δ18OTRC from the reference simulations with those obtained from two experiments. First, to isolate the contribution of the source water on δ18OTRC, we set the relative humidity (hair) and δ18OV constant using the average values of hair and δ18OV obtained from the reference simulations. Second, to isolate the contribution of the isotopic enrichment of the leaf water during transpiration on δ18OTRC, we set δ18O in xylem water (δ18OXW) constant using the average value of the reference simulation. We then compared the reference and experimental simulations using the coefficient of determination (R2).

2.6 Validation of the snow model

The snow module was validated using a split-sample approach. At both of our study sites, we divided the period of SNDP observations (1979–2013) into half-periods: 1979–1996 and 1997–2013. We then used the SNDP observations in each half-period to calibrate the snow module, following the same procedure described in Sect. 2.2 for the whole period of observations. Using the snow parameter obtained from each half-period, we simulated the whole period and compared the simulated snow pile with the snow observations using the NNSE and KGE metrics. To account for gaps in the observational record, we performed this comparison on the mean annual cycle of the signal (calculated by averaging the series for every day of the year) which was also smoothed with a 30 d spline.

3 Results

3.1 Model validation

Using the split-sample calibration method, MAIDENiso was able to simulate SNDP data for the full period that compared well with the SNDP observations at each site, as indicated by the two NNSE and KGE metrics (Table 3). The obtained values were well above the thresholds that establish the model as a better predictor than the mean of the observations, with thresholds of NNSE=0.5 and KGE=0.41. In addition, the calibrations using half-periods produced similar results to the calibration using the full period, with the first half-period (1979–1996) producing slightly higher values than the second (1997–2013). The most notorious difference was the KGE obtained for the second half-period in Caniapiscau, which showed lower values (0.49) than the first half-period (0.69) and the full period (0.66), but it was still within the values needed to be considered a good model result.

Table 3Normalized Nash–Sutcliffe efficiency (NNSE) and Kling–Gupta efficiency (KGE) for the mean annual cycles of SNDP simulated by calibrating the model with the full period and the two half-periods, at the Tungsten and Caniapiscau sites.

Download Print Version | Download XLSX

3.2 Snow calibration and impact on hydrology

MAIDENiso simulated SNDP and snow density using NARR meteorological data at our two sites, which we compared with the SNDP product from the NARR dataset. To compare with the direct observations of SWE at the Caniapiscau site, we used the simulated snow density to transform the SNDP (both that simulated by MAIDENiso and the NARR product) into SWE, which are shown in Fig. 3.

Figure 3Snow water equivalent (SWE) averaged over 1979–1997 at the (a) Tungsten and (b) Caniapiscau sites, extracted from the NARR data (black) and simulated by MAIDENiso (red) using NARR meteorology. The solid line indicates the average of the same day of the year (DOY) during this period, and the shadows indicate the 2σ variability. Direct observations of SWE at Caniapiscau, taken at discrete intervals, are shown as blue dots.


The MAIDENiso simulations reproduced the temporal change of the snow pile and showed a similar pattern to the real SWE observations collected at the Caniapiscau site (Fig. 3b). In contrast, the NARR-based SWE estimates showed higher values of the snowpack during the winter months and offsets in the timing of snow accumulation and melting. The discrepancies between the NARR-based SWE data and the MAIDENiso SWE simulations could arise from a mismatch between the NARR meteorology (used to drive the model) and the NARR's snow-pile data, where, according to the available documentation, the latter was artificially increased to match other sources (Mesinger et al.2006). Therefore, our SWE simulations made by MAIDENiso using as inputs NARR meteorological data were in better agreement with the direct observations of SWE at the Caniapiscau site than the SWE data obtained directly from the NARR dataset.

The calibration process converged and constrained the values of ra (resistance to water vapour transfer) well at both sites, as shown in Table 4. We obtained a value of ra almost twice as large at Tungsten (87.35 s m−1) than at Caniapiscau (47.88 s m−1).

Table 4Calibration parameters for the snow pile (for the version with snow exclusively) and δ18OTRC. Parameters, units, prior range, and posterior range (with parameter value in the plausible block) for both MAIDENiso versions and both sites.

n/a – not applicable

Download Print Version | Download XLSX

Figure 4 shows the combined runoff and drainage simulated by MAIDENiso with and without including the snow module and the observations of river discharge in the Caniapiscau basin (scaled by the area of the basin). The observations showed a peak in river discharge between May and July, corresponding to the melting of the snow accumulated during winter. In the simulations computed using the original version of MAIDENiso (without snow module), the outflux of the model resembled the pattern of precipitation during the year, because all incoming precipitation was considered liquid and did not show any peak. Conversely, the simulations produced using our new MAIDENiso version (with snow module) did not have any outflux during winter, when all water is in solid state, and they reproduced more accurately the peak of water outflux during spring melting. Overall, the timing of the spring discharge was well reproduced by the MAIDENiso version with snow, while the original version was unable to simulate this peak. This improvement was confirmed by the NNSE between the observed and modelled river discharge at Caniapiscau, which was lower for the model without snow (NNSE=0.1) than for the new model with snow (NNSE=0.45). The KGE for the observed and modelled river discharge at Caniapiscau improved only from 1.04 (without snow) to 0.57 (with snow), which indicates that our modelled river discharge (instant additions of drainage and runoff from all of the basin) can still be further improved.

Figure 4Water outflux (drainage + runoff) simulated for 1979–1997 in the Caniapiscau site by MAIDENiso (smoothed over a 10 d period), without (red) and with (blue) the snow module. The black line shows the discharge from observations from the Caniapiscau basin between 1979–1997, scaled with the area of the basin. Solid lines indicate the average of the same DOY during this period, and shadows indicate the 2σ variability.


3.3 GPP calibration

GPP was calibrated twice at each station: first for the original version of MAIDENiso and second for the new version with snow. Both versions were able to predict GPP observations 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.

Figure 5GPP at the Uaf site (a) and the EOBS site (b). We show observations from flux towers (black) and MAIDENiso simulations without (red) and with (blue) snow processes. Independent optimizations are run for the two versions of MAIDENiso. A t test determined that the GPP simulated by the model with and without snow are statistically not different; for this reason, the red line is overlapped by the blue line.


A Welch t test determined that the GPP simulated by both versions of the model were statistically not different at the two sites of EOBS (t=0.04, p value=0.97) and Uaf (t=0.47, p value=0.63). This lack of influence of the snow module in the GPP was expected considering that while the snow module increases the availability of water in spring, the trees in our sites are not limited by water availability. This result also guarantees that the effects of the snow module on the δ18O outputs that we are testing below are not affected by GPP by means of Eq. (B4).

3.4 Effects of the snow module on the δ18O outputs

Both model versions (MAIDENiso with and without snow) reproduced the mean level of the δ18OTRC series (Fig. 6). This was expected because (1) the calibration process maximized the coincidence between observed and simulated δ18OTRC, and (2) the biochemical and kinetic fractionation parameters ϵ0 and ϵk could compensate the mean level of δ18OTRC in Eq. (B1) for differences in δ18OXW. Regarding the agreement between inter-annual variations, the observed and simulated δ18OTRC were not significantly correlated when snow was absent in the model (Fig. 6), but they were significantly correlated when MAIDENiso included snow (r=0.57 in Tungsten and r=0.52 in Caniapiscau, p<0.01; Fig. 6). In the case of Caniapiscau, none of the versions of the model were able to simulate the amplitude of the variability of the observed δ18OTRC, which makes it difficult to appreciate the improvement in the correlation between simulated and observed δ18OTRC when adding snow to the model in Fig. 6b. To facilitate the interpretation of Fig. 6a and b, we have standardized the observed and simulated δ18OTRC series (transformed to mean 0 and standard deviation 1) in Fig. 6c and d.

Figure 6The δ18OTRC observed (black) and simulated by the versions of MAIDENiso without (red) and with (blue) snow, using the calibrated parameters for (a, c) Tungsten and (b, d) Caniapiscau. Top panels (a, b) show the raw δ18OTRC series; bottom panels (c, d) show the standardized δ18OTRC series. Top-right inner panels show the scatter diagrams of the simulated (by both versions of the model, red for the version without snow and blue for the version with snow) and observed values of δ18OTRC, with dashed lines showing the linear regression models. Correlations are identical for the raw series (top panels a and b) and the standardized series (bottom panels c and d).


The distribution of the optimized parameters controlling δ18OTRC can help to assess how the model compensated for the absence of snow (Table 4 and Fig. 7). In our simulations, adding snow induced an increase in the dampening factor f0 at both sites, suggesting that the signal of the source water on δ18OTRC was stronger than without considering snow.

Figure 7Posterior probability density distributions of the parameters controlling δ18OTRC at the Tungsten and Caniapiscau sites, for the model without (red) and with (blue) snow.


To cast light on the effect of the snow module on δ18O, we compared the δ18O values from various parts of the model (precipitation, xylem water, discharge water, leaf water and TRC) for both versions (with and without snow) of MAIDENiso at the two study sites (Fig. 8). The δ18O of the precipitation (δ18OP, Fig. 8a) is the original source of the isotopic signals in the other parts of the model and it matches well to the IsoGSM data we used for its calibration, although it is consistently lower than the monthly data of the OIPC (Online Isotopes in Precipitation Calculator) at both sites. The snow directly impacted the δ18O of the source (xylem) water, δ18OXW (Fig. 8b). Without the snow module, δ18OXW followed closely the δ18OP signal shown in Fig. 8a, with a small delay due to the isotopic mixing in the soil. The δ18OXW signal was slightly enriched with respect to δ18OP due to isotopic fractionation associated with soil and canopy evaporation. In contrast, when snow was present in the model, the soil absorbed melted water from the snow pile in spring, which was enriched due to fractionation during sublimation in winter. The δ18OXW values were higher in the version with snow. The difference in δ18OXW values between the two versions reduced in time due to the infiltration and mixing of the enriched summer precipitation. Isotopic composition of discharge water (δ18Odis, Fig. 8c) follows a similar pattern to xylem water but with a delay of about 2 months and lower amplitudes in the seasonal variations.

Figure 8Simulated δ18O at different stages of the water cycle for the period 1979–2003 for the sites of Tungsten (left) and Caniapiscau (right). (a) The δ18O in precipitation (δ18OP) simulated by MAIDENiso using the NARR meteorology (black), monthly data from the Online Isotopes in Precipitation Calculator (OIPC, red) and mean from the IsoGSM dataset (green). (b) The δ18O in the xylem water (δ18OXW), (c) the δ18O in the discharge water (δ18Odis), (d) the δ18O in the leaf (δ18Oleaf) and (e) in the tree-ring cellulose (δ18OTRC) for the model without snow (red) and with snow (blue). Shadows indicate the 2σ variability for the same DOY within the 1979–2003 period. Vertical dashed lines in (d) and (e) indicate the start (budburst) and end of the growth season. Note that isotopic calculations are still made outside of this period, but no water is absorbed by the tree. Panels (d) and (e) include the t scores from a Welch t test, which show that the curves obtained for two versions of the model are different.


The daily δ18O in the leaf (δ18Oleaf) and the TRC calculated with Eqs. (B1) and (B2) are shown in Fig. 8d and e. The existence of the snow layer induced a higher δ18O at the leaf level, especially in spring, due to the higher δ18OXW level. The mean level of δ18OTRC did not change after adding snow, despite the enrichment of δ18Oleaf and δ18OXW, because it was compensated by the lower values of the biochemical fractionation ϵ0. A Welch t test was used to test the hypothesis that the curves obtained for two versions of the model were different, which was confirmed in all figures.

3.5 Relative influence of xylem water and leaf-level processes to the δ18OTRC signature

Finally, we investigated the relative contributions from the source water through δ18OXW and the leaf transpiration enrichment through δ18Oleaf on the δ18OTRC time series (Fig. 9, peak values in Table 5). At both sites, the leaf water δ18O isotopic enrichment had a stronger influence on δ18OTRC than the δ18O variability of xylem source water, as shown by the higher variance explained by the experiment that simulated δ18OTRC considering only the effect of δ18O leaf water enrichment indicated by a higher coefficient of determination (R2).

Figure 9Density distribution of the coefficients of determination (R2) between the reference simulations and the water source (xylem) experiments (solid line, δ18OV and hair set constant), and the leaf water enrichment experiments (dashed line, δ18OXW set constant) for the model without snow (red) and with snow (blue). Data are shown for (a) Tungsten and (b) Caniapiscau.


Table 5Mode of the probability density function for the coefficients of determination (R2) in Fig. 9 between the reference simulations and the water source (xylem) experiments and the leaf water enrichment experiments for the model without snow and with snow and for the two sites of Tungsten and Caniapiscau.

Download Print Version | Download XLSX

The addition of snow increased the R2 for both types of experiments but more importantly for the xylem source water experiment. As a consequence, the difference between the R2 of xylem and leaf experiments became smaller, although leaf transpiration still explained higher R2 at both sites. This was in agreement with the increase in f0 seen in Table 4 for both sites, pointing to an important influence of snow on the source water and on δ18OTRC.

4 Discussion

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 paragraphs, we discuss the impacts of the snow module addition on the different components of the model (i.e. the hydrological, photosynthetic and δ18O 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.

4.1 Improvements of the hydrological module

The calibration of the snow module at our sites yielded a value of ra almost twice as large in Tungsten than in Caniapiscau. Because of the lack of data for daily wind speed, we chose to implement ra as a constant parameter. The higher value of ra at Tungsten implies that snow sublimated at a slower rate, likely associated with the fact that the average wind speed during winter at this site was smaller than at Caniapiscau, in accordance to the interpretation of ra (Lawrence et al.2019; Blanken and Black2004).

The skills of MAIDENiso to reproduce the hydrological cycle improved with the implementation of the snow module. Because of the accumulation and melting of snow, the new version of MAIDENiso is now able to simulate the observed peak of river discharge in early spring, while the previous version without the snow module could not simulate any peak (Fig. 4). The magnitude of this peak cannot be compared directly with observations, because downscaling the river discharge by the size of the basin is not enough to make a direct comparison, as we also need to consider the following. (1) The water outflux simulated by MAIDENiso is the surface runoff (which is incorporated immediately to the streams) plus the subterranean drainage (which takes a longer time to reach the stream), which creates a time difference between outflux sources within the same spatial point. (2) The outflux from different points of the basin takes different times to reach the main stream of the basin. (3) The outflux over the whole area of the Caniapiscau basin is not necessarily identical. A routing model (Oki et al.1999; Southworth et al.2007; Miller et al.1994) could be used to calculate the delay and flow to the main stream from across the whole basin for both types of water outflows. The incorporation of a routing model in MAIDENiso would allow us to produce an estimate of streamflow for a basin, allowing for direct comparison with river discharge observations.

4.2 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. The parameters controlling GPP that we obtained showed very similar posterior distributions and values in the plausible block (Table A1 and Figs. A1 and A2). These similarities indicate that, at our study sites, photosynthesis was not very sensitive to additional water from snowmelt, suggesting that radial tree growth was not limited by water availability. This 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 Running2006; D'Orangeville et al.2018). However, different results could be found in sites where trees are more dependent on water derived 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 δ18OTRC variations, because it allows us to discard GPP as a possible cause for the differences observed between the two model versions.

4.3 Effects of snow on xylem water, leaf and tree-ring cellulose δ18O

Following the approach proposed by Lavergne et al. (2017), we produced yearly δ18OTRC time series by weighting the daily values with the GPP, assuming that C allocation to the stem is proportional to GPP. MAIDENiso has a module for the allocation of available C to the different parts of the tree, which provides an alternative and more realistic way of calculating yearly δ18OTRC. However, the calibration of this module ideally requires observations of the same units as the product of the allocation module, i.e. C mass per unit of stand basal area allocated to the stem. Although tree-ring width (TRW) data were available for both sites, their use was complicated as TRW observations represent just a portion of the total C allocation of the entire tree and do not offer an intra-annual C allocation resolution to constrain the simulations. Therefore, the use of GPP to weight the daily δ18OTRC was the best option for this particular study. The model without snow produced depleted values of δ18OXW, because it lacked the enrichment effect of evaporative fractionation in snow. This was also reflected in the δ18Odis signal, with depleted values for the model without snow. This contradicts studies that show that stream water does not usually show signs of evaporative fractionation (Evaristo et al.2015). However, the discharge in MAIDENiso is simply the addition of the runoff and the drainage from the upper soil layers, and this lacks the complexity of interactions with deeper groundwater that the real soil has. Therefore, while orientational, δ18Odis is only useful as a possible input to a more complex groundwater model.

The addition of the snow module corrects for an important overcompensation effect that stemmed from the unrealistic representation of hydrology in the previous version. The model without snow produced depleted values of δ18OXW, which the model had to compensate through the δ18O parameters (higher values of ϵ0 and ϵk and lower contributions of xylem water through lower f0). The calibration of the δ18O processes for the two versions of MAIDENiso and the two sites yielded significant differences in the optimized parameters, both in the distribution of the optimal blocks and the values in the plausible blocks (Fig. 7 and Table 4). The dampening factor f0, which controls the direct contribution of the source water to the δ18OTRC signal, was significantly higher (especially in Tungsten) after adding snow. The calibration of the model without snow converged to a value of f0≈0.32, with the posterior distributions pushing toward the lower prior limit of 0.3 (Table 4 and Fig. 7), which suggests that the calibration procedure would have converged towards a smaller value if it had been allowed. Adding snow increased the dampening factor to f0=0.48 in Tungsten and f0=0.43 in Caniapiscau, in agreement with the range of f0=0.4–0.5 reported in previous studies (Roden and Ehleringer2000; Saurer et al.1997; Sternberg et al.1986; Yakir1992). Lavergne et al. (2017) obtained a dampening factor of f0=0.41 in Quebec using the original model as the parameters from Eq. (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 δ18OTRC.

The calibrated value for the biochemical fractionation ϵ0 was different at the two sites, ranging with snow to without snow from 27.41 ‰ to 27.85 ‰ in Tungsten and from 24.15 ‰ to 24.48 ‰ in Caniapiscau (Table 4). The ϵ0 values were slightly higher at both sites when the model lacked snow, which suggests that the calibration compensates for consistently lower values of δ18OXW and/or δ18Oleaf in Eq. (B1) to adjust the mean δ18OTRC to the observations. The kinetic fractionation ϵk obtained also differed strongly between sites, with snow to without snow from 11.8 ‰ to 13 ‰ at Tungsten and from 22.8 ‰ to 23.4 ‰ at Caniapiscau (Table 4). The ϵk was set to 26.5 ‰ by Farquhar et al. (1989), but it can vary over a larger range (Buhay et al.1996). Lavergne et al. (2017) obtained a value of ϵk=17.20 ‰ for Quebec, with a similar posterior distribution that the one obtained here.

Our results also showed that the leaf 18O enrichment due to transpiration has a stronger influence on δ18OTRC than the isotopic composition of the source (xylem) water, both in Tungsten and Caniapiscau, suggesting that it is the main driver of δ18OTRC variations. These results are in agreement with Lavergne et al. (2017) findings and reflect the strong effect of vapour pressure deficit on δ18Oleaf in Quebec. Nevertheless, the δ18OTRC signature also had a strong imprint of the source water signal as recently reported for the Tungsten δ18OTRC record that shared the same large-scale atmospheric patterns than spring–summer δ18OP (Field et al.2021).

The addition of the snow module to MAIDENiso therefore frees the calibration process from having to overcompensate for the artificially depleted δ18OXW values during the growing season (Fig. 8b). As our results have shown, this significantly increased the correlation between the observed and simulated δ18OTRC compared to the version without snow (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 δ18OTRC indirectly through the source water.

Overall, the improvements found in the δ18OTRC simulations at both sites indicate that snow plays a critical role in δ18O of the source water and thus on the final signature of δ18OTRC. Even if the addition of snow would not had resulted in a significant improvement of the correlation between the simulated and observed δ18OTRC, 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. Process-based 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 calibrated under favourable conditions and then used for different datasets (Guiot et al.2014). 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.

4.4 Implications for future studies

Based on our results and comparison with other studies, we can conclude that the snow module predicted more realistic 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 δ18Oleaf enrichment signal on δ18OTRC 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 English2015). 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 studies that did not take into account the additional effects of snowpack dynamics on δ18OTRC, e.g. Lavergne et al. (2017). Now, MAIDENiso can simulate more reliable interactions between the coupled water and carbon cycles and tree physiological mechanisms in cold environments. Our findings will contribute to reduce uncertainties in the predictions of the response of forest productivity to hydrological changes, leading to better forward predictions that can eventually be used to reconstruct seasonal and long-term hydroclimatic variations.

An inverse modelling approach has previously been developed and tested using MAIDENiso to reconstruct paleoclimate from tree-ring data in the Fontainebleau Forest, France (Boucher et al.2014). However, this exercise was restricted to the reconstruction of meteorological variables during summer and to regions where the tree-ring proxies were not significantly affected by winter meteorology. The inclusion of snow in the model opens new possibilities for reconstructing hydroclimate in cold regions, considering that the new version of MAIDENiso produces simulation of δ18OTRC that account for snow-related processes.

Suitable regions for the application of MAIDENiso in future studies include high-mountain regions, now that the model has a working snow module. Regions with snow-dominated winters and dry summers, such as the southwestern USA or some Mediterranean sites, can also be of interest for future studies with MAIDENiso, as trees in these sites can be more dependent on water derived from snowmelt than the sites used in the present study. The application of MAIDENiso to any site is possible provided that there is sufficient meteorological data to drive the model and local information to calibrate the model parameters (GPP, snow, TRW and TRC stable isotopes). We expect MAIDENiso to be applied more broadly in high-latitude and high-altitude environments in the near future.

5 Conclusions

In this paper we presented the new snow module incorporated into MAIDENiso, which consists of new hydrological calculations of snow dynamics and a thermal module. Our results show how this snow module improves the simulation of outputs associated with the hydrological cycle at cold and high-latitude sites without affecting simulations from the carbon cycle component. These findings were expected as 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 Linder2000). The simulations of the new version of MAIDENiso reproduce the observed δ18OTRC better than the original snowless version of MAIDENiso. Based on the development presented here, the potential for the application of MAIDENiso is notably increased.

Appendix A: Photosynthesis model

GPP (gCm-2d-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 following De Pury and Farquhar (1997) as explained in Misson (2004). Daily Vcmax (Vcmaxi) is modelled as

(A1) Vcmax i = Vmax 1 + exp ( Vb ( Tday i - Vip ) ) .

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.

The stomatal conductance for carbon (µmolm-2s-1) is calculated using the Leuning et al. (1995) model, modified by Gea-Izquierdo et al. (2015) to incorporate soil water stress:

(A2) g sc = g 0 + g 1 A n ( C a - Γ ) ( 1 + VPD / VPD 0 ) θ g P atm ,

where g0=0µmolm-2s-1 and g1=10µmolm-2s-1 are fitted parameters representing the residual conductance as the net assimilation rate (An) approaches zero and the slope of the function, respectively. Patm is the atmospheric pressure (Pa). Ca is the atmospheric CO2 pressure (Pa). Γ is the CO2 compensation point in the absence of dark respiration (Pa), which is calculated following Bernacchi et al. (2001). VPD is the vapour pressure deficit (kPa), and VPD0 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:

(A3) θ g = 1 1 + exp ( soilb ( SWC i - soilip ) ) .

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 account by the model. This is done by replacing Tday in Eq. (A1) by the temperature transformation S, defined as

(A4) d S i d i = Tday i - S i τ ,

where τ is a parameter representing the number of days needed by the tree to adapt the photosynthesis to changing temperatures.

There are a total of six undetermined parameters that control GPP production in MAIDENiso in Eqs. (A1), (A3) and (A4). These parameters were calibrated at the two eddy covariance flux stations described in Sect. 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.

Table A1Calibration parameters for the GPP module. Parameters, units, prior range and posterior range are shown (with parameter value in the plausible block) for both MAIDENiso versions and both flux towers.

n/a – not applicable

Download Print Version | Download XLSX

Figure A1Posterior probability density distributions of the parameters controlling GPP at the Uaf site for the model without (red) and with (blue) snow.


Figure A2Posterior probability density distributions of the parameters controlling GPP at the EOBS site for the model without (red) and with (blue) snow.


Appendix B: Isotopic model

MAIDENiso keeps track of the stable isotopic composition of oxygen (δ18O) in all the water/ice pools and fluxes of the hydrological model (Fig. 2). The isotopic module calculates fractionation from evaporation (from soil and canopy water) and transpiration at leaf level to produce an isotopic oxygen signature in TRC (δ18OTRC). This is based on the Danis et al. (2012) formulation of the Craig–Gordon model (Craig and Gordon1965):

(B1) δ 18 O TRC = ( 1 - f 0 ) δ 18 O leaf + f 0 δ 18 O XW + ϵ 0 ,

with δ18O at leaf level being


Here, f0 (unitless) is the dampening factor reflecting the exchange of the oxygen atoms between sucrose and xylem water during the synthesis of cellulose in the xylem cells of the tree rings, typically within a range of 0.4–0.5 (Roden and Ehleringer2000; Saurer et al.1997; Sternberg et al.1986; Yakir1992). ϵ0 is the biochemical fractionation due to oxygen exchange between water and the carbonyl groups (C=O) in the organic molecules, undetermined but expected in a range of 24 ‰–30 ‰ (DeNiro and Epstein1979; Farquhar et al.1998). ϵ is the equilibrium fractionation due to the change of phase of water from liquid to vapour at leaf temperature (fixed at 21.4 C, which is the temperature threshold for maximum carbon assimilation), with a value of 9.65 ‰ (Helliker and Richter2008). ϵk is the kinetic fractionation due to the diffusion of vapour 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 as it can vary over larger ranges (Buhay et al.1996). hair is the relative humidity, which is estimated in MAIDENiso from the daily air temperature and the dew point temperature (Running et al.1987). δ18OV and δ18OXW are the δ18O of vapour and xylem (source) water, respectively. δ18OV is calculated from the δ18O of precipitation (δ18OP) and the fractionation due to the phase change from liquid water to vapour at mean air temperature, ϵTair (Horita and Wesolowski1994):

(B3) δ 18 O V = δ 18 O P - ϵ T air .

The δ18OTRC time series produced through Eq. (B1) are daily, while the δ18OTRC 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 (GPPd):

(B4) δ 18 O TRC , y = d δ 18 O TRC , d GPP d d GPP d .
Code and data availability

The MAIDENiso code is available in the Zenodo repository. Note that there are two versions of the code, corresponding to the model with snow (, Hermoso de Mendoza2021a) and the model without snow (, Hermoso de Mendoza2021b). In addition, a university website has been created (, last access: 1 December 2021) for the MAIDEN model, where a technical description of the model and access to different model versions will be available. The meteorological input files and parameter files needed to run MAIDENiso and the observational data are available in the Zenodo repository (, Hermoso de Mendoza2021c).

Author contributions

IHdM implemented the new snow module in MAIDENiso, performed the simulations and analyses, and wrote the first draft of the manuscript. FG and AL helped in setting up the new module in the original version of the model. LAH, RF and EB provided the δ18OTRC chronology. EB, LAH, FG and AL contributed to the design of the study (analyses to perform and structure of the paper) and interpretation of the results. All authors contributed to improve the article and guided the simulations and analyses.

Competing interests

The contact author has declared that neither they nor their co-authors have any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


This research has been supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) PERSISTENCE project (grant number RDC 485475 - 15 to Etienne Boucher) and the US National Science Foundation (NSF) (grant numbers PLR-1504134, PLR-1603473, AGS- 1502150 and OISE-1743738 to Laia Andreu-Hayles). The PERSISTENCE project is a collaborative research grant that involves the participation of Hydro-Québec, Manitoba Hydro and the Ouranos consortium. We thank Luc Perreault and Dominique Tapsoba from Hydro-Québec for providing the SWE data at the Caniapiscau basin used here. Aliénor Lavergne was supported by a Marie Sklodowska-Curie Individual Fellowship under the European Union's Horizon 2020 Research and Innovation Programme (grant agreement number: 838739 ECAW-ISO). Fabio Gennaretti was supported by the Ministère des Forêts, de la Faune et des Parcs (MFFP; contract number 142332177-D) and NSERC (Alliance grant number ALLRP 557148-20). We are thankful to Wei Huang and Jean-Francois Hélie respectively, from the Stable Isotope Laboratory of Lamont-Doherty Earth and GEOTOP (UQAM), for their support with isotopic measurements.

Financial support

This research has been supported by the Natural Sciences and Engineering Research Council of Canada (NSERC) PERSISTENCE project (grant number RDC 485475 - 15 to Etienne Boucher), and the US National Science Foundation (NSF) (grant numbers PLR-1504134, PLR-1603473, AGS-1502150 and OISE-1743738 to Laia Andreu-Hayles).

Review statement

This paper was edited by Tomomichi Kato and reviewed by Vladimir Shishov and one anonymous referee.


Barnhart, T. B., Molotch, N. P., Livneh, B., Harpold, A. A., Knowles, J. F., and Schneider, D.: Snowmelt rate dictates streamflow, Geophys. Res. Lett., 43, 8006–8016,, 2016. a

Bergeron, O., Margolis, H. A., Black, T. A., Coursolle, C., Dunn, A. L., Barr, A. G., and Wofsy, S. C.: Comparison of carbon dioxide fluxes over three boreal black spruce forests in Canada, Glob. Change Biol., 13, 89–107, doi10.1111/j.1365-2486.2006.01281.x, 2007 a

Berghuijs, W., Woods, R., and Hrachowitz, M.: A precipitation shift from snow towards rain leads to a decrease in streamflow, Nat. Clim. Change, 4, 583–586,, 2014. a

Beria, H., Larsen, J. R., Ceperley, N. C., Michelon, A., Vennemann, T., and Schaefli, B.: Understanding snow hydrological processes through the lens of stable water isotopes, Wiley Interdisciplinary Reviews: Water, 5, e1311,, 2018. a, b

Bernacchi, C., Singsaas, E., Pimentel, C., Portis Jr, A., and Long, S. P.: Improved temperature response functions for models of Rubisco-limited photosynthesis, Plant Cell Environ., 24, 253–259,, 2001. a

Blanken, P. and Black, T.: The canopy conductance of a boreal aspen forest, Prince Albert National Park, Canada, Hydrol. Process., 18, 1561–1578,, 2004. a, b

Boisvenue, C. and Running, S. W.: Impacts of climate change on natural forest productivity–evidence since the middle of the 20th century, Glob. Change Biol., 12, 862–882,, 2006. a

Boivin, S. and Bégin, Y.: Development of a black spruce (Picea mariana) shoreline stand in relation to snow level variations at Lake Bienville in northern Quebec, Can. J. Forest Res., 27, 295–303, 1997. a

Boucher, É., Guiot, J., Hatté, C., Daux, V., Danis, P.-A., and Dussouillez, P.: An inverse modeling approach for tree-ring-based climate reconstructions under changing atmospheric CO2 concentrations, Biogeosciences, 11, 3245–3258,, 2014. a, b

Boucher, E., Nicault, A., Arseneault, D., Bégin, Y., and Karami, M.: Decadal variations in Eastern Canada's taiga wood biomass production forced by ocean-atmosphere interactions, Sci. Rep., 7, 2457,, 2017. a

Buhay, W., Edwards, T., and Aravena, R.: Evaluating kinetic fractionation factors used for ecologic and paleoclimatic reconstructions from oxygen and hydrogen isotope ratios in plant water and cellulose, Geochim. Cosmochim. Ac., 60, 2209–2218,, 1996. a, b

Cernusak, L. A. and English, N. B.: Beyond tree-ring widths: stable isotopes sharpen the focus on climate responses of temperate forest trees, Tree Physiol., 35, 1–3,, 2015. a

Coulthard, B. L., Anchukaitis, K. J., Pederson, G. T., Cook, E., Littell, J., and Smith, D. J.: Snowpack signals in North American tree rings, Environ. Res. Lett., 16, 034037,, 2021. a, b

Craig, H. and Gordon, L. I.: Deuterium and oxygen 18 variations in the ocean and the marine atmosphere, Consiglio nazionale delle richerche, Laboratorio de geologia nucleare Pisa, 1965. a

Danis, P.-A., Hatté, C., Misson, L., and Guiot, J.: MAIDENiso: a multiproxy biophysical model of tree-ring width and oxygen and carbon isotopes, Can. J. Forest Res., 42, 1697–1713,, 2012. a, b, c, d

De Pury, D. and Farquhar, G.: Simple scaling of photosynthesis from leaves to canopies without the errors of big-leaf models, Plant Cell Environ., 20, 537–557,, 1997. a

DeNiro, M. J. and Epstein, S.: Relationship between the oxygen isotope ratios of terrestrial plant cellulose, carbon dioxide, and water, Science, 204, 51–53,, 1979. a

D'Orangeville, L., Houle, D., Duchesne, L., Phillips, R. P., Bergeron, Y., and Kneeshaw, D.: Beneficial effects of climate warming on boreal tree growth may be transitory, Nat. Commun., 9, 1–10,, 2018. a

Drew, D. M., Downes, G. M., and Battaglia, M.: CAMBIUM, a process-based model of daily xylem development in Eucalyptus, J. Theor. Biol., 264, 395–406,, 2010. a

Du, J., He, Z., Yang, J., Chen, L., and Zhu, X.: Detecting the effects of climate change on canopy phenology in coniferous forests in semi-arid mountain regions of China, Int. J. Remote Sens., 35, 6490–6507,, 2014. a

Evaristo, J., Jasechko, S., and McDonnell, J. J.: Global separation of plant transpiration from groundwater and streamflow, Nature, 525, 91–94,, 2015. a

Farquhar, G., Hubick, K., Condon, A., and Richards, R.: Carbon isotope fractionation and plant water-use efficiency, in: Stable isotopes in ecological research, Springer, 21–40,, 1989. a, b

Farquhar, G., Barbour, M., and Henry, B.: Interpretation of oxygen isotope composition of leaf material, in: Stable isotopes: Integration of Biological, Ecological and Geochemical Processes, Garland Science, ISBN: 9781003076865, 27–62, 1998. a

Farquhar, G. D., von Caemmerer, S., and Berry, J. A.: A biochemical model of photosynthetic CO2 assimilation in leaves of C3 species, Planta, 149, 78–90,, 1980. a

Field, R., Andreu-Hayles, L., D'Arrigo, R., Oelkers, R., Luckman, B., Morimoto, D., Boucher, E., Gennaretti, F., Hermoso, I., Lavergne, A., and Lévesque, M.: Tree-ring Cellulose δ18O Records Similar Large-scale Climate Influences as Precipitation δ18O in the Northwest Territories of Canada, Clim. Dynam., 58, 759–776,, 2021. a, b

Fritts, H. C., Vaganov, E. A., Sviderskaya, I. V., and Shashkin, A. V.: Climatic variation and tree-ring structure in conifers: empirical and mechanistic models of tree-ring width, number of cells, cell size, cell-wall thickness and wood density, Clim. Res., 1, 97–116,, 1991. a

Gaucherel, C., Campillo, F., Misson, L., Guiot, J., and Boreux, J.-J.: Parameterization of a process-based tree-growth model: comparison of optimization, MCMC and particle filtering algorithms, Environ. Modell. Softw., 23, 1280–1288,, 2008. a

Gea-Izquierdo, G., Guibal, F., Joffre, R., Ourcival, J. M., Simioni, G., and Guiot, J.: Modelling the climatic drivers determining photosynthesis and carbon allocation in evergreen Mediterranean forests using multiproxy long time series, Biogeosciences, 12, 3695–3712,, 2015. a, b, c

Gennaretti, F., Gea-Izquierdo, G., Boucher, E., Berninger, F., Arseneault, D., and Guiot, J.: Ecophysiological modeling of photosynthesis and carbon allocation to the tree stem in the boreal forest, Biogeosciences, 14, 4851–4866,, 2017a. a, b, c, d, e, f, g, h

Gennaretti, F., Huard, D., Naulier, M., Savard, M., Bégin, C., Arseneault, D., and Guiot, J.: Bayesian multiproxy temperature reconstruction with black spruce ring widths and stable isotopes from the northern Quebec taiga, Clim. Dynam., 0, 4107–4119,, 2017b. a

Girardin, M. P., Raulier, F., Bernier, P. Y., and Tardif, J. C.: Response of tree growth to a changing climate in boreal central Canada: a comparison of empirical, process-based, and hybrid modelling approaches, Ecol. Model., 213, 209–228,, 2008. a

Guiot, J., Boucher, E., and Gea-Izquierdo, G.: Process models and model-data fusion in dendroecology, Frontiers in Ecology and Evolution, 2, 52,, 2014. a, b

Gupta, H. V., Kling, H., Yilmaz, K. K., and Martinez, G. F.: Decomposition of the mean squared error and NSE performance criteria: Implications for improving hydrological modelling, J. Hydrol., 377, 80–91,, 2009. a

Helliker, B. R. and Richter, S. L.: Subtropical to boreal convergence of tree-leaf temperatures, Nature, 454, 511–514,, 2008. a

Hermoso de Mendoza, I., Boucher, E., Gennaretti, F., and Lavergne, A.: MAIDENiso (Version v2021), Zenodo [code],, 2021a. a

Hermoso de Mendoza, I., Boucher, E., Gennaretti, F., and Lavergne, A.: MAIDENiso-snowless (Version v2021), Zenodo [code],, 2021b. a

Hermoso de Mendoza, I., Boucher, E., Andreu-Hayles, L., Gennaretti, F., and Lavergne, A.: Meteorology, MAIDENiso parameter files and observational data for Tungsten, Caniapiscau, and associated GPP stations, Zenodo [data set],, 2021c. a

Hölttä, T., Mäkinen, H., Nöjd, P., Mäkelä, A., and Nikinmaa, E.: A physiological model of softwood cambial growth, Tree Physiol., 30, 1235–1252,, 2010. a

Horita, J. and Wesolowski, D. J.: Liquid-vapor fractionation of oxygen and hydrogen isotopes of water from the freezing to the critical temperature, Geochim. Cosmochim. Ac., 58, 3425–3437,, 1994. a

Huang, R., Zhu, H., Liang, E., Grießinger, J., Wernicke, J., Yu, W., Hochreuther, P., Risi, C., Zeng, Y., Fremme, A., Sodemann, H., and Bräuning, A.: Temperature signals in tree-ring oxygen isotope series from the northern slope of the Himalaya, Earth Planet. Sc. Lett., 506, 455–465,, 2019. a

Hutchinson, M. F., McKenney, D. W., Lawrence, K., Pedlar, J. H., Hopkinson, R. F., Milewska, E., and Papadopol, P.: Development and testing of Canada-wide interpolated spatial models of daily minimum–maximum temperature and precipitation for 1961–2003, J. Appl. Meteorol. Clim., 48, 725–741,, 2009. a

Jarvis, P. and Linder, S.: Constraints to growth of boreal forests, Nature, 405, 904–905, 2000. a

Keeling, C. D., Bacastow, R. B., Bainbridge, A. E., Ekdahl Jr, C. A., Guenther, P. R., Waterman, L. S., and Chin, J. F.: Atmospheric carbon dioxide variations at Mauna Loa observatory, Hawaii, Tellus, 28, 538–551,, 1976. a

Knauer, J., Werner, C., and Zaehle, S.: Evaluating stomatal models and their atmospheric drought response in a land surface scheme: A multibiome analysis, J. Geophys. Res.-Biogeo., 120, 1894–1911,, 2015. a

Knoben, W. J. M., Freer, J. E., and Woods, R. A.: Technical note: Inherent benchmark or not? Comparing Nash–Sutcliffe and Kling–Gupta efficiency scores, Hydrol. Earth Syst. Sci., 23, 4323–4331,, 2019. a

Kurita, N., Yoshida, N., Inoue, G., and Chayanova, E. A.: Modern isotope climatology of Russia: A first assessment, J. Geophys. Res.-Atmos., 109, D03104,, 2004. a

Lavergne, A., Gennaretti, F., Risi, C., Daux, V., Boucher, E., Savard, M. M., Naulier, M., Villalba, R., Bégin, C., and Guiot, J.: Modelling tree ring cellulose δ18O variations in two temperature-sensitive tree species from North and South America, Clim. Past, 13, 1515–1526,, 2017. a, b, c, d, e, f, g

Lawrence, D. M., Fisher, R. A., Koven, C. D., Oleson, K. W., Swenson, S. C., Bonan, G., Collier, N., Ghimire, B., van Kampenhout, L., Kennedy, D., Kluzek, E., Lawrence, P. J., Li, F., Li, H., Lombardozzi, D., Riley, W. J., Sacks, W. J., Shi, M., Vertenstein, M., Wieder, W. R., Xu, C.,, Ali, A. A., Badger, A. M., Bisht, G., van den Broeke, M., Brunke, M. A., Burns, S. P., Buzan, J., Clark, M., Craig, A., Dahlin, K., Drewniak, B., Fisher, J. B., Flanner, M., Fox, A. M., Gentine, P., Hoffman, F., Keppel-Aleks, G., Knox, R., Kumar, S., Lenaerts, J., Leung, L. R., Lipscomb, W. H., Lu, Y., Pandey, A., Pelletier, J. D., Perket, J., Randerson, J. T., Ricciuto, D. M., Sanderson, B. M., Slater, A., Subin, Z. M., Tang, J., Thomas, R. Q., Val Martin, M., and Zeng, X.: The Community Land Model version 5: Description of new features, benchmarking, and impact of forcing uncertainty, J. Adv. Model. Earth Sy., 11, 4245–4287,, 2019. a, b, c

Leuning, R., Kelliher, F. M., De Pury, D., and Schulze, E.-D.: Leaf nitrogen, photosynthesis, conductance and transpiration: scaling from leaves to canopies, Plant Cell Environ., 18, 1183–1200,, 1995. a

Li, D., Wrzesien, M. L., Durand, M., Adam, J., and Lettenmaier, D. P.: How much runoff originates as snow in the western United States, and how will that change in the future?, Geophys. Res. Lett., 44, 6163–6172,, 2017. a

Li, G., Harrison, S. P., Prentice, I. C., and Falster, D.: Simulation of tree-ring widths with a model for primary production, carbon allocation, and growth, Biogeosciences, 11, 6711–6724,, 2014. a

Loescher, H., Hanson, C., and Ocheltree, T.: The psychrometric constant is not constant: A novel approach to enhance the accuracy and precision of latent energy fluxes through automated water vapor calibrations, J. Hydrometeorol., 10, 1271–1284,, 2009. a

McCabe, G. J. and Wolock, D. M.: Recent declines in western US snowpack in the context of twentieth-century climate variability, Earth Interact., 13, 1–15,, 2009. a

Meredith, M., Sommerkorn, M., Cassotta, S., Derksen, C., Ekaykin, A., Hollowed, A., Kofinas, G., Mackintosh, A., Melbourne-Thomas, J., Muelbert, M., Ottersen, G., Pritchard, H., and Schuur, E. A. G.: Polar Regions, chap. 3, IPCC Special Report on the Ocean and Cryosphere in a Changing Climate, ISBN-10: 1009157973, 2019. a

Mesinger, F., DiMego, G., Kalnay, E., Mitchell, K., Shafran, P. C., Ebisuzaki, W., Jović, D., Woollen, J., Rogers, E., Berbery, E. H., B. Ek, M., Fan, Y., Grumbine, R., Higgins, W., Li, H., Lin, Y., Manikin, G., Parrish, D., and Shi, W.: North American regional reanalysis, B. Am. Meteorol. Soc., 87, 343–360,, 2006. a, b, c

Miller, J. R., Russell, G. L., and Caliri, G.: Continental-scale river flow in climate models, J. Climate, 7, 914–928,<0914:CSRFIC>2.0.CO;2, 1994. a

Misson, L.: MAIDEN: a model for analyzing ecosystem processes in dendroecology, Can. J. Forest Res., 34, 874–887,, 2004. a, b, c, d, e, f

Nash, J. E. and Sutcliffe, J. V.: River flow forecasting through conceptual models part I–A discussion of principles, J. Hydrol., 10, 282–290, 1970. a

Naulier, M., Savard, M. M., Bégin, C., Marion, J., Arseneault, D., and Bégin, Y.: Carbon and oxygen isotopes of lakeshore black spruce trees in northeastern Canada as proxies for climatic reconstruction, Chem. Geol., 374, 37–43, 2014. a

Naulier, M., Savard, M. M., Bégin, C., Gennaretti, F., Arseneault, D., Marion, J., Nicault, A., and Bégin, Y.: A millennial summer temperature reconstruction for northeastern Canada using oxygen isotopes in subfossil trees, Clim. Past, 11, 1153–1164,, 2015a. a

Naulier, M., Savard, M. M., Bégin, C., Marion, J., Nicault, A., and Bégin, Y.: Temporal instability of isotopes–climate statistical relationships–A study of black spruce trees in northeastern Canada, Dendrochronologia, 34, 33–42, 2015b. a

Nicault, A., Boucher, E., Bégin, C., Guiot, J., Marion, J., Perreault, L., Roy, R., Savard, M. M., and Bégin, Y.: Hydrological reconstruction from tree-ring multi-proxies over the last two centuries at the Caniapiscau Reservoir, northern Québec, Canada, J. Hydrol., 513, 435–445,, 2014. a

Nicault, A., Boucher, E., Tapsoba, D., Arseneault, D., Berninger, F., Bégin, C., DesGranges, J.-L., Guiot, J., Marion, J., Wicha, S., and Bégin, Y.: Spatial analysis of black spruce (Picea mariana (Mill.) BSP) radial growth response to climate in northern Québec–Labrador Peninsula, Canada, Can. J. Forest Res., 45, 343–352,, 2015. a

Nossent, J. and Bauwens, W.: Application of a normalized Nash-Sutcliffe efficiency to improve the accuracy of the Sobol'sensitivity analysis of a hydrological model, Math. Comput. Simulat., 55, 271–280, 2011. a

Oki, T., Nishimura, T., and Dirmeyer, P.: Assessment of annual runoff from land surface models using Total Runoff Integrating Pathways (TRIP), J. Meteorol. Soc. Jpn. Ser. II, 77, 235–255,, 1999. a

Ols, C., Trouet, V., Girardin, M. P., Hofgaard, A., Bergeron, Y., and Drobyshev, I.: Post-1980 shifts in the sensitivity of boreal tree growth to North Atlantic Ocean dynamics and seasonal climate, Global Planet. Change, 165, 1–12, 2018. a

Perkins, D. L. and Swetnam, T. W.: A dendroecological assessment of whitebark pine in the Sawtooth–Salmon River region, Idaho, Can. J. Forest Res., 26, 2123–2133,, 1996. a

Peters, W., Jacobson, A. R., Sweeney, C., Andrews, A. E., Conway, T. J., Masarie, K., Miller, J. B., Bruhwiler, L. M., Pétron, G., Hirsch, A. I., Worthy, D. E. J., van der Werf, G. R., Randerson, J. T., Wennberg, P. O., Krol, M. C., and Tans, P. P.: An atmospheric perspective on North American carbon dioxide exchange: CarbonTracker, P. Natl. Acad. Sci. USA, 104, 18925–18930,, 2007. a

Peterson, D. W. and Peterson, D. L.: Effects of climate on radial growth of subalpine conifers in the North Cascade Mountains, Can. J. Forest Res., 24, 1921–1932,, 1994. a

Porter, T. J., Pisaric, M. F., Field, R. D., Kokelj, S. V., Edwards, T. W., deMontigny, P., Healy, R., and LeGrande, A. N.: Spring-summer temperatures since AD 1780 reconstructed from stable oxygen isotope ratios in white spruce tree-rings from the Mackenzie Delta, northwestern Canada, Clim. Dynam., 42, 771–785,, 2014. a

Rathgeber, C., Nicault, A., Kaplan, J. O., and Guiot, J.: Using a biogeochemistry model in simulating forests productivity responses to climatic change and [CO2] increase: example of Pinus halepensis in Provence (south-east France), Ecol. Model., 166, 239–255,, 2003. a, b

Roden, J. S. and Ehleringer, J. R.: Hydrogen and oxygen isotope ratios of tree ring cellulose for field-grown riparian trees, Oecologia, 123, 481–489,, 2000. a, b

Running, S. W., Nemani, R. R., and Hungerford, R. D.: Extrapolation of synoptic meteorological data in mountainous terrain and its use for simulating forest evapotranspiration and photosynthesis, Can. J. Forest Res., 17, 472–483,, 1987. a

Rutter, N., Essery, R., Pomeroy, J., Altimir, N., Andreadis, K., Baker, I., Barr, A., Bartlett, P., Boone, A., Deng, H., Dou- ville, H., Dutra, E., Elder, K., Ellis, C., Feng, X., Gelfan, A., Goodbody, A., Gusev, Y., Gustafsson, D., Hellström, R., Hirabayashi, Y., Hirota, T., Jonas, T., Koren, V., Kuragina, A., Lettenmaier, D., Li, W. P., Luce, C., Martin, E., Nasonova, O., Pumpanen, J., Pyles, R. D., Samuelsson, P., Sandells, M., Schädler, G., Shmakin, A., Smirnova, T. G., Stähli, M., Stöckli, R., Strasser, U., Su, H., Suzuki, K., Takata, K., Tanaka, K., Thompson, E., Vesala, T., Viterbo, P., Wiltshire, A., Xia, K., Xue, Y., and Yamazaki, T.: Evaluation of forest snow processes models (SnowMIP2), J. Geophys. Res.-Atmos., 114, D06111,, 2009. a

Saurer, M., Aellen, K., and Siegwolf, R.: Correlating δ13C and δ18O in cellulose of trees, Plant Cell Environ., 20, 1543–1550,, 1997. a, b

Shishov, V. V., Tychkov, I. I., Popkova, M. I., Ilyin, V. A., Bryukhanova, M. V., and Kirdyanov, A. V.: VS-oscilloscope: a new tool to parameterize tree radial growth based on climate conditions, Dendrochronologia, 39, 42–50,, 2016. a, b

Southworth, F., Peterson, B., and Lambert, B.: Development of a regional routing model for strategic waterway analysis, Transp. Res. Record, 1993, 109–116,, 2007. a

St. George, S., Meko, D. M., Girardin, M.-P., MacDonald, G. M., Nielsen, E., Pederson, G. T., Sauchyn, D. J., Tardif, J. C., and Watson, E.: The tree-ring record of drought on the Canadian Prairies, J. Climate, 22, 689–710,, 2009. a

Sternberg, L. D. S., Deniro, M. J., and Savidge, R. A.: Oxygen isotope exchange between metabolites and water during biochemical reactions leading to cellulose synthesis, Plant Physiol., 82, 423–427,, 1986. a, b

Stigter, E. E., Litt, M., Steiner, J. F., Bonekamp, P. N., Shea, J. M., Bierkens, M. F., and Immerzeel, W. W.: The importance of snow sublimation on a Himalayan glacier, Front. Earth Sci.-Cryosphere, 6, 108,, 2018. a

Szejner, P., Wright, W. E., Belmecheri, S., Meko, D., Leavitt, S. W., Ehleringer, J. R., and Monson, R. K.: Disentangling seasonal and interannual legacies from inferred patterns of forest water and carbon cycling using tree-ring stable isotopes, Glob. Change Biol., 24, 5332–5347,, 2018.  a

Ueyama, M., Iwata, H., and Harazono, Y.: AmeriFlux US-Uaf University of Alaska, Fairbanks, Ver. 8.5,, 2021. a

Vaganov, E. A., Hughes, M. K., and Shashkin, A. V.: Growth dynamics of conifer tree rings: images of past and future environments, vol. 183, Springer Science & Business Media,, 2006. a

van Kampenhout, L., Lenaerts, J. T., Lipscomb, W. H., Sacks, W. J., Lawrence, D. M., Slater, A. G., and van den Broeke, M. R.: Improving the representation of polar snow and firn in the Community Earth System Model, J. Adv. Model. Earth Sy., 9, 2583–2600,, 2017. a

Woodhouse, C. A.: A 431-yr reconstruction of western Colorado snowpack from tree rings, J. Climate, 16, 1551–1561,<1551:AYROWC>2.0.CO;2, 2003. a

Yakir, D.: Variations in the natural abundance of oxygen-18 and deuterium in plant carbohydrates, Plant Cell Environ., 15, 1005–1020,, 1992. a, b

Yoshimura, K., Kanamitsu, M., Noone, D., and Oki, T.: Historical isotope simulation using reanalysis atmospheric data, J. Geophys. Res.-Atmos., 113, D19108,, 2008. a

Zhu, H., Huang, R., Asad, F., Liang, E., Bräuning, A., Zhang, X., Dawadi, B., Man, W., and Grießinger, J.: Unexpected climate variability inferred from a 380 year tree-ring earlywood oxygen isotope record in the Karakoram, Northern Pakistan, Clim. Dynam., 57, 701–715,, 2021. a

Short summary
We modify the numerical model of forest growth MAIDENiso by explicitly simulating snow. This allows us to use the model in boreal environments, where snow is dominant. We tested the performance of the model before and after adding snow, using it at two Canadian sites to simulate tree-ring isotopes and comparing with local observations. We found that modelling snow improves significantly the simulation of the hydrological cycle, the plausibility of the model and the simulated isotopes.