Hysteresis of tropical forests in the 21st century

Study area and period

Our study area is the tropics between 15°N–35°S⁶. We divided the study area into three continents and studied them separately: South America, Africa, and Australasia. Australasia includes Australia and southeast Asia, but excludes southern India. Our results are generated on 0.25° spatial resolution. We classify a cell as forest if it contained at least 50% tree cover (‘forest cover’ in this manuscript) in 1999 according to the dataset from ref. ⁴¹. The moisture recycling simulations were carried out for 2003–2014 (‘recent climate’), for which a consistent set of input data was available (see also ref. ¹¹). ‘Late 21st century’ refers to 2071–2100.

Local-scale forest hysteresis

Previous research has shown that tropical forests may have local-scale tipping points at certain mean annual rainfall levels, but are also affected by the seasonality of that rainfall^5,6,8. Local-scale tipping points for forest were determined using tree cover data following a method from ref. ⁶. Using potential analysis⁴², an empirical stability landscape (as in Fig. 1a) is constructed based on spatial distributions of tree cover against environmental variables such as mean annual rainfall for each continent separately. For each value of the environmental variable, the probability density of tree cover was determined using the MATLAB function ksdensity with a bandwidth of 5%. We applied Gaussian weights to the environmental variable with a standard deviation of 0.05 times the length of the axis of the environmental variable. Local maxima of the resulting probability density function are interpreted as stable states, where we ignored local maxima below a threshold value of 0.004. We used Landsat tree cover data for 2000 on 30 m resolution downloaded for every 0.01°⁴³. We masked out human-used areas, water bodies, and bare ground using the ESA GlobCover land cover dataset for 2009 on 300 m resolution (values 11–30 and ≥190). From the resulting dataset we randomly sampled one million locations for each continent and used them to construct the stability landscapes⁶ against mean annual rainfall and average MCWD. MCWD is the cumulative difference between evapotranspiration and rainfall using monthly averages of those fluxes calculated for each calendar year⁴⁴. It is set to zero when monthly rainfall exceeds monthly evapotranspiration and becomes more negative with an increasing water deficit. Following ref. ¹¹, for both mean annual rainfall and MCWD, we took monthly data from the GLDAS 2.0 dataset⁴⁵ for 1970–1999 so the 30-year period leading up to the land-cover sample (for the year 2000) was used.

Forest evapotranspiration

To estimate the fraction of evapotranspiration attributable to forest cover we used the large-scale hydrological model PCR-GLOBWB, run at 0.5° resolution⁴⁶. Per grid cell, the model simulates evapotranspiration for a range of land-cover types. Here, we are specifically interested in the evapotranspiration of forests, or ‘tall natural vegetation’ in PCR-GLOBWB. Note that we here account for both forest transpiration and canopy interception evaporation instead of, as in ref. ¹¹, only transpiration.

PCR-GLOBWB computes the water balance in two soil layers and a groundwater layer. Soil type, fractional area of saturated soil, and the spatiotemporal distribution of groundwater depth are accounted for (see refs. ^46,47). It includes six land-cover types, with spatially varying parameters⁴⁶: tall and short natural vegetation, pasture, rainfed crops, and paddy and non-paddy irrigated crops. The model was forced with WATCH Forcing Data ERA-Interim precipitation, temperature, and reference potential evapotranspiration for 1979–2014⁴⁸. We used monthly evapotranspiration output of PCR-GLOBWB, implying that we assume that forest component of evapotranspiration remains equal within each month. For detailed model descriptions and validation, we refer to earlier studies^11,46,49,50.

Atmospheric moisture tracking

As an essential step in estimating the forest–rainfall feedback, we determined where the moisture from enhanced evapotranspiration precipitates again by using atmospheric moisture tracking. The method for atmospheric moisture tracking resembles that in ref. ¹¹. Apart from the expansion of the study area to the entire tropics, a notable difference is that we here used ERA5 reanalysis data rather than ERA-Interim, meaning that the simulations were based on finer resolution input data (i.e. on 0.25° instead of 0.75° for spatial resolution, and 1 h instead of 3 h for temporal resolution). ERA5 has better performance than ERA-Interim regarding wind fields and rainfall, especially in the tropics^51,52,53. Below we summarize the method (see also ref. ¹¹).

