Partitioning of canopy and soil CO₂ fluxes in a pine forest at the dry timberline across a 13-year observation period

Partitioning carbon fluxes is key to understanding the process underlying ecosystem response to change. This study used soil and canopy fluxes with stable isotopes (¹³C) and radiocarbon (¹⁴C) measurements in an 18 km², 50-year-old, dry (287 mm mean annual precipitation; nonirrigated) Pinus halepensis forest plantation in Israel to partition the net ecosystem's CO₂ flux into gross primary productivity (GPP) and ecosystem respiration (R_e) and (with the aid of isotopic measurements) soil respiration flux (R_s) into autotrophic (R_sa), heterotrophic (R_h), and inorganic (R_i) components. On an annual scale, GPP and R_e were 655 and 488 g C m⁻², respectively, with a net primary productivity (NPP) of 282 g C m⁻² and carbon-use efficiency (CUE = NPP ∕ GPP) of 0.43. R_s made up 60 % of the R_e and comprised 24±4 %R_sa, 23±4 %R_h, and 13±1 %R_i. The contribution of root and microbial respiration to R_e increased during high productivity periods, and inorganic sources were more significant components when the soil water content was low. Comparing the ratio of the respiration components to R_e of our mean 2016 values to those of 2003 (mean for 2001–2006) at the same site indicated a decrease in the autotrophic components (roots, foliage, and wood) by about −13 % and an increase in the heterotrophic component (R_h∕R_e) by about +18 %, with similar trends for soil respiration (R_sa∕R_s decreasing by −19 % and R_h∕R_s increasing by +8 %, respectively). The soil respiration sensitivity to temperature (Q₁₀) decreased across the same observation period by 36 % and 9 % in the wet and dry periods, respectively. Low rates of soil carbon loss combined with relatively high belowground carbon allocation (i.e., 38 % of canopy CO₂ uptake) and low sensitivity to temperature help explain the high soil organic carbon accumulation and the relatively high ecosystem CUE of the dry forest.

1 Introduction

The annual net storage of carbon in the land biosphere, known as net ecosystem production (NEP), is the balance between carbon uptake during gross primary productivity (GPP) and carbon loss during growth, maintenance respiration by plants (i.e., autotrophic respiration, R_a), and decomposition of litter and soil organic matter (i.e., heterotrophic respiration, R_h; Bonan, 2008). The difference between GPP and R_a expresses the net primary production (NPP) and is the net carbon uptake by plants that can be used for new biomass production. Measurements from a range of ecosystems have shown that total plant respiration can be as large as 50 % of GPP (e.g., Etzold et al., 2011) and together with R_h comprises total ecosystem respiration (R_e; $R_{\mathrm{e}}=R_{\mathrm{a}}+R_{\mathrm{h}}$). The partitioning of the ecosystem carbon fluxes can therefore be summarized as

$$\begin{array}{}\text{(1)}& \mathrm{GPP}=\mathrm{NPP}+R_{\mathrm{a}}=\mathrm{NEP}+R_{\mathrm{h}}+R_{\mathrm{a}}.\end{array}$$

Earlier campaign-based measurements carried out by Maseyk et al. (2008a) and Grünzweig et al. (2009) in the semiarid Pinus halepensis (Aleppo pine) Yatir Forest indicated that GPP at this site was lower than among temperate coniferous forests (1000–1900 g C m⁻² yr⁻¹) but within the range estimated for Mediterranean evergreen needleleaf and boreal coniferous forests (Falge et al., 2002; Flechard et al., 2019b) and had a high carbon-use efficiency (CUE) of 0.4 (CUE = NPP ∕ GPP; DeLucia et al., 2007). The total flux of CO₂ released from the ecosystem (R_e) can be partitioned into aboveground autotrophic respiration (i.e., foliage and sapwood, R_f) and soil CO₂ flux (R_s). R_s, in turn, is a combination of three principal components and can be further partitioned into the components originating from roots or rhizospheres and mycorrhizas (i.e., R_sa), from carbon respired during the decomposition of dead organic matter by soil microorganisms and macrofauna (R_h; Bahn et al., 2010; Kuzyakov, 2006), and from pedogenic or anthropogenic acidification of soils containing CaCO₃ (R_i; Joseph et al., 2019; Kuzyakov, 2006), which is expressed as

$$\begin{array}{}\text{(2)}& R_{\mathrm{e}}=R_{\mathrm{s}}+R_{\mathrm{f}}=\left[R_{\mathrm{sa}}+R_{\mathrm{h}}+R_{\mathrm{i}}\right]+R_{\mathrm{f}}.\end{array}$$

Previously published results show that the contribution of R_sa and R_h to R_s ranges from 24 % to 65 % and from 29 % to 74 %, respectively, in forest soils in different biomes and ecosystems (Binkley et al., 2006; Chen et al., 2010; Flechard et al., 2019a; Frey et al., 2006; Hogberg et al., 2009; Subke et al., 2011). Some studies reported significant proportions of abiotic contribution to R_s, ranging between 10 % and 60 % (Martí-Roura et al., 2019; Ramnarine et al., 2012; Joseph et al., 2019). However, most of these experiments were performed in boreal, temperate, or subtropical forests, and there is a general lack of information on water-limited environments, such as dry Mediterranean ecosystems. Using both ¹³C and CO₂ ∕ O₂ ratios also showed that abiotic processes, such as CO₂ storage, transport, and interactions with sediments, can influence R_s measurements at such sites (Angert et al., 2015; Carmi et al., 2013). Furthermore, root-respired CO₂ can also be dissolved in the xylem water and carried upward with the transpiration stream (Etzold et al., 2013).

