Introduction 

Modern climatic changes are leading to significant air temperature increase in the Eurasian subarctic at even faster rates compared to the global average ([33], [14]). Siberian forests respond to both the timing and magnitude of changes in soil moisture as well as to soil temperature, permafrost thickness and distribution, which is in turn affected by snow and vegetation cover, soil texture and geothermal heat flux as well as atmospheric dryness ([5], [4], [8], [24]). Therefore, the fate of trees growing on permafrost undergoing climatic changes is of great interest, owing to the important role of permafrost in these ecosystems: (i) as additional water source ([43], [5], [6], [47]) for trees during droughts; and (ii) for the large amounts of carbon stored in these soils, which essentially contribute to the global carbon and water budget ([5]).

Increasing temperatures at higher latitudes, especially in areas with low precipitation inputs such as parts of Yakutia (200-250 mm year^-1), enhance the risks of tree decline and mortality in boreal forests ([39], [4], [6], [44]). A number of studies report a pronounced increase of seasonal permafrost thaw depth in western Siberia ([43], [30], [34], [16]). Sugimoto et al. ([43]) reported, that the active soil layer depth at the central Yakutia can reach up to 1.4 m making thawed permafrost water available for Larix gmelinii Rupr. Moreover, the authors demonstrated that during wet summers in central Yakutia Larix gmelinii Rupr. can use rain water, while thawed permafrost water can be used by trees mainly during drought summers. Fyodorov-Davydov et al. ([16]) investigated the spatial and temporal trends in active soil layer for northern Yakutia and its relation to landscape and climate variables. The soil layer thaws during June-August and freezes again in autumn, defying the length of the short growing season. The few studies in Siberian regions show significant differences between sites in hydrological and geomorphological characteristics, which should be taken into account. To our knowledge, there are no direct seasonal or long-term measurements of the thawing depth of the active soil layer for northeastern Yakutia, which is the driest forested area in the world. Precipitation and thawed permafrost water are particularly crucial for trees growing in this region with severe temperature limitations and low amounts of precipitation ([2], [40], [39]).

Based on isotopic differences between precipitation and permafrost, variations of stable oxygen isotopes (δ¹⁸O) may help to reveal water relations in deep-rooted plants ([37]). Oxygen isotopes in plant organic matter are influenced by variations in the isotopic composition of the water source, which is closely related to the precipitation and soil water (after modification by evaporation at the soil surface and due to the shallow active layer depth).

In this study, we aim to determine the source water of Larix cajanderi Mayr. and Pinus sylvestris L. trees from the Siberian subarctic using measured and modeled δ¹⁸O variation in annual tree-ring cellulose at two contrasting study sites (northeastern and central Yakutia).

  Materials and Methods 

Study sites

Both study sites are located in the broad permafrost zone in northeastern Siberia: in northeastern Yakutia (NE - 69° 22′ N, 148° 25′ E, ~ 250 m a.s.l.) and Central Yakutia (CE - 62° 14′ N, 129° 37′ E, ~ 220 m a.s.l.) in Russia.

The mean winter temperatures (December-February) averaged for the period from 1948 to 2004 are -36 °C and -38 °C and in summer (June-July) +12 °C and +19 °C for NE and CE, respectively. However, minimum temperature range during winter can vary between -38 °C and -46 °C, while maximum day time summer temperatures range from +27 °C to +42 °C. Larix cajanderi Mayr. and Pinus sylvestris L. trees are the dominant tree species in NE and CE respectively. The growing period is rather short for trees at both sites and lasts around 60 days ([3]).

Low precipitation inputs are typical for both regions (200-250 mm year^-1) and most of the precipitation falls in summer. The main difference between the two sites are soil conditions. Larix cajanderi Mayr. trees in NE are mainly growing on stones, rocks and mixtures of clay ([40]), while Pinus sylvestris L. in CE are growing on mixtures of sand and clay soils ([32]). The root system is also different between these species as they adapt to site-specific edaphic conditions. Larix cajanderi Mayr. roots are shallower than those of Pinus sylvestris L., which tend to be slightly deeper and more extensive in spread.