We used a Lagrangian method of moisture tracking that is based on previous studies^11,54,55,56 that track parcels of evaporated moisture forward through the atmosphere to their subsequent precipitation location. Moisture particles that enter the atmosphere are assigned a random location within the 0.25° grid cell and random starting height in the atmosphere scaled with the humidity profile, and their trajectories are then tracked through the atmosphere. The trajectories are forced with the three-dimensional ERA5 reanalysis estimates of wind speed and direction, which were linearly interpolated at every time step of 0.25 h. Water particles in the atmosphere have an equal probability of raining out, regardless of vertical position. Rainfall A (mm per time step) at location x,y and time t that has evaporated from any location of release in any cell is ref. ⁵⁶

where P is rainfall in mm per time step, W_parcel is the water in the tracked parcel in mm, E_source is its fraction of water that evapotranspired from the source, and TPW is the precipitable water in the atmospheric water column in mm. Every time step, the amount of water in the parcel is updated based on evapotranspiration ET into the parcel and rainfall P from it:

The fraction of water in the parcel that has evapotranspired from the source then becomes

Thus, the amount of water that was tracked from the source location decreases with precipitation along its trajectory. Parcels were followed until either less than 5% of its original amount was left in the atmosphere, or the tracking time was 30 days. Any moisture remaining in the parcel when the trajectories end is assumed to rain out over non-land areas, thus not contributing to our analysis. We analysed each continent separately for reasons of computability. However, small moisture flows between forests in different continents can be expected, as has been simulated for flows from Africa to the Amazon⁵⁷. Over all land points, ERA5 hourly evapotranspiration is linearly interpolated to every 0.25-h time step. The moisture flow m_ij in mm per month linking evapotranspiration in cell i to rainfall in cell j where (left[ {x,y} right] , {it{epsilon }} , j) over the course of a given month becomes

where ET_i,t is the evapotranspiration in mm per time step, and W_i,t is the tracked amount of water from source cell i at time step t.

At continental scales, evaporated moisture can rain down and re-evaporate multiple times. This also means that forest evapotranspiration can enhance rainfall multiple times. We accounted for this ‘cascading moisture recycling’ following refs. ^11,14, in which the rainfall attributed to an upwind forest source is subsequently tracked upon re-evaporation. After six re-evapotranspiration cycles, cascading moisture recycling has decreased to practically zero¹¹. Therefore, following ref. ¹¹, seven iterations of cascading moisture recycling were performed.

Hysteresis experiments

We determined the hysteresis of tropical forests through a series of iterative runs; each one started either from a fully forested continent or from a fully deforested continent. We simulated the hypothetical mean annual rainfall levels across the (tropical part of the) continent given this initial condition, that is, rainfall without any forest evapotranspiration or rainfall including evapotranspiration from an entirely forested continent. Next, based on the empirical bifurcation diagrams (i.e. nonforest, bistable forests, and stable forests in each continent are determined based on the bistability ranges shown in Supplementary Figs. 1−3), we determined either the minimal distribution of tropical forest (in case of a no-forest initial condition, based on the higher end of the bistability range) or the maximal distribution (in case of a fully forested initial condition, based on the lower end of the bistability range) at these rainfall levels. Thus, in the simulations with an empty initial condition, only stable forests (green in Fig. 1) would regrow; in those with a full initial condition, both stable and bistable forests (green and yellow in Fig. 1) would disappear. Because the resulting new distribution of forest would generate different levels of rainfall, the procedure was repeated with the respective forest distribution as initial condition. This occurred until rainfall levels had (practically) stabilized between iterative runs (up to three iterations).

We assumed that no other vegetation type replaces the forest in order to show the theoretically possible distributions of tropical forests. This may lead to an overestimation of the actual effects of forest on rainfall, especially if forests would be replaced by highly transpiring crops⁵⁸. Furthermore, land-cover changes will alter wind patterns and therefore the expected coupling between forests through evapotranspiration and rainfall⁵⁹. Fossil fuel emissions not only alter the climate, but the emitted CO₂ also fertilizes plants. This increases trees’ water-use efficiency, reducing their water demand, but it also increases biomass production⁶⁰. The effects of this CO₂ fertilization on the water cycle might be small⁶¹, but its net effects on tropical forest hysteresis remains uncertain.