Rates of the soil–atmosphere CO₂ flux (R_s) have been altered owing to global climatic change, particularly through changes in soil temperature (T_s) and soil moisture (SWC; Bond-Lamberty and Thomson, 2010; Buchmann, 2000; Carvalhais et al., 2014; Hagedorn et al., 2016; Zhou et al., 2009), which could account for 65 %–92 % of the variability of R_s in a mixed deciduous forest (Peterjohn et al., 1994). Soil moisture impacts on R_s have been observed in arid and Mediterranean ecosystems, where T_s and SWC are negatively correlated (e.g., Grünzweig et al., 2009). CO₂ efflux generally increases with increasing soil temperatures (Frank et al., 2002), which can produce positive feedback on climate warming (Conant et al., 1998), converting the biosphere from a net carbon sink to a carbon source (IPCC, 2014). A range of empirical models have been developed to relate R_s rate and temperature (Balogh et al., 2011; Lellei-Kovács et al., 2011), and the most widely used models rely on the Q₁₀ approach (Bond-Lamberty and Thomson, 2010), which quantifies the sensitivity of R_s to temperature and can integrate it with physical processes, such as the rate of O₂ diffusion into and CO₂ diffusion out of soils and the intrinsic temperature dependency of enzymatic processes (Davidson and Janssens, 2006). Soil moisture (SWC) may be of greater importance than temperature in influencing R_s in water-limited ecosystems (Hagedorn et al., 2016; Grünzweig et al., 2009; Shen et al., 2008). In general, the R_s rate increases with the increase of SWC at low levels but decreases at high levels of SWC (Deng et al., 2012; Hui and Luo, 2004; Jiang et al., 2013). Several studies highlight the sensitivity of carbon fluxes in semiarid Mediterranean ecosystems to the irregular seasonal and interannual distribution of rain events (Poulter et al., 2014; Ross et al., 2012). While R_s is generally constrained by low SWC during summer months, abrupt and large soil CO₂ pulses have been observed after rewetting the dry soil (Matteucci et al., 2015).

The objectives were twofold: first, to obtain detail on partitioning of the carbon fluxes in a semiarid pine forest to help explain the high productivity and carbon use efficiency recently reported for this ecosystem (Qubaja et al., 2019) and provide process-based information to assess the carbon sequestration potential of such a semiarid afforestation system; and second, to combine this 2016 study with the results of a similar one at the same site in 2003 (mean values for 2001–2006; Grünzweig et al., 2007, 2009) to obtain a long-term perspective across 13 years on soil respiration and its partitioning. We hypothesized that the high carbon-use efficiency of the dry-forest ecosystem is associated with high belowground carbon allocation and relatively low decomposition rates and that the long-term trend associated with warming may be suppressed by the dry conditions.

2 Materials and methods

2.1 Site description

The Yatir Forest (31^∘20^′49^′′ N, 35^∘03^′07^′′ E; 650 m a.s.l.) is situated in the transition zone between subhumid and arid Mediterranean climates (Fig. S1 in the Supplement) on the edge of the Hebron mountain ridge. The ecosystem is a semiarid pine afforested area established in the 1960s and covering approximately 18 km². The average air temperatures for January and July are 10.0 and 25.8 ^∘C, respectively. Mean annual potential evapotranspiration (ET) is 1600 mm, and mean annual precipitation is 287 mm. Only winter (December to March) precipitation occurs in this region, creating a distinctive wet season, while summer (June to October) is an extended dry season. There are short transition periods between seasons, with a wetting season (i.e., fall) and a drying season (i.e., spring). The forest is dominated by Aleppo pine (Pinus halepensis Mill.), with smaller proportions of other pine species and cypress and little understory vegetation. Tree density in 2007 was 300 trees ha⁻¹; mean tree height was 10.0 m; diameter at breast height (DBH) was ∼15.9 cm, and the leaf area index (LAI) was ∼1.5. The native background vegetation was sparse shrubland, which is dominated by the dwarf shrub Sarcopoterium spinosum (L.) Spach, with patches of herbaceous annuals and perennials reaching a total vegetation height of 0.30–0.50 m (Grünzweig et al., 2003, 2007). The root density range is 30–80 roots m⁻² at the upper 0.1 m soil depth, falling to the minimum value (∼0 roots m⁻²) at 0.7 m soil depth (Preisler et al., 2019). Biological soil crust (BSC) is evident in the forest but is less than in the surrounding shrub by ∼40 % (Gelfand et al., 2012).

The soil at the research site is shallow (20–40 cm), reaching only 0.7–1.0 m; the stoniness fraction for the soil depth (0–1.2 m) is 15 %–60 %, and the rock cover of the surface ranges between 9 % and 37 %, as recently described in detail (Preisler et al., 2019); the soil is eolian-origin loess with a clay–loam texture (31 % sand, 41 % silt, and 28 % clay; density is 1.65±0.14 g cm⁻³) overlying chalk and limestone bedrock. Deeper soils (up to 1.5 m) are sporadically located at topographic hollows. While the natural rocky hill slopes in the region are known to create flash floods, the forested plantation reduces runoff dramatically to less than 5 % of annual rainfall (Shachnovich et al., 2008). Groundwater is deep (> 300 m), reducing the possibility of groundwater recharge due to negative hydraulic conductivity or of water uptake by trees from the groundwater.

2.2 Flux and meteorological measurements

An instrumented eddy covariance (EC) tower was erected in the geographical center of Yatir Forest, following the EUROFLUX methodology (Aubinet et al., 2000). The system uses a three-dimensional (3-D) sonic anemometer (Omnidirectional R3, Gill Instruments, Lymington, UK) and a closed path LI-7000 CO₂∕H₂O gas analyzer (LI-COR Inc., Lincoln, NE, USA) to measure the evapotranspiration flux (ET) and net CO₂ flux (NEE). EC flux measurements were used to estimate the annual scale of NEP by integrating half-hourly NEE values. The long-term operation of our EC measurement site (since 2000; see Rotenberg and Yakir, 2010) provides continuous flux and meteorological data with about 80 % coverage, which are subjected to U^* nighttime correction and quality control, and gap filling is based on the extent of the missing data, as recently described in more detail in Tatarinov et al. (2016). A site-specific algorithm was used for flux partitioning into R_e and GPP. Daytime ecosystem respiration (R_e−d; in µmol m⁻² s⁻¹) was estimated based on measured nighttime values (R_e−n; i.e., when the global radiation was < 5 W m⁻²), averaged for the first 3 half hours of each night. The daytime respiration for each half hour was calculated according to Eq. (3) (Maseyk et al., 2008a; Tatarinov et al., 2016):

$$\begin{array}{}\text{(3)}& R_{\mathrm{e}-d}=R_{\mathrm{e}-n}\left({\mathit{\alpha }}_{\mathrm{1}}{\mathit{\beta }}_{\mathrm{s}}^{{\mathrm{dT}}_{\mathrm{s}}}+{\mathit{\alpha }}_{\mathrm{2}}{\mathit{\beta }}_{\mathrm{w}}^{{\mathrm{dT}}_{\mathrm{a}}}+{\mathit{\alpha }}_{\mathrm{3}}{\mathit{\beta }}_{\mathrm{f}}^{{\mathrm{dT}}_{\mathrm{a}}}\right),\end{array}$$

where β_s, β_w, and β_f are coefficients that correspond to soil, wood, and foliage, respectively; dT_s and dT_a are soil and air temperature deviations from the values at the beginning of the night; and α1, α2, and α3 are partitioning coefficients fixed at 0.5, 0.1, and 0.4, respectively. The β_s, β_w, and β_f coefficients were calculated as follows: β_s values were based on Q₁₀ from the Grünzweig et al. (2009) study at the same site, where β_s = 2.45 for wet soil (i.e., SWC in the upper 30 cm above 20 % vol.); β_s = 1.18 for dry soil (i.e., SWC in the upper 30 cm equal to or below 20 % vol.); β_f=3.15–0.036 T_a; and ${\mathit{\beta }}_{\mathrm{w}}=\mathrm{1.34}+\mathrm{0.46}\exp (-\mathrm{0.5}((\mathrm{DoY}-\mathrm{162})/\mathrm{66.1}{)}^{\mathrm{2}})$, where DoY is the day of the hydrological year starting from 1 October. Finally, GPP was calculated as GPP = NEE–R_e. Negative values of NEE and GPP indicated that the ecosystem was a CO₂ sink.

Half-hourly auxiliary measurements used in this study included photosynthetic activity radiation (PAR; mol m⁻² s⁻¹), vapor pressure deficit (VPD; kPa), wind speed (m s⁻¹), and relative humidity (RH; %), with additional measurements as described elsewhere (Tatarinov et al., 2016). Furthermore, the soil microclimatology half-hourly measurements were measured and calculated with soil chamber measurements, using the LI-8150-203 (LI-COR, Lincoln, NE), as described below, namely air temperature (T_a; ^∘C) and relative humidity at 20 cm above the soil surface and soil temperature (T_s; ^∘C) at a 5 cm soil depth using a soil temperature probe, as well as volumetric soil water content (SWC₀₋₁₀; m³ m⁻³) in the upper 10 cm of the soil near the chambers, using the ThetaProbe model ML2x (Delta-T Devices Ltd., Cambridge, UK), which was calibrated to the soil composition based on the manufacturer's equations.

2.3 Soil CO₂ fluxes

Soil CO₂ fluxes (R_s) were measured with automated non-steady-state systems, using 20 cm diameter opaque chambers and a multiplexer to allow for simultaneous control of several chambers (LI-8150, -8100-101, -8100-104; LI-COR, Lincoln, NE). The precision of CO₂ measurements in the chambers' air is ±1.5 % of the measurements' range (0–20 000 ppm). The chambers were closed on preinstalled PVC collars of 20 cm diameter, allowing for a short measurement time (i.e., 2 min), and positioned away from the collars for the rest of the time. Data were collected using a system in which air from the chambers was circulated (2.5 L min⁻¹) through an infrared gas analyzer (IRGA) to record CO₂ (µmol CO₂ mol⁻¹ air) and H₂O (mmol H₂O mol⁻¹ air) concentrations in the system logger (1 s⁻¹). Gap filling of missing data due to technical problems (i.e., 27 % of the data across the study period between November 2015 and October 2016) was based on the average diurnal cycle of each month.

The rates of soil CO₂ flux, R_s (µmol CO₂ m⁻² s⁻¹), were calculated from chamber data using a linear fit of change in the water-corrected CO₂ mole fraction using Eq. (4) as follows:

$$\begin{array}{}\text{(4)}& R_{\mathrm{s}}={\displaystyle \frac{\mathrm{dC}}{\mathrm{d}t}}\cdot {\displaystyle \frac{vP}{s{T}_{\mathrm{a}}R}},\end{array}$$

where dC∕dt is the rate of change in the water-corrected CO₂ mol fraction (µmol CO₂ mol⁻¹ air s⁻¹), v is the system volume (m³), P is the chamber pressure (Pa), s is the soil surface area within the collar (m²), T_a is the chamber air temperature (K), and R is the gas constant (J mol⁻¹ K⁻¹). A measurement period of 2 min was used, based on preliminary tests to obtain the most linear increase of CO₂ in the chambers with the highest R².

Soil CO₂ fluxes in the experimental plot were measured between November 2015 and October 2016 by means of three measurement chambers using 21 collars grouped in seven sites in the forest stand, with three locations (i.e., three collars) per site, based on different distances from the nearest tree (D_t). The collars were inserted 5 cm into the soil. Data were recorded on a half-hourly basis (48 daily records). The three chambers were rotated between the seven sites every 1–2 weeks to cover all sites and to assess spatial and temporal variations.

Upscaling of the collar measurements to plot-scale soil CO₂ flux was carried out by grouping collars based on three locations (i.e., under trees (< 1 m from nearest tree; UT), in gaps between trees (1–2.3 m; BT), and in open areas (> 2.3 m; OA)), with one chamber taking measurements at each location, and estimating the fractional areas (∅) of the three locations based on mapping the sites according to the distances noted above, as previously done by Raz-Yaseef et al. (2010):

$$\begin{array}{}\text{(5)}& {\displaystyle }{R}_{\mathrm{s}}={R_{\mathrm{s}}}_{\mathrm{OA}}\cdot {\mathit{\varnothing }}_{\mathrm{OA}}+{R_{\mathrm{s}}}_{\mathrm{BT}}\cdot {\mathit{\varnothing }}_{\mathrm{BT}}+{R_{\mathrm{s}}}_{\mathrm{UT}}\cdot {\mathit{\varnothing }}_{\mathrm{UT}},\text{(6)}& {\displaystyle }{\mathit{\varnothing }}_{\mathrm{OA}}+{\mathit{\varnothing }}_{\mathrm{BT}}+{\mathit{\varnothing }}_{\mathrm{UT}}=\mathrm{1}.\end{array}$$

The annual scale of R_s was derived from the upscaled chamber measurements (Eq. 5) based on daily records (48 half-hourly values) of spatially upscaled R_s.

Estimating the temperature sensitivity of R_s (Q₁₀) was performed as described by Davidson and Janssens (2006) using a first-order exponential equation (see also Xu et al., 2015):

$$\begin{array}{}\text{(7)}& R_{\mathrm{s}}=a{e}^{b{T}_{\mathrm{s}}},\end{array}$$

where R_s represents the half-hourly spatially upscaled time series of soil respiration flux (µmol m⁻² s⁻¹), T_s (^∘C) is soil temperature at a 5 cm depth (upscaled spatially and temporally using the same method as for R_s), and a and b are fitted parameters. The b values were used to calculate the Q₁₀ value according to the following equation:

$$\begin{array}{}\text{(8)}& Q_{\mathrm{10}}=e^{\mathrm{10}b}.\end{array}$$

2.4 Soil CO₂ flux partitioning

The determination of different sources of soil CO₂ efflux was based on linear mixing models (Lin et al., 1999) to estimate proportions of three main sources (autotrophic, heterotrophic, and abiotic), using isotopic analysis of soil CO₂ profiles and soil incubation data from eight campaigns (January to September) during 2016, according to Eqs. (9)–(11). Partitioning of the monthly R_s values into components was done using a three-end-member triangular model for interpreting the δ¹³C and Δ¹⁴C values of CO₂ flux; the three-end-member triangular corners are the autotrophic (R_sa), heterotrophic (R_h), and abiotic (R_i) sources of R_s. The δ¹³C and Δ¹⁴C isotope signatures of monthly R_s locate it inside the triangle (Fig. S2):

$$\begin{array}{}\text{(9)}& {\displaystyle }{\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{R_{\mathrm{s}}}=f_{\mathrm{sa}}\cdot {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{sa}}+f_{\mathrm{h}}\cdot {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{h}}+f_{\mathrm{i}}\cdot {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{i}},\text{(10)}& {\displaystyle }{\mathrm{\Delta }}^{\mathrm{14}}{\mathrm{C}}_{R_{\mathrm{s}}}=f_{\mathrm{sa}}\cdot {\mathrm{\Delta }}^{\mathrm{14}}{\mathrm{C}}_{\mathrm{sa}}+f_{\mathrm{h}}\cdot {\mathrm{\Delta }}^{\mathrm{14}}{\mathrm{C}}_{\mathrm{h}}+f_{\mathrm{i}}\cdot {\mathrm{\Delta }}^{\mathrm{14}}{\mathrm{C}}_{\mathrm{i}},\text{(11)}& {\displaystyle }\mathrm{1}=f_{\mathrm{sa}}+f_{\mathrm{h}}+f_{\mathrm{i}},\end{array}$$

where f indicates the fraction of total soil flux (e.g., $f_{\mathrm{h}}=R_{\mathrm{h}}/R_{\mathrm{s}}$), while subscripts sa, h, and i indicate autotrophic, heterotrophic, and inorganic components, respectively. The three-equations system was used to solve the three unknown f fractions of the total soil flux based on empirical estimates of the isotopic end-members. Additionally, δ¹³C and Δ¹⁴C are the stable and radioactive carbon isotopic ratios, where δ¹³C = [([¹³C∕¹²C]_sample∕ [¹³C∕¹²C]_reference) −1] ⋅1000 ‰, and the reference is the Vienna international standard (VPDB). Radiocarbon data are expressed as Δ¹⁴C in parts per thousand or per mil (‰), which is the deviation of a sample ¹⁴C ∕¹²C ratio relative to the OxI standard in 1950 (see Taylor et al., 2015), that is, Δ¹⁴C = [([¹⁴C∕¹²C]_sample∕ (0.95⋅ [¹⁴C∕¹²C]_reference⋅exp [$(y-\mathrm{1950})/\mathrm{8267}$])) −1] ⋅1000 ‰, where y is the year of sample measurements.

The δ¹³C${}_{R_{\mathrm{s}}}$ was estimated monthly using the Keeling plot approach (Figs. S3 and S4; Pataki et al., 2003; Taneva and Gonzalez-Meler, 2011). Soil air was sampled using closed-end stainless-steel tubes (6 mm diameter) perforated near the tube bottom at four depths (30, 60, 90, and 120 cm). Samples of soil air were collected in pre-evacuated 150 mL glass flasks with high-vacuum valves, the dead volume in the tubing and flask necks having been purged with soil air using a plastic syringe equipped with a three-way valve.

Note that the Keeling plot approach is based on the two-end-member mixing model (see Review of Pataki et al., 2003), which often does not hold in soils because of variations in the δ¹³C values of source material with depth (see a recent example in Joseph et al., 2019). However, probably because of the very dry conditions at our study site, no change in δ¹³C with depths in the root zone is observed (±0.1 ‰ across the 35 cm depth profiles; Fig. S5), providing an opportunity to avoid this caveat; we must also conclude of course that the variations among the contributions of R_sa, R_h, and R_i do not change significantly with depth, permitting the use of the single set of isotopic signatures in Table 2. The soil CO₂ samplings carried out therefore represented predominantly the mixing of atmospheric CO₂ with a single integrated soil source signal, consistent with the Keeling plot approach.

The autotrophic (δ¹³C_sa) end-member was estimated based on incubations during the sampling periods of excised roots, following Carbone et al. (2008). Fine roots (< 2 mm diameter) were collected, rinsed with deionized water, and incubated for 3 h in 10 mL glass flasks connected with Swagelok Ultra-Torr tee fittings to 330 mL glass flasks equipped with Louwers high-vacuum valves. The flasks were flushed with CO₂-free air at room temperature close to field conditions. The CO₂ was allowed to accumulate to at least 2000 ppm (∼2 h).

The heterotrophic (δ¹³C_h) end-member was estimated as in Taylor et al. (2015), and, similar to the root-incubation experiment, soil samples from the top 5 cm of the litter layer or 10 cm below the soil surface were collected, and roots were carefully removed to isolate heterotrophic components. Root-free soils were placed in 10 mL glass flasks and allowed to incubate for 24 h before being transferred to evacuated 330 mL glass flasks. The inorganic source (δ¹³C_i) end-member was estimated using 1 g of dry soil (ground to pass through a 0.5 mm mesh) placed in a 10 mL tube with a septum cap; then, 12 mL of 1 M HCl was added to dissolve the carbonate fraction, and the fumigated CO₂ withdrawn from each tube was collected using a 10 mL syringe and injected into a 330 mL evacuated flask for isotopic analysis.

Radiocarbon estimates were based on the work of Carmi et al. (2013) at the same site, adjusted to the measured atmospheric ¹⁴C values during the study period (49.5 ‰; Carmi et al., 2013). The Δ¹⁴C_sa and Δ¹⁴C_h end-members were estimated based on the assumption that they carry the ¹⁴C signatures of 4 and 8.5 years, respectively, older than the ¹⁴C signature of the atmosphere at the time of sampling, based on mean ages previously estimated (Graven et al., 2012; Levin et al., 2010; Taylor et al., 2015). The ratio Δ¹⁴C_i was obtained from Carmi et al. (2013). Monthly values of Δ¹⁴C${}_{R_{\mathrm{s}}}$ were obtained using the linear equation of the regression line of the measured δ¹³C values of R_sa, Rsh, and R_i and the corresponding estimated Δ¹⁴C values (Fig. S2) and monthly δ¹³C values of R_s.

2.5 Isotopic analysis

Isotopic analysis followed the methodology described in Hemming et al. (2005). The δ¹³C of CO₂ in the air was analyzed using a continuous-flow mass spectrometer connected to a 15-flask automatic manifold system. An aliquot of 1.5 mL of air was expanded from each flask into a sampling loop on a 15-position valve (Valco, Houston, TX, USA). CO₂ was cryogenically trapped from the air samples using helium as a carrier gas; it was then separated from N₂O with a Carbosieve G (Sigma-Aldrich) packed column at 70 ^∘C and analyzed on a Europa 20-20 isotope ratio mass spectrometer (IRMS; Sercon, Crewe, UK). The δ¹³C results were quoted in parts per thousand (‰) relative to the VPDB international standard. The analytical precision was 0.1 ‰. To measure [CO₂], an additional 40.0 mL subsample of air from each flask was expanded into mechanical bellows and then passed through an infrared gas analyzer (LI-6262; LI-COR, Lincoln, NE, USA) in an automated system. The precision of these measurements was 0.1 ppm. Flasks filled with calibrated standard air were measured with each batch of 10 sample flasks; five standards were measured per 10 samples for δ¹³C analyses and four standards per 10 samples for [CO₂] analyses.

Organic matter samples were dried at 60 ^∘C and milled using a Wiley mill fitted with a size 40 mesh, and soil samples were ground in a pestle and mortar. Soils containing carbonates were treated with 1 M hydrochloric acid. Between 0.2 and 0.4 mg of each dry sample was weighed into tin capsules (Elemental Microanalysis Ltd., Okehampton, UK), and the δ¹³C of each was determined using an elemental analyzer linked to a Micromass Optima IRMS (Manchester, UK). Three replicates of each sample were analyzed, and two samples of a laboratory working standard cellulose were measured for every 12 samples. Four samples of the acetanilide (Elemental Microanalysis Ltd.) international standard were used to calibrate each run, and a correction was applied to account for the influence of a blank cup. The precision was 0.1 ‰.

2.6 Total belowground carbon allocation (TBCA)

TBCA (g C m⁻² yr⁻¹) was calculated following Giardina and Ryan (2002) for the study year (November 2015–October 2016) as follows:

$$\begin{array}{}\text{(12)}& \mathrm{TBCA}=R_{\mathrm{s}}-R_{\mathrm{alp}}+\mathrm{\Delta }{\mathrm{C}}_{\mathrm{soil}},\end{array}$$

where R_alp is the annual aboveground litter production between November 2014 and October 2015, and ΔC_soil is the annual change in belowground total soil organic C. Litter production, not measured during the present study, was estimated based on values obtained by Masyk et al. (2008) for 2000–2006 (56 g C m⁻² yr⁻¹) and assumed to have increased in the study period (2014–2015) proportionally to the measured increase in leaf area index (LAI; 1.31 to 1.94; i.e., R_alp= [$(\mathrm{1.94}\cdot \mathrm{56})/\mathrm{1.31}$] = 83 g C m⁻² yr⁻¹). For herbaceous litter production, three plots of 25 m² were randomly selected in 2002 and harvested at the end of the growing season, total fresh biomass was weighed, and subsamples were used to determine dry weight and C content. Grünzweig et al. (2007) found that herbaceous litter production was close to the average rainfall for the specific year; this method was adapted in the current study for the period between November 2014 and October 2015. Since aboveground litter (R_alp; the sum of tree litter and herbaceous litter production) of a given year was mainly produced during that year but decayed during the following hydrological year, TBCA was based on the current year's R_s (2015–2016) and the previous year's R_alp (2014–2015). ΔC_soil was set constant as the average annual belowground carbon increase since afforestation (Qubaja et al., 2019).

2.7 Statistical analyses

Two-way ANOVA tests were performed at a significance level set at p=0.001 to detect significant effects of locations (OA, BT, and UT), sites, and their interactions on R_s and metrological parameters. Pearson correlation analysis (r) was used to detect the correlation between R_s and meteorological parameters. To quantify spatiotemporal variability in R_s, the coefficient of variation (CV %) was calculated as [(STDEV ∕ Mean) ⋅ 100 %]. Heterogeneity was considered weak if CV % ≤ 10 %, moderate if 10 % < CV % ≤ 100 %, and strong if CV % > 100 %. All the analyses were performed using MATLAB software, Version R2017b (MathWorks, Inc., MA, USA).

Table 1 Annual mean of half-hourly values across locations (OA, open area; BT, between trees; UT, under tree) in seven sites in the forest during the study period of soil respiration flux rates (R_s) together with the soil water content at 10 cm depth (SWC), minimum distances from nearby tree (D_t), soil temperature at 5 cm depth (T_s), and air temperature (T_a) and relative humidity (RH) at the soil surface (numbers in parentheses indicate ±SE).

^* Correlation is significant at the 0.05 level (two-tailed). ${}^{**}$ Correlation is significant at the 0.01 level (two-tailed).

3 Results

3.1 Spatial variations

The spatial variations in R_s across locations (distance from nearest tree) and sites (across the study area) are reported in Table 1, together with other measured variables. The results indicated an overall mean R_s value of 0.8±0.1 µmol m⁻² s⁻¹, with distinct values for the three locations. R_s was greater at UT locations than at the BT and OA locations by a factor of ∼2. The spatial variability among the locations was also apparent in the R_s daily cycle (Fig. 1), with clear differences between the wet season (November to April), when the UT showed consistently higher R_s values than at other locations by a factor of about 1.6, and the dry season, when the equivalent values differed by a factor of approximately 2.6. Note that the daily peak in R_s remained at midday in both the wet and dry seasons. Overall, the 21 collars showed moderate variations (CV = 55 %; Table 1); R_s was negatively correlated with distance from trees (D_t; $r=-\mathrm{0.62}$; p<0.01) and with soil and air temperatures (T_s and T_a; $r=-\mathrm{0.45}$; p<0.05) and positively correlated with soil water content and relative humidity (SWC and RH; r=0.50; p<0.05). The inverse correlation between R_s and distance from the nearest tree could be useful in considering the expected decline in stand density due to thinning and mortality (e.g., associated with a drying climate). For a first approximation, the results indicate that decreasing from the present stand density of 300 to 100 trees ha⁻¹ and the resulting increase in mean distance among trees could result in decreasing ecosystem R_s by 11 %.

Figure 1 Representative diurnal cycles of soil respiration (R_s; using soil chambers across locations: open area, OA; between trees, BT; under trees, UT) and sites in panels (a) and (b), of net ecosystem exchange (NEE; canopy-scale eddy covariance) and gross primary production (GPP) and ecosystem respiration (R_e) and its partitioning to soil respiration (R_s) and aboveground tree respiration (R_f) in panels (c) and (d), during the wet (November–April) and dry (May–October) periods. Based on half-hourly values over the diurnal cycle; shaded areas indicate ±SE; R_f was estimated as the residual as $R_{\mathrm{f}}=R_{\mathrm{e}}-R_{\mathrm{s}}$ and is presented as a dashed line.

3.2 Temporal dynamics

On the diurnal timescale, CO₂ fluxes showed typical daily cycles (Fig. 1). As expected, on average, all CO₂ fluxes were higher during the wet period compared to the dry season by a factor of ∼2. However, R_s and R_e peaked around midday in both the wet and dry seasons, while the more physiologically controlled NEE and GPP showed a shift from midday (around 11:00–14:00 LT) to early morning (08:00–11:00 LT) in the dry season, with a midday depression and a secondary afternoon peak (Fig. 1d).

Table 2 The δ¹³C and Δ¹⁴C signature of soil respiration (R_s) and its partitioning into autotrophic (R_sa), heterotrophic (R_h), and abiotic (R_i), together with the relative contribution of each to the soil and ecosystem respiration for Yatir Forest during eight campaigns of measurements from January to September 2016 (numbers in parentheses indicate ±SE) in comparison to results obtained previously in the same forest (2001–2006 mean values). The monthly contribution of R_sa, R_h, and R_i to R_s or R_e is presented in Fig. 3a and b, respectively.

¹ Measured in the present study. ² Measured by Carmi et al. (2013). ³ Calculated based on the measured atmospheric value by Carmi et al. (2013). ⁴ Calculated based on the best-fit regression equation in Fig. S2.

The temporal variations across the seasonal cycle are reported in Fig. 2, based on monthly mean values and exhibiting sharp differences between the wet and dry seasons. As previously observed in this semiarid site, all CO₂ fluxes peak in early spring between March and April. The corresponding high-resolution data are reported in Fig. S6, which show also that the high winter (February) R_s rates were associated with clear days when photosynthetic active radiation (PAR) increased with air temperature, T_a. These data also show that, following rainy days, daily R_s values could reach 6.1 µmol m⁻² s⁻¹ (i.e., in the UT microsite; data not shown), although the average was 1.1±0.2 µmol m⁻² s⁻¹ during the wet period, which diminished by ∼55 % in the dry season to mean daily values of 0.5±0.1 µmol m⁻² s⁻¹. In spring (April), all CO₂ fluxes peaked during the crossover trends of decreasing soil moisture content and increasing temperature and PAR (Fig. S6).

Figure 2 Seasonal trends of monthly mean values during the research period of (a) the fluxes of net ecosystem exchange (NEE), gross primary production (GPP), and ecosystem respiration (R_e) and its components, soil respiration (R_s) and aboveground tree respiration (R_f), and monthly mean of precipitation (P), and monthly mean of key environmental parameters; (b) soil water content at the top 10 cm (SWC₀₋₁₀) and soil temperature at 5 cm (T_s); and (c) vapor pressure deficit (VPD) and photosynthetic activity radiation (PAR). R_f is obtained from R_e-R_s. Vertical dotted lines indicate the winter, spring, summer, and fall seasons.

The temporal variations in the half-hourly values of R_s reflected changes in soil moisture at 0–5 cm depth and PAR (r=0.5 and 0.2, respectively; p<0.01) and negative correlations with T_s and RH (r=0.2 and 0.1, respectively; p<0.01). The variations in the integrated R_s showed a CV of 71 %, with the temporal variations dominated strongly by PAR (CV > 100 %), moderately by SWC (CV ∼85 %), and weakly by RH (CV ∼40 %; correlations and CV values were not included in figures and tables). Repeating the models applied by Grünzweig et al. (2009), the potential climatic factors that best predicted daily R_s shifted from SWC and PAR in the dry season to T_s and PAR in the wet season (Table S2). These equations explained 43 % and 70 % of the variation in R_s in the dry and wet seasons, respectively (Table S2). A reasonable forecast of the temporal variations in R_s (µmol m⁻² s⁻¹) at half-hourly values (R²=0.60, p<0.0001) was obtained based on SWC₀₋₁₀ and T_s values across the entire seasonal cycle, based on

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{R}_{\mathrm{s}}& =\mathrm{0.05126}\cdot \exp (\mathrm{0.04274}\cdot T_{\mathrm{s}}\\ & +\mathrm{28.51}\cdot \mathrm{SWC}-\mathrm{74.44}\cdot {\mathrm{SWC}}^{\mathrm{2}}).\end{array}\end{array}$$

At the ecosystem scale, R_e was characterized by high fluxes in the wet season and peak values of ∼2.4 µmol m⁻² s⁻¹ in February to April (Fig. 2; Table S1). R_e fluxes rapidly decreased after the cessation of rain and reached the lowest values in the fall (September to October), with mean dry-period values of 0.5±0.1 µmol m⁻² s⁻¹ (Fig. 2, Table S1). GPP had a mean value of $-\mathrm{1.8}\pm \mathrm{0.4}$ µmol m⁻² s⁻¹, and daily NEE had a mean value of $-\mathrm{0.5}\pm \mathrm{0.3}$ µmol m⁻² s⁻¹ (Table S1 and Fig. S6), with the same seasonality for both (Fig. 2).

Figure 3 (see also Table 2) summarizes the seasonal variations in R_s and R_e partitioning. The monthly R_sa and R_h were not significantly different but were significantly different from R_i (p<0.05). The R_sa∕R_s ratios ranged from 0.32 to 0.46, the largest contribution occurring in early spring from February to April. The R_h∕R_s fraction ranged between 0.33 and 0.45, being the highest during the wet season. The R_i∕R_s fraction – the fraction of inorganic sources from the total soil respiration – ranged from 0.09 to 0.35, the largest contribution being in the driest period. The mean relative contributions of these components to R_s over the sampling campaigns are presented in Fig. 3a, but, on average, soil biotic fluxes were higher than abiotic fluxes by a factor of ∼4. R_e partitioning showed an average increase in R_f∕R_e from 25 % in the wet season to 54 % in the dry season and a decline in R_s∕R_e from 75 % to 46 % on average from the wet to the dry season, respectively, which reflected a seasonal change of R_f in the wet season to peak values in the dry season (Fig. 3b). Both the highest and lowest R_s fractions (∼0.74 and nearly 0.34) along the seasonal cycle were associated with low total R_e fluxes, that is, in the fall before the R_f peak in the spring and in the summer, when physiological controls limited water loss (Fig. 2).

Figure 3 (a) Seasonal variations in the relative contribution of soil autotrophic (R_sa), heterotrophic (R_h), and abiotic (R_i) components to R_s, and (b) seasonal variations in the relative contribution of soil autotrophic (R_sa), heterotrophic (R_h), abiotic (R_i), and foliage and stem respiration (R_f is obtained from R_e−R_s) components to ecosystem respiration (R_e) during eight campaigns (January–September) in 2016. The contributions were estimated with linear mixing models using δ¹³C and Δ¹⁴C of soil respiration (R_s) and a soil CO₂ profile method at 0 to 120 cm soil depth. Vertical dotted lines indicate the winter, spring, summer, and fall seasons. These results confirmed earlier estimates of Grünzweig et al. (2009) and Maseyk et al. (2008a).

Table 3 Mean annual values of ecosystem respiration (R_e), its components and associated ratios, net ecosystem exchange (NEE; from eddy covariance), net primary productivity (NPP), gross primary productivity (GPP), carbon-use efficiency (CUE), leaf area index (LAI), and ratio of total belowground carbon allocation (TBCA) to GPP (TBCA / GPP) in the present study (mean of November 2015 to October 2016) and in comparison to results obtained previously in the same forest (2001–2006 mean values). R_i, R_h, R_sa, R_s, R_l and R_w denote abiotic, heterotrophic, soil autotrophic, soil, foliage, and wood CO₂ flux, respectively. Q₁₀ is derived during the two studies for the wet and dry seasons.

¹ SWC ≥0.2 (m³ m⁻³). ² SWC < 0.2 [m³ m⁻³]. ³ The mean of GPP used by Grünzweig et al. (2009) to estimate the TBCA ∕ GPP ratio was 834 g m⁻² yr⁻¹.

3.3 Annual scale

On an annual timescale, estimates of CO₂ flux components based on EC measurements resulted in annual values of GPP, NPP, R_e, and NEP of 655, 282, 488, and 167 g C m⁻² yr⁻¹, respectively (Tables 3 and S1). On average across the measurement period, R_s was the main CO₂ flux to atmosphere, making up 60±6 % of R_e (295±4 g C m⁻² yr⁻¹; Tables 3 and S1), and R_f was another significant component accounting for 40±6 % of R_e (Fig. 3b), which reflected the low-density (300 trees ha⁻¹) nature of the semiarid forest. As indicated above, R_e partitioning showed a decrease in R_s∕R_e and an increase in R_f∕R_e from winter to summer, which is clearly apparent in Fig. 3b. On an annual scale, during the study period, estimates of R_f, R_sa, R_h, and R_i values were 194±36, 119±21, 115±20, and 61±6 g C m⁻² yr⁻¹, respectively. These rates of respiration fluxes translated at the ecosystem scale to R_e ∕ GPP of ∼75 %, lower than observed in other ecosystems (Table S3) and leading, in turn, to high ecosystem CUE of 0.43.

Using the site records of nearly 20 years, long-term trends in GPP, NPP, R_e, and NEP were examined. Soil respiration and its partitioning could not be similarly monitored continuously, but combining the present results with the 2001–2006 values obtained by Grünzweig et al. (2009) and Maseyk et al. (2008a) provided a basis for estimating the long-term trends in soil respiration. Notably, no clear or significant trend over time was observed in any of the canopy-scale continuously monitored fluxes, but, because of relatively large interannual variations, associated mainly with those in precipitation (see Qubaja et al., 2020), it is likely that the relative contributions of the different fluxes, expressed as ratios in Table 3, provide a more robust perspective of the long-term temporal changes in the ecosystem functioning. The results presented in Table 3 reflect the long-term growth of the forest, with a relatively large increase in LAI, but the TBCA remained around 40 %. The results also indicated little change in the total soil respiration, R_s, component, (as a fraction of R_e or GPP) but a general shift from the autotrophic components to the heterotrophic component (i.e., R_h). This was reflected in the decreasing ratio of the autotrophic components (i.e., R_sa, R_l, and R_w) and the increasing ratio of R_h to R_e (Table 3) across the 13-year observation period (2003 to 2016).

4 Discussion

Partitioning ecosystem carbon fluxes and long-term observational studies are key to understanding ecosystem carbon dynamics and their response to change. Overall, the results support our research hypothesis that the observed high CUE at our site is at least partly due to the characteristics of the carbon flux partitioning that can be associated with the semiarid conditions. Compared to other sites and climates (see comparative compilation in Table S3 in the Supplement), the results reflect several key points: (1) relatively high belowground allocation; (2) low soil respiration in general and low heterotrophic respiration in particular; (3) combining the results for 2016 and those of our earlier study offered a long-term perspective across 13 years, indicating that the low R_s in this ecosystem is robust and exhibits reduced sensitivity to temperature; and (4) there is a general long-term shift from autotrophic to heterotrophic respiration.

Comparing CO₂ fluxes in this forest with fluxes in a range of European forests showed that mean NEP in the semiarid forest (167 g C m⁻² yr⁻¹) is similar to the mean NEP in other European forests (150 g C m⁻² yr⁻¹; FLUXNET).

Carbon partitioning belowground (TBCA / GPP) was relatively high (∼38 %), with little change across the long-term observation period. It is, however, within the range of mean values for forests in different biomes (Litton et al., 2007). High belowground allocation helps explain the high rate of SOC (soil organic carbon) accumulation observed over the period since afforestation (Grünzweig et al., 2007; Qubaja et al., 2019). Note that, irrespective of the soil carbon accumulation, the abiotic component to the CO₂ flux seems to be significant in dry environments (Table 3) and in particular in the dry seasons, when biological activities drastically decrease (Kowalski et al., 2008; Lopez-Ballesteros et al., 2017; Serrano-Ortiz et al., 2010; Martí-Roura et al., 2019). The results show that considering the abiotic effects on estimating soil respiration and, in turn, on estimating the carbon budget in dry calcareous soils can play an important part in estimating soil and ecosystem respiration fluxes (Angert et al., 2015; Roland et al., 2012).

The soil CO₂ efflux in the semiarid forest (295 g C m⁻² yr⁻¹) is at the low end of R_s values across the range of climatic regions, from 50 to 2750 g C m⁻² yr⁻¹ (Adachi et al., 2017; Chen et al., 2014; Grünzweig et al., 2009; Hashimoto et al., 2015). This is clearly lower than the mean R_s value for global evergreen needle forests, which is estimated at 690 g C m⁻² yr⁻¹ (Chen et al., 2014), and between estimates for desert scrub and Mediterranean woodland (224–713 g C m⁻² yr⁻¹; Raich and Schlesinger, 1992) or for Mediterranean forests (561–1,015 g C m⁻² yr⁻¹; Casals et al., 2011; Luyssaert et al., 2007; Matteucci et al., 2015; Misson et al., 2010; Rey et al., 2002; Rodeghiero and Cescatti, 2005). The mean instantaneous rate of R_s, 0.8 µmol m⁻² s⁻¹, is also in the range reported for unmanaged forest and grassland in the dry Mediterranean region (0.5 and 2.1 µmol m⁻² s⁻¹; Correia et al., 2012).

The observed low R_s values were associated with a relatively high fraction of autotrophic and a lower fraction of heterotrophic respiration. The mean annual-scale R_sa∕R_s ratio of 0.40 was at the high end of the global range of 0.09 to 0.49 (Chen et al., 2014; Hashimoto et al., 2015). In contrast, heterotrophic respiration showed an annual-scale R_h∕R_s ratio of 0.39±0.02 (Table 2 and Fig. 3), which is lower than the estimated global mean R_h∕R_s value of 0.56 (Hashimoto et al., 2015) but within the range of Mediterranean region forest, which varies between 0.29 and 0.77 (Casals et al., 2011; Luyssaert et al., 2007; Matteucci et al., 2015; Misson et al., 2010; Rey et al., 2002; Rodeghiero and Cescatti, 2005). The relatively low annual scale of the heterotrophic respiration to R_s is consistent with the dry soil over much of the year in this forest (Figs. 2 and S6) and the observed low decomposability of plant detritus and the high mean SOC accumulation rate (Grünzweig et al., 2007).

The long-term perspective from the 13-year observation period indicates emerging trends that can be a basis for assessing the effects of forest age and the evident increase in LAI (Table 3) and changes in environmental conditions (generally warming and drying; see, e.g., Lelieveld et al., 2012). Here, because comparing the noncontinuous data from the present (2016) and earlier (2001–2006) studies is sensitive to the large interannual variations in the ecosystem activities and fluxes (Qubaja et al., 2019), we focused on the more robust changes in the ratio of the respiration components to the overall fluxes (R_e; Table 3). This shows a shifting trend from the autotrophic components to the heterotrophic, with little change in the contribution of R_s to the overall efflux. The ratios of R_sa, R_l, and R_w to R_e tended to decrease by about 13 %, while that of R_h increased by about 18 %; similar trends were seen in soil respiration, with R_sa∕R_s decreasing by −19 % and R_h∕R_s increasing by +8 % (Table 3). The relatively low R_s under conditions of high temperature in the semiarid ecosystem implies reduced sensitivity of respiration to temperature. This is partly imposed by low SWC conditions during extended parts of the year (Grünzweig et al., 2009; cf. Rey et al., 2002; Xu and Qi, 2001). Accordingly, R_s showed greater sensitivity to T_s in the wet period, but, during the 8–9 months of the year when SWC was below ∼0.2 m³ m⁻³, R_s varied predominantly with water availability. The long-term perspective reported in Table 3 indicates a further decrease in temperature sensitivity, with mean Q₁₀ values in the dry season decreasing from 1.6 to 1.1. These estimated Q₁₀ values are generally consistent with published values for different ecosystems (1.4 to 2.0; Hashimoto et al., 2015; Zhou et al., 2009) and with low values under low SWC (Reichstein et al., 2003; Tang et al., 2005). This is also consistent with soil warming experiments by 0.76 ^∘C in Mediterranean ecosystems, which decreased the R_s by 16 % and Q₁₀ by 14 % (Wang et al., 2014). Note also that the low temperature sensitivity in the dry season is likely to be related to reduced microbial activity but may also involve downregulation of plant activity (Maseyk et al., 2008a) and drought-induced dormancy of shallow roots (Schiller, 2000). Finally, we also note that the greater importance of moisture availability in influencing respiration is clearly apparent from the observed relationships of R_s and R_h to mean annual precipitation (MAP) in European evergreen needle forests (Fig. S8; see also Grünzweig et al., 2007), which are not observed with respect to mean annual temperature.