A new snow module improves predictions of the isotope-enabled MAIDENiso forest growth model
- 1Centre de Recherche sur la dynamique du système Terre (GEOTOP), Université du Québec Montréal (UQAM), Montréal, Quebec, H2X 3R9, Canada
- 2Centre d'études nordiques (CEN), Université de Laval, Québec City, Quebec, G1V 0A6, Canada
- 3Institut de Recherche sur les Forêts (IRF), Université du Québec en Abitibi-Témiscamingue (UQAT), Amos, Quebec, J9T 2L8, Canada
- 4Carbon Cycle Research Group, Space and Atmospheric Physics, Physics Department, Imperial College London, London, SW7 2AZ, United Kingdom
- 5NASA Goddard Institute for Space Studies, Applied Physics and Applied Mathematics, Columbia University, New York, NY, USA
- 6Tree-Ring Laboratory, Lamont-Doherty Earth Observatory, Columbia University, Palisades, NY, 10964, USA
- 7Ecological and Forestry Applications Research Centre (CREAF), Bellaterra (Cerdanyola del Vallés), Barcelona, Spain
- 8Catalan Institution for Research and Advanced Studies (ICREA), Pg. Lluís Companys 23, Barcelona, Spain
Correspondence: Ignacio Hermoso de Mendoza (email@example.com)
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.
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égin, 1997). Snowmelt contributes to mitigate the 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 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; Woodhouse, 2003; 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 (Misson, 2004), 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 (Misson, 2004) 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 (Misson, 2004; 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.1 MAIDENiso model
2.1.1 Original model
MAIDENiso (Misson, 2004; 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.
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).
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 Wolock, 2009). 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:
where Δ () is the gradient of the saturation vapour pressure curve, γ () 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 Black, 2004). 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.
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 (). 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.
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 , or the GPP parameters control the amount of carbon produced, which in turn affects . 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 .
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 . 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 . Two different approaches can be used to infer . The first and most direct way is to use daily values of 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 from an observed dataset to obtain a linear regression model for daily values of based on temperature and precipitation:
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 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 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).
2.4 Tree-ring δ18O, GPP and snow data
We used published 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; http://fluxnet.ornl.gov/site/269 (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 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 Sutcliffe, 1970), which is equivalent to a coefficient of determination:
where and 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 Bauwens, 2011):
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 (Knoben et al., 2019).
To estimate the effect of the new snow module on predictions of , we compared the parameters influencing obtained by independent calibrations.
We also investigated the relative contributions to 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 from the reference simulations with those obtained from two experiments. First, to isolate the contribution of the source water on , we set the relative humidity (hair) and constant using the average values of hair and obtained from the reference simulations. Second, to isolate the contribution of the isotopic enrichment of the leaf water during transpiration on , we set δ18O in xylem water () 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.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.
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.
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).
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.
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.
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 series (Fig. 6). This was expected because (1) the calibration process maximized the coincidence between observed and simulated , and (2) the biochemical and kinetic fractionation parameters ϵ0 and ϵk could compensate the mean level of in Eq. (B1) for differences in . Regarding the agreement between inter-annual variations, the observed and simulated 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 , which makes it difficult to appreciate the improvement in the correlation between simulated and observed 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 series (transformed to mean 0 and standard deviation 1) in Fig. 6c and d.
The distribution of the optimized parameters controlling 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 was stronger than without considering 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 (, 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, (Fig. 8b). Without the snow module, followed closely the signal shown in Fig. 8a, with a small delay due to the isotopic mixing in the soil. The signal was slightly enriched with respect to 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 values were higher in the version with snow. The difference in values between the two versions reduced in time due to the infiltration and mixing of the enriched summer precipitation. Isotopic composition of discharge water (, Fig. 8c) follows a similar pattern to xylem water but with a delay of about 2 months and lower amplitudes in the seasonal variations.
The daily δ18O in the leaf () 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 level. The mean level of did not change after adding snow, despite the enrichment of and , 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 signature
Finally, we investigated the relative contributions from the source water through and the leaf transpiration enrichment through on the time series (Fig. 9, peak values in Table 5). At both sites, the leaf water δ18O isotopic enrichment had a stronger influence on than the δ18O variability of xylem source water, as shown by the higher variance explained by the experiment that simulated considering only the effect of δ18O leaf water enrichment indicated by a higher coefficient of determination (R2).
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 .
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 Black, 2004).
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 Running, 2006; 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 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 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 . 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 was the best option for this particular study. The model without snow produced depleted values of , because it lacked the enrichment effect of evaporative fractionation in snow. This was also reflected in the 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, 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 , 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 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 Ehleringer, 2000; Saurer et al., 1997; Sternberg et al., 1986; Yakir, 1992). 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 .
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 and/or in Eq. (B1) to adjust the mean 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 than the isotopic composition of the source (xylem) water, both in Tungsten and Caniapiscau, suggesting that it is the main driver of variations. These results are in agreement with Lavergne et al. (2017) findings and reflect the strong effect of vapour pressure deficit on in Quebec. Nevertheless, the signature also had a strong imprint of the source water signal as recently reported for the Tungsten record that shared the same large-scale atmospheric patterns than spring–summer (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 values during the growing season (Fig. 8b). As our results have shown, this significantly increased the correlation between the observed and simulated 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 indirectly through the source water.
Overall, the improvements found in the simulations at both sites indicate that snow plays a critical role in δ18O of the source water and thus on the final signature of . Even if the addition of snow would not had resulted in a significant improvement of the correlation between the simulated and observed , 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 enrichment signal on 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 studies that did not take into account the additional effects of snowpack dynamics on , 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 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.
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 Linder, 2000). The simulations of the new version of MAIDENiso reproduce the observed better than the original snowless version of MAIDENiso. Based on the development presented here, the potential for the application of MAIDENiso is notably increased.
GPP () 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
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.
where g0=0 and g1=10 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:
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
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.
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 (). This is based on the Danis et al. (2012) formulation of the Craig–Gordon model (Craig and Gordon, 1965):
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 Ehleringer, 2000; Saurer et al., 1997; Sternberg et al., 1986; Yakir, 1992). ϵ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 Epstein, 1979; 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 Richter, 2008). ϵ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). and are the δ18O of vapour and xylem (source) water, respectively. is calculated from the δ18O of precipitation () and the fractionation due to the phase change from liquid water to vapour at mean air temperature, (Horita and Wesolowski, 1994):
The time series produced through Eq. (B1) are daily, while the 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):
The MAIDENiso code is available in the Zenodo repository. Note that there are two versions of the code, corresponding to the model with snow (https://doi.org/10.5281/zenodo.5597877, Hermoso de Mendoza, 2021a) and the model without snow (https://doi.org/10.5281/zenodo.5598076, Hermoso de Mendoza, 2021b). In addition, a university website has been created (https://dendro-eco.uqat.ca/maiden/, 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 (https://doi.org/10.5281/zenodo.5599091, Hermoso de Mendoza, 2021c).
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 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.
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.
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).
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, https://doi.org/10.1002/2016GL069690, 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, https://doi.org/10.1038/NCLIMATE2246, 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, https://doi.org/10.1002/wat2.1311, 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, https://doi.org/10.1111/j.1365-3040.2001.00668.x, 2001. a
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, https://doi.org/10.1111/j.1365-2486.2006.01134.x, 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, https://doi.org/10.5194/bg-11-3245-2014, 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, https://doi.org/10.1038/s41598-017-02580-9, 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, https://doi.org/10.1016/0016-7037(96)00073-7, 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, https://doi.org/10.1093/treephys/tpu115, 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, https://doi.org/10.1088/1748-9326/abd5de, 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, https://doi.org/10.1139/x2012-089, 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, https://doi.org/10.1111/j.1365-3040.1997.00094.x, 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, https://doi.org/10.1126/science.204.4388.51, 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, https://doi.org/10.1038/s41467-018-05705-4, 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, https://doi.org/10.1016/j.jtbi.2010.02.013, 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, https://doi.org/10.1080/01431161.2014.955146, 2014. 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, https://doi.org/10.1007/978-1-4612-3498-2_2, 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
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, https://doi.org/10.1007/s00382-021-05932-4, 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, https://doi.org/10.3354/cr001097, 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, https://doi.org/10.1016/j.envsoft.2008.03.003, 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, https://doi.org/10.5194/bg-12-3695-2015, 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, https://doi.org/10.5194/bg-14-4851-2017, 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, https://doi.org/10.1007/s00382-017-3565-5, 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, https://doi.org/10.1029/2010JG001287, 2008. a
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, https://doi.org/10.1016/j.jhydrol.2009.08.003, 2009. 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], https://doi.org/10.5281/zenodo.5599091, 2021c. 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, https://doi.org/10.1016/0016-7037(94)90096-5, 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, https://doi.org/10.1016/j.epsl.2018.11.002, 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, https://doi.org/10.1175/2008JAMC1979.1, 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, https://doi.org/10.3402/tellusa.v28i6.11322, 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, https://doi.org/10.1002/2015jg003114, 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, https://doi.org/10.5194/hess-23-4323-2019, 2019. 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, https://doi.org/10.5194/cp-13-1515-2017, 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, https://doi.org/10.1029/2018MS001583, 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, https://doi.org/10.1111/j.1365-3040.1995.tb00628.x, 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, https://doi.org/10.1002/2017GL073551, 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, https://doi.org/10.5194/bg-11-6711-2014, 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, https://doi.org/10.1175/2009JHM1148.1, 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, https://doi.org/10.1175/2009EI283.1, 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, https://doi.org/10.1175/BAMS-87-3-343, 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, https://doi.org/10.1175/1520-0442(1994)007<0914:CSRFIC>2.0.CO;2, 1994. a
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, https://doi.org/10.5194/cp-11-1153-2015, 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, https://doi.org/10.1016/j.jhydrol.2014.03.054, 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, https://doi.org/10.1139/cjfr-2014-0080, 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, https://doi.org/10.2151/jmsj1965.77.1B_235, 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, https://doi.org/10.1139/x26-241, 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, https://doi.org/10.1073/pnas.0708986104, 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, https://doi.org/10.1139/x94-247, 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, https://doi.org/10.1007/s00382-013-1674-3, 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, https://doi.org/10.1016/S0304-3800(03)00161-3, 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, https://doi.org/10.1007/s004420000349, 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, https://doi.org/10.1139/x87-081, 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, https://doi.org/10.1029/2008JD011063, 2009. a
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, https://doi.org/10.1016/j.dendro.2015.10.001, 2016. a, b
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, https://doi.org/10.1175/2008JCLI2441.1, 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, https://doi.org/10.1104/pp.82.2.423, 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, https://doi.org/10.3389/feart.2018.00108, 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, https://doi.org/10.1111/gcb.14395, 2018. 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, https://doi.org/10.1007/3-540-31298-6, 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, https://doi.org/10.1002/2017MS000988, 2017. a
Woodhouse, C. A.: A 431-yr reconstruction of western Colorado snowpack from tree rings, J. Climate, 16, 1551–1561, https://doi.org/10.1175/1520-0442(2003)016<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, https://doi.org/10.1111/j.1365-3040.1992.tb01652.x, 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, https://doi.org/10.1029/2008JD010074, 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, https://doi.org/10.1007/s00382-021-05736-6, 2021. a