For display of Fig. 2, the resolution of rainfall values was increased by a factor of 2, to 0.125° and smoothed using a two-dimensional Gaussian kernel with a standard deviation of 0.5°.

Climate-change scenario

As the estimate of late 21st-century rainfall conditions, we used the rainfall output from the SSP5-8.5 scenario simulations by seven available CMIP6 models⁶². These models are BCC-CSM2-MR, CanESM5, CNRM-CM6-1, CNRM-ESM2, IPSL-CM6A-LR, MRI-ESM2.0, and UKESM1.0-LL. We took the mean across the model outputs for the annual rainfall values for 2071–2100. The scenario is considered a severe climate-change scenario. Because the moisture tracking model is forced with atmospheric reanalysis data, we assumed that (forest-induced) moisture flows in the scenario are the same as in the period of our atmospheric simulations (2003–2014).

We assumed that a tipping point from a forested to a nonforested state occurs when mean annual rainfall in a forested cell (forest cover ≥ 50%) is currently (2003–2014) above the lower tipping point (Supplementary Figs. 1–3), but is reduced to below the lower tipping point in the climate-change scenario. Similarly, a tipping point from a nonforested to a forested state occurs when mean annual rainfall in a nonforested cell (forest cover < 50%) is currently (2003–2014) below the upper tipping point (Supplementary Figs. 1–3), but is increased to above the upper tipping point in the climate-change scenario.

To explore whether the CMIP6 models project a change in moisture transport for the late 21st century, we compared the vertically integrated eastward and northward moisture fluxes (in kg m⁻¹ s⁻¹) for 35°S–35°N for 2015–2020, which is the start of the simulation runs, and 2095−2100, the end of the runs. We did this for the same seven models and SSP as mentioned above.

Validation and sensitivity analyses

We conducted a number of additional analyses regarding model validation and uncertainties. We compared our evapotranspiration product GLDAS to estimates from FLUXCOM. Instead of using climate forcing data, FLUXCOM merges energy flux measurements from FLUXNET eddy covariance towers with remote sensing⁶³. Thus, it provides an independent as possible comparison with GLDAS. Over all, the two products agree well with a concordance correlation r² = 0.69 across the tropics (Supplementary Fig. 15). This correspondence is lower when we consider forested areas only (r² = 0.26; Supplementary Fig. 16). Especially at relatively low values of monthly evapotranspiration they differ, where FLUXCOM tends to produce higher estimates of (forest) evapotranspiration than GLDAS. Positive and negative differences exist throughout the tropics, but especially in Africa, FLUXCOM estimates higher evaporation levels than GLDAS (Supplementary Fig. 17). Underestimations of evaporation by GLDAS would imply that changes in forest cover may have larger effects than we currently account for, but systematic bias in flux measurement data might also be responsible⁶³.

We assess the sensitivity of forest hysteresis on each continent to a number of variables. For these sensitivity analyses we performed our atmospheric simulations for 2003 only. We did this for: (1) the share that forest cover contributes to evapotranspiration, using 80, 90, 100, 110, and 120% of the estimated levels used in the main analyses. (2) The values of the bifurcation points, where we simultaneously changed both the lower and upper bifurcation point by −200, −100, 0, 100, and 200 mm per year. (3) The mixing strength of atmospheric moisture along the vertical moisture column. This was shown to be the most important source of uncertainty in Lagrangian atmospheric moisture tracking¹⁶. Here, we applied three levels of atmospheric mixing: low, in which moisture gets assigned a new random vertical location every 120 h; medium, used in the main analyses, in which moisture gets assigned a new location every 24 h; and high, where mixing occurs every hour. These specific analyses were done on 0.5° instead of 0.25°. (4) The CMIP6 climate model, where we estimated the hysteresis for each of the used models separately.

All data analyses were carried out in MATLAB R2019a. Figure 2 was made using Matplotlib 2.2.5.

Reporting summary

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

Source: Ecology - nature.com