Local climate data

To assess the drivers of δ¹⁸O variation in tree-ring cellulose for measured and modeled datasets, we used daily maximum and minimum temperature, daily precipitation and daily relative humidity data from the Chokurdach weather station (70° 30′ N, 148° 08′ E, ~ 250-300 m a.s.l.) located close to the north-eastern site (NE) for the period from 1948 to 2004. For Central Yakutia (CE) we used data from the weather station in Yakutsk (62° 02′ N, 129° 44′ E, ~ 220 m a.s.l.) for the period from 1966 to 2004.

Annual δ¹⁸O in tree-ring cellulose chronologies

As a reference, we used δ¹⁸O values of Larix cajanderi Mayr. cellulose chronologies, which were constructed based on individually analysed annual tree rings from four trees from the NE site ([40]) for the period from 1880 to 2004 and from five trees for the period from 1900 to 2013 from the CE sites ([44]). These measured δ¹⁸O chronologies were used for comparative analyses with the LPX-Bern modeled data of the source water, δ¹⁸O of cellulose, and thaw depths for the common period (from 1901 to 2004) for both sites.

Cellulose samples from both study sites (NE and CE) were measured at the Stable Isotope Research Facility of the Paul Scherrer Institute (Villigen, Switzerland). For measurements of the oxygen isotope ratios from NE, cellulose samples were pyrolysed to CO at 1080 °C ([38]) in the continuous flow mode with an elemental analyzer (Carlo Erba EA-1110^®, Milan, Italy) linked to an isotope-ratio mass-spectrometer (Delta S^®, Finnigan MAT, Bremen, Germany). The CE samples were pyrolyzed using PYRO-cube^® at 1420 °C (Elementar, Hanau, Germany) and analyzed with a Delta Plus XL^® Mass spectrometer (Thermo Finnigan, Bremen, Germany). An earlier study showed that both systems yielded similar precisions (± 0.2‰.) and the values from the two instruments were highly correlated (better than the measurement precisions - [49]). The isotopic values were expressed in the δ notation relative to the international standard (eqn. 1):

where R_sample is the molar fraction of ¹⁸O/¹⁶O ratio of the sample and R_standard of the Vienna Standard Mean Ocean Water for oxygen.

Estimation of the δ¹⁸O of source water using RLE model

Roden-Lin-Ehleringer (RLE) model was applied to determinate the δ¹⁸O of source water over time. For snowmelt and thawed permafrost water as well as rainwater we employed ranges between -28 and -25‰ and -6 to -11‰, respectively ([25], Prokushkin, personal communication). An earlier study by Saurer et al. ([37]) showed that soil water is a mixture of precipitation and thawed permafrost water, and that evaporation plays a minor role in this region. Barometric pressure was obtained as an estimate of 95 kPa from elevation tables in Pearcy et al. ([35]). Fan et al. ([12]) reported that values for stomatal conductance in leaves of Salix boganidensis from lowland Indigirka site in the northwestern Yakutia could range between 0.4 and 0.6 mol m^-2 s^-1. Streit et al. ([42]) reported that the stomatal conductance of Larix decidua needles was higher compared to Pinus mugo and varied from 0.2 to 0.6 mol m^-2 s^-1 during July at the high-elevated site in Stillberg, Davos (2180 m a.s.l.), Switzerland. To our knowledge, there are no direct measurements for Larix cajanderi Mayr. or Pinus sylvestris for Yakutia. For stomatal conductance we used value of 0.3 mol m^-2 s^-1, assuming a low stomatal conductance for conifer needles from a northern tree line in the cold Yakutia region.

The δ¹⁸O of xylem cellulose (δ¹⁸O_cx) was predicted as a function of both the isotopic ratio of the substrate sucrose and of the medium water at the site of cellulose synthesis, root, or xylem as part of the tree-ring ([36]). The δ¹⁸O in tree-ring cellulose is calculated as (eqn. 2):

where f₀ is the fraction of exchange with leaf water (δ¹⁸O_wx) and estimated to be 0.42 for oxygen from the best fit relationship derived from controlled experiments ([36]). The subscript cx is xylem cellulose, wx is xylem water and lw is leaf water at the site of sucrose synthesis. The parameter ε₀ is a biochemical fractionation associated with the water-carbonyl oxygen interactions (27 ‰ - [41], [27]).

We modeled leaf water as (eqn. 3):

where δ¹⁸O_lw is the isotope ratio at the site of evaporation, δ¹⁸O_source is a mixture of soil and precipitation water used by trees, ε_e = 1.0096, which is the temperature dependent equilibrium fractionation due to the phase change from liquid water to vapour ([29]), ε_k = 1.0285, which is the kinetic fractionation ([31]) due to diffusion of vapour into unsaturated air, δ¹⁸O_vapour is the isotopic composition of atmospheric water vapour, and e_a/e_i the ratio of the ambient to the leaf internal vapour pressure (kPa - [9], [10]).

The RLE model can predict source water δ¹⁸O values for both the NE and CE study sites by the following equation (eqn. 4):

where ƒ_{0 =} 0.42, which is fraction of exchange with medium water; ε_oh and ε_oa = 27, which is heterotrophic and autotrophic fractionation factor during cellulose synthesis; δ¹⁸O_lwm is the modeled leaf water ([46]).

Simulating δ¹⁸O in LPX Bern

The Land surface Processes and eXchanges (LPX-Bern v. 1.3) model ([22]) is a dynamic global vegetation model and features a process-based representation of the coupled terrestrial nitrogen, carbon and water cycle. Here only a natural vegetation, consisting of ten Plant Functional Types (PFTs) competing for resources and confined by bioclimatic limits, is considered. The δ¹⁸O isotope enabled version of the LPX-Bern model is introduced and discussed in detail by Keel et al. ([21]), which helps to reduce uncertainties in the interpretation of raw measured δ¹⁸O chronologies for investigating different temporal scales of the oxygen isotope patterns in tree rings. The δ¹⁸O_cell (stem cellulose isotopic composition) was calculated based on the eqn. 4.

The δ¹⁸O_lw is calculated according to eqn. 3, but modified by including the Péclet effect ([13]) using the modeled rate of transpiration (eqn. 5):

where P = EL/CD; L is mixing path length (m); E is the rate of transpiration (mol m^-2 s^-1); C is the concentration of water 5.55 · 10⁴ mol m^-3; D is the diffusivity of H₂¹⁸O in water 2.66 · 10^-9 m² s^-1; P is the Péclet number ([20]); δ¹⁸O_{lw Péclet} is the Péclet corrected leaf water enrichment; δ¹⁸O_lw from the eqn. 3 is the enrichment of the leaf water at the surface of evaporation.

In this study, the model runs for both study sites with a single grid cell configuration from 1901 to 2004. Monthly temperature, precipitation, cloud cover and number of wet days data are extracted from the 0.5° latitude/longitude grid cells closest to the study cells locations of the global CRU-TS3.22 climate dataset ([19]). Internally, the meteorological data are interpolated to a daily time-step, except for precipitation, where a stochastic weather generator is applied to compute daily precipitation according to Gerten et al. ([17]). The temperature and precipitation data correspond well (r = 0.96; p < 0.05) with the local weather stations (see above) from 1966 to 2004. Similarly, nitrogen deposition is supplied to the model ([26]). Furthermore, the oxygen module requires input values for the δ¹⁸O of soil water and water vapor, as well as relative humidity data, which are extracted from global transient isotope-enabled simulations of the coupled atmosphere-land surface model ECHAM5-JSBACH ([18]). Since the simulation only spans the period 1960-2012, earlier years are approximated by the values of the year 1960. Identical atmospheric CO₂ concentrations are prescribed to both study sties from a global reconstruction ([11], [15], [28]).

The LPX-Bern δ¹⁸O simulations are useful for investigating different temporal scales of the oxygen isotope patterns in tree-ring, which help to reduce uncertainties in the interpretation of raw measured δ¹⁸O chronologies.

Climatological and statistical analyses

Pearson correlation coefficients were calculated between local weather station data and the CRU-TS3.22 climate dataset were calculated ([19]). Statistical calculations were performed in licensed version of the software Statistica^® v. 13.3 (StatSoft, Tulsa, OK, USA).

  Results 

Modeled and measured δ¹⁸O tree-ring cellulose data

Using the LPX-Bern ([21]) and RLE ([36]) models, we obtained modeled δ¹⁸O tree-ring cellulose values for the NE (Fig. 1a) and CE (Fig. 1b) study sites in Yakutia.

Fig. 1 - The δ¹⁸O in tree-ring cellulose measured and modeled using LPX-Bern and RLE models. Modeled thaw depth as output parameter of the LPX model is presented for the (a) northeastern (NE) and (b) Central (CE) sites in Yakutia. Maximal seasonal thaw depth (LPX-TD_max) simulated by LPX-Bern model is presented in black lines for both sites (a, b).

  Enlarge/Shrink   Download   Full Width  Open in Viewer

Modeled LPX-Bern thaw soil depth varied much stronger for the NE site (in the range from 50-100 cm) compared to the modeled soil depth at the CE site (50-70 cm). Modeled values and data for δ¹⁸O in cellulose agree in so far that the tree-ring cellulose is more depleted in the heaver ¹⁸O isotope at the NE site than at the CE.

Measured δ¹⁸O data from the NE and CE are significantly correlated with RLE (r = 0.92 and r = 0.74, p < 0.05) and LPX-Bern (r = 0.59 and r = 0.46, p < 0.05), respectively for both studied regions. On average the offset between modeled RLE and measured δ¹⁸O data for both NE and CE is low (less than 1‰).

The LPX-Bern modeled thaw depth (maximum values) correlates significantly (r = -0.29, p < 0.05) with measured NE δ¹⁸O data only. The δ¹⁸O measured data from the NE and CE correlated significantly between each other (r = 0.20, p < 0.05) for the common period from 1901 to 2004, although the distance between the two sites is about 860 km. Modeled data for thaw permafrost soil depth at NE and CE also correlate significantly (r = 0.39, p < 0.05).

Modeling of source water

The δ¹⁸O leaf water in the RLE is calculated based on the isotopic composition of leaf water at the site of evaporation, which is not the same as bulk leaf water due to the contribution of unevaporative vein water and tissues ([45]) and the “Péclet effect” (see [13] which describe the convective / diffusive mixing between evaporatively enriched and un-enriched source water from the xylem).

A strong increase in estimated isotopic composition of the source water after the 1990s for NE was revealed, while an overall declining trend for CE was found (Fig. 2). Correlation coefficients between source water δ¹⁸O and July precipitation were significant at NE (r = -0.27, p < 0.05) and for May precipitation (r = -0.28, p < 0.05) at the CE site. However, the negative correlation with July precipitation was weak.

Fig. 2 - RLE modeled source water δ¹⁸O for (a) larch (red line) vs. July precipitation (black bars) and (b) pine trees (blue line) vs. May precipitation (blue bars).

  Enlarge/Shrink   Download   Full Width  Open in Viewer

A decrease in July precipitation at the NE site was detected after the year 2000, which is consistent with the modeled δ¹⁸O source water values.

  Discussion and conclusion 

Simulations of thawing depth and δ¹⁸O values in tree-ring cellulose, using the output of the RLE model as a proxy for source water over time, indicate that thawed permafrost water is available for trees from NE in addition to precipitation water. Based on this analysis we conclude that thawed permafrost water could be an important water source for trees, especially for periods when the amount of precipitation is lower than the evaporative water loss. Our model simulation is in line with earlier findings by Sugimoto et al. ([43]) for the central part of Yakutia, confirming that trees from the CE could use thawed permafrost water during dry events. Yet the question remains whether Larix cajanderi Mayr. from the NE site can utilize thawed permafrost water. There was also an indication of source water enrichment (¹⁸O) increased over time. This could be caused by rising temperatures resulting in higher evapotranspiration due to higher vapor pressure deficit (VPD) during drought conditions ([48]), but may also be the result of a diminishing influence of thawed permafrost water, which has a more negative isotope ratio than rain water ([37]).

Precipitation or snowmelt can be stored in soils or lost as run off, depending on permafrost depth. Source water for Pinus sylvestris L. trees from the CE site did not indicate a strong increase of thawed permafrost use as an additional water source as it was found at the NE. Opposite to the NE, trees from the CE have limited water availability caused by drought conditions in this area.

Our study showed that trees from NE and CE might have access to additional thawed permafrost water during dry summer periods, which is in line with the study by [43]. Due to differences in the soil structure, active thaw soil depth and root systems of trees at two Siberian sites, Larix cajanderi Mayr. trees cannot access water deeper than 50 cm, while Pinus sylvestris L. in Central Yakutia can access water as deep as 80 cm. This is opposite to the finding of Sugimoto et al. ([43]), who showed possible water access up to 1.4 m. These differences can be explained by different soil structures (sand in central Yakutia and rocks in northeastern Yakutia) as well as with species-specific rooting distributions, resulting in different rooting depths. Thus, trees have different access to water sources depending on their root distribution, i.e., surface rooting system (Larix cajanderi Mayr.) vs. roots proliferated to deeper soil layers (Pinus sylvestris L. - [40], [23], [7], [37], [24]).

We could show that the application of different oxygen isotope models ([21], [36]) helps to reveal the impact of important hydro-climatic factors such as the interplay between water sources; e.g., precipitation, soil moisture and frozen soil water and soil thawing depth, influencing tree-growth and physiology in an ecosystem under permafrost conditions. The discrepancies between modeled and measured data are in part due to different parameterizations, but more so due to different model designs for different scales, from whole plant to ecosystems. So far, direct seasonal measurement for the thawed permafrost depth or measurement for stomatal conductance at the remote study site, like NE are rare due to the harsh climatic conditions leading to costly and complex fieldwork at this region. Therefore, the modeling approach can be a powerful tool for the investigation of the influence of a variety of climatic factors on Siberian forest ecosystem water relations that impact on isotope fractionations in tree-ring cellulose from remote Siberian sites.

Furthermore, we must keep in mind that the heterogeneous soil structure will result in a great variation of δ¹⁸O values of the xylem water, especially with regard to time gaps between the isotopic signals of precipitation and what we find in xylem water. Depending on the depth of the active soil layer (where trees can access water) and its structure, these time gaps can be considerable ([1]). Such structural variations are difficult to take into account in any modeling approach. Further studies on adjustments and harmonization between LPX-Bern and RLE models and their parameterization based on a wider pool of measured data will provide the basis to enhance the relevance of these models to a powerful analytical tool. Although each model has its strengths and weaknesses, making use of their complementarity and strengths facilitates the enhancement of our understanding of the source water availability and its the complexity for trees in these boreal ecosystems under permafrost dynamics. For instance, the use of a process-based model, such as the LPX-Bern ([21]), allows investigating δ¹⁸O signals in tree rings as a proxy for hydrological and climate dynamics on a regional or global scale. Furthermore, factorial simulations (where each parameter is modified specifically) are a potential avenue to further improve our understanding of the linkage between the vegetation cover, climatic drivers and permafrost dynamics.

  Acknowledgments 

This work was supported by the Russian Science Foundation (RSF) (19-14188000 28), EU-FP7 Marie Curie Fellowship (EU-ISOTREC 235122, 909122). FJ and SL acknowledge support by the Swiss National Science Foundation (# 200020_172476). We thank Martin Werner for providing the soil water and humidity δ¹⁸O data from the ECHAM5-JSBACH model.

  Author Contributions 

OCh did the measurements and data interpretations for the NE Yakutia, and applied the RLE model to both studied sites. GT performed the experiments at the CE Yakutia site. SL modeled δ¹⁸O and thaw depth for both NE and CE sites. All authors significantly contributed to the writing of the manuscript.

  Conflicts of Interest 

The authors declare no conflict of interest.

  References