Large variations in afforestation-related climate cooling and warming effects across short distances

Study site description

The study was carried out at the edge of an arid region in mature plantations dominated by P. halepensis, of age of 40-50 years (Pinus halepensis), and their adjacent non-forested ecosystems. The sites were distributed across a climatic gradient from arid and semi-arid to dry sub-humid (Fig. 1, Supplementary Table S1). Three selected paired sites included: (1) An arid site at Yatir forest (annual precipitation, P: 280 mm; aridity index, AI: 0.18; elevation: 650 m; light brown Rendzina soil, and forest density: 300 trees ha⁻¹), where a permanent flux tower has been operating since the year 2000 (http://fluxnet.ornl.gov). Note that an AI of 0.2 marks the boundary between arid and semi-arid regions. Yatir, with AI = 0.18, is formally within an arid zone, but on the edge of a semi-arid zone. (2) An intermediate semi-arid site in Eshtaol forest (P: 480 mm; AI: 0.37; elevation: 350 m; light brown Rendzina soil, and forest density: 450 trees ha⁻¹). (3) A dry-subhumid site in northern Israel at the Birya forest (P: 770 mm; AI: 0.64; elevation: 755 m; dark brown Terra-Rossa and Rendzina soil, and forest density: 600 trees ha⁻¹). Non-forest ecosystems were sparse dwarf shrublands, dominated by Sarcopoterium spinosum in a patchy distribution with a wide variety of herbaceous species, mostly annuals, that grew in between the shrubs during winter to early spring, and then dried out. All non-forested sites had been subjected to livestock grazing (exposing soils). Finally, an additional site that was characterized as Oak-forest vegetation was added. The site was dominated by two oak species, Quercus calliprinos and Quercus ithaburensis (P: 540 mm; AI: 0.4; see previous publication for more details on the oak site⁵³). All sites were under high solar radiation load, with an annual average of approximately 240 Wm⁻² in the arid region and only 3% lower in the northern site in the dry-subhumid region (Table 1).

Mobile laboratory

Measurements were conducted on a campaign basis using a mobile lab with a flux measurement system at all sites except the Yatir forest, where the permanent flux tower was used (http://fluxnet.ornl.gov;⁵⁴). Repeated campaigns of approximately two weeks at each site, along the seasonal cycle, were undertaken during 4 years of measurements, 2012–2015 (a total of 6-7 campaigns per site, evenly distributed between the seasons) the 4 years of measurements were found to be representative of previous 70 years of precipitation record (Supplementary Figure S1 and Table S2). The mobile lab was housed on a 12-ton 4 × 4 truck with a pneumatic mast with an eddy-covariance system and provided the facility for any auxiliary and related measurements. Non-radiative flux measurements were undertaken using an eddy-covariance system to quantify CO₂, sensible heat (H), and latent heat (LE) fluxes using a 3D sonic anemometer (R3, Gill Instruments, Hampshire, UK) and an enclosed-path CO₂/H₂O infrared gas analyzer (IRGA; LI-7200, LI-COR). Non-radiative flux measurements were accompanied by meteorological sensors, including air temperature (Ta), relative humidity (RH), and pressure (Campbell Scientific Inc., Logan, UT, USA), radiation fluxes of net solar- and net long-wave radiation (SWnet and LWnet, respectively), and photosynthetic radiation sensors (Kipp & Zonen, Delft, Holland). Raw EC data and the data from the meteorological sensors were collected using a computer and a CR3000 logger (Campbell Sci., Logan, UT, USA), respectively. The EC system was positioned at the center of each field site with the location and height aimed at providing sufficient ‘fetch’ of relatively homogeneous terrains. For detailed information on the use of the mobile lab and the following data processing of short and long-term fluxes see previous publications^29,55,56.

Data processing

Mean 30-min fluxes (CO₂, LE, and H) were computed using Eddy-pro 5.1.1 software (LiCor, Lincoln, Nebraska, USA). Quality control of the data included a spike removal procedure. A linear fit was used for filling short gaps (below three hours) of missing values due to technical failure. Information about background meteorological parameters, including P, Ta, RH, and global radiation (Rg), was collected from meteorological stations (standard met stations maintained by the Israel Meteorological Service, https://ims.gov.il/en). The data were obtained at half-hourly time resolution and for a continuous period of 15 years since 2000.

Estimating continuous fluxes using the flux meteorological algorithm

Estimation of the flux-based annual carbon and radiation budgets was undertaken using the short campaign measurements as a basis to produce a continuous, seasonal, annual, and inter-annual scale dataset of ecosystem fluxes (flux meteorological algorithm). The flux meteorological algorithm method was undertaken based on the relationships between measured fluxes (CO₂, LE, H, SWnet, and LWnet) and meteorological parameters (Ta, RH, Rg, VPD, and transpiration deficit, a parameter that correlated well with soil moisture, see main text and supplementary material of previous publication²⁹. A two-step multiple stepwise regression was established, first between the measured fluxes (H, LE, and the ecosystem net carbon exchange) and the meteorological parameters measured by the mobile lab devices, and then between the two meteorological datasets (i.e., the variables measured by the Israel meteorological stations) for the same measurement times. Annual fluxes were calculated for the combined dataset of all campaigns at each site using the following generic linear equation:

where, y is the ecosystem flux of interest, the daily average for radiative fluxes (LWnet and SWnet), non-radiative fluxes (H and LE), and daily sum for net ecosystem exchange (NEE), a and ({b}_{i}) are parameters, and ({x}_{i}) is Ta, RH, Rg, vapor pressure deficit, or transpiration deficit. The meteorological variables (({x}_{i})) were selected by stepwise regression, with ({b}_{i}=0) when a specific ({x}_{i}) was excluded.

Based on this methodology, ecosystem flux data were extrapolated to the previous 7–15 years (since 2000 in the dry-subhumid and arid sites, since 2004 in the semi-arid sites, and since 2008 in the Oak-forest site) using all the available continuous meteorological parameters from the meteorological stations associated with our field sites. The long-term annual sums and means of extrapolated ecosystem fluxes were averaged for multi-year means of each site for the period of available extrapolated data. In the first, previously published phase of this study²⁹, the extrapolation method was extensively tested, including simulation experiments at the arid forest site, where continuous flux data from the 20 years old permanent flux tower were available. Five percent of the daily data were selected by bootstrap (with 20 repetitions), a stepwise regression was performed for this sample, and then, the prediction of fluxes using Eq. 1 above was performed for the entire observation period. In the current phase of the study additional flux measurements are included with the R² coefficients for the additional measurements ranging between 0.4-0.9, see Supplementary Table S3.

The aridity index of the Oak-forest was in between those of the semi-arid and dry-subhumid Pine-forests (0.4 compared to 0.37 and 0.67, respectively). Therefore, to compare the Oak-forest with Pine-forest and non-forest sites, the average results from the semi-arid and dry-subhumid paired sites were used.

Radiative forcing and carbon equivalence equations

To compare the changes in the carbon and radiation budgets caused by forestation, we adopted the approach of Myhre et al. ³⁰, and used Eq. 2:

where land-use changes in radiative forcing (RF_ΔC) are calculated based on the CO₂ reference concentration, C₀ (400 ppm for the measured period of study), and ΔC, which is the change in atmospheric CO₂ in ppm, with a constant radiative efficiency (RE) value of 5.35. Here, ΔC is calculated based on the annual net ecosystem productivity (NEP; positive carbon gain by the forest, which is identical to net ecosystem exchange (NEE), the negative carbon removal from the atmosphere) as the difference between forested and non-forested ecosystems (ΔNEP) multiplied by a unit conversion constant:

where, k is a unit conversion factor, from ppm to g C (k = 2.13 × 10⁹), calculated as the ratio between the air molar mass (M_a = 28.95; g mol⁻¹), to carbon molar mass (M_c = 12.0107; g mol⁻¹), and total air mass (m_a = 5.15 × 10⁹; g).

Etminan et al. ⁵⁷ introduced an updated approach to calculate the RE as a co-dependent of the change in CO₂ concentration and atmospheric N₂O:

where, (triangle {{{{{rm{C}}}}}}) is the change in atmospheric CO₂ in ppm resulting from the forestation, as calculated in Eq. 3, (bar{{{{{{rm{N}}}}}}}) is the atmospheric N₂O concentration in ppb (323), and the coefficients a₁, b₁, and c₁ are −2.4 × 10⁻⁷ Wm⁻²ppm⁻¹, 7.2 × 10⁻⁴ Wm⁻² ppm⁻¹, and −2.1 × 10⁻⁴ Wm⁻²ppb⁻¹, respectively.

Combining Eqs. 2 and 4 with an airborne fraction of (zeta =0.44)⁵⁸, we obtain Eq. 5:

Next, the annual average radiative forcing due to differences in radiation flux was calculated as follows:

where, ΔR is the difference between forest and non-forest reflected short-wave or emitted long-wave radiation (ΔSWnet and ΔLWnet, respectively), assuming that the atmospheric incoming solar and thermal radiation fluxes are identical for the two, normalized by the ratio of the forest area (({A}_{F})) to the Earth’s area (({A}_{E}=5.1times {10}^{14},{m}^{2})).

As forest conversion mostly has a lower albedo, the number of years needed to balance (‘Break even time’) the warming effect of changes in radiation budget by the cooling effect of carbon sequestration is calculated by combining Eqs. 5 and 6:

The multiyear averages of NEP for each of the three paired sites (forest and non-forest) were then modeled over a forest life span of 80 years. This was done based on a logarithmic model, modified for dryland, which takes as an input the long-term averages of NEP ((overline{{NEP}})) as in Eq. 3:

where annual carbon gain at time t (NEP_t) is a function of the multiyear average carbon gain ((overline{{{{{{rm{NEP}}}}}}})), forest age (t), and growth rate (b). Parameter b is a constant (b = −0.17) and is calculated based on the global analysis of Besnard et al. ⁵⁹, limiting the data to only dryland flux sites⁶⁰. Note that this analysis indicates NEP reaching a steady state and was used here to describe the initial forest growth phase, while growth analyses indicate that carbon sequestration peaks after about 80 years, followed by a steep decline⁵⁰. This is consistent with the time scale for forest carbon sequestration considered here.

In contrast to the one-year differences presented in Table 1 (ΔNEP), the net sequestration potential ((triangle)SP) for each of the paired sites was calculated as the accumulated ecosystem ΔNEPt along with forest age ((triangle) is the difference between forest and non-forest sites):

The ΔSP growth model was compared with previously published data of long-term carbon stock changes in arid forests (i.e., cumulative NEP over 50 years since forest establishment, t = 50), demonstrating agreement within ± 10%²⁷.

For comparison with previous studies, the carbon emission equivalent of shortwave forcing (EESF) was calculated using an inverse version of Eqs. 5 and 6 based on Betts¹:

where, C₀ is the reference atmospheric CO₂ concentration (400 ppm, the average atmospheric concentration for the past decade), RF_ΔR is the multiyear average change in radiative forcing as a result of the change in surface albedo (Eq. 6 W m⁻²), RE is the radiative efficiency (Eq. 4, W m⁻²), ζ is the airborne fraction (0.44 as in Eq. 5), and k is a conversion factor, from ppm to g C (2.13 × 10⁹ as in Eq. 3). Equation 10 was also used to calculate the emission equivalent of longwave forcing (EELF) with the RF_ΔR as the multiyear average change in radiative forcing as a result of the change in net long-wave radiation (ΔLWnet).

Finally, the net equivalent change in carbon stock due to both the cooling effect of carbon sequestration and the warming effect due to albedo change (net equivalent stock change; NESC), was calculated by a simple subtraction:

A comparison of the ΔSP (Eq. 9), the EESF (Eq. 10), and NESC (Eq. 11) with the same metrics as those used in other studies^1,14,61 was done when appropriate. An exception was made for Arora & Montenegro (2011), where only carbon stock changes (ΔSP) were available in carbon units, and biogeophysical (BGP) and biogeochemical (BGC) effects were expressed as temperature changes. To overcome this metric difference, we converted the biogeophysical to carbon equivalent units (EESF + EELF) by multiplying the carbon stock changes (ΔSP) by the ratio between the BGP and BGC effects on temperature (EESF + EELF = ΔSP × BGP/BGC).

Statistical and data analyses

The paired t-test was used to compare multi-annual averages of all variables between forested and adjacent non-forested sites and between sites across the climatic gradient. The variables of interest that were detected for their significant differences were albedo, net radiation and its longwave and shortwave components, latent heat fluxes, sensible heat fluxes, and net ecosystem productivity. All statistical and data analyses were performed using R 3.6.0 (R Core Team, 2020)⁶².

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Source: Ecology - nature.com