Rising ecosystem water demand exacerbates the lengthening of tropical dry seasons

Climate and land cover data

Our study of tropical dry season dynamics required climatic variables with high temporal resolution (i.e., daily) and full coverage of tropic regions. To reduce uncertainties associated with the choice of precipitation (P) and evapotranspiration (Ep or E) datasets, we used an ensemble of eight precipitation products, three reanalysis-based products for Ep, and one satellite-based land E product. These precipitation datasets were derived four gauge-based or satellite observation (CHIRPS⁵⁸, GPCC⁵⁹, CPC-U⁶⁰ and PERSIANN-CDR⁶¹), three reanalyses (ERA-5⁶², MERRA-2⁶³, and PGF⁶⁴) and a multi-source weighted ensemble product (MSWEP v2.8⁶⁵). The potential evapotranspiration (Ep) was calculated using the FAO Penman–Monteith equation⁶⁶ (Eqs. (1, 2)), which requires meteorological inputs of wind speed, net radiation, air temperature, specific humidity, and surface pressure. We derived these meteorological variables from the three reanalysis products (ERA-5, MERRA-2, and GLDAS-2.0⁶⁷). Since PGF reanalysis lacked upward short- and long-wave radiation output and thus net radiation, we used available meteorological outputs from GLDAS-2.0 instead, which was forced entirely with the PGF input data.

Where Ep is the potential evapotranspiration (mm day⁻¹). R_n is net radiation at the surface (MJ m⁻² day⁻¹), T is mean daily air temperature at 2 m height (°C), ({u}_{2}) is wind speed at 2 m height (m s⁻¹), ((,{e}_{s}-{e}_{a})) is the vapor pressure deficit of the air (kPa), ({RH}) is the relative air humidity near surface (%), ∆ is the slope of the saturation vapor pressure-temperature relationship (kPa °C⁻¹), γ is the psychrometric constant (kPa °C⁻¹), G is the soil heat flux (MJ m⁻² day⁻¹, is often ignored for daily time steps G ≈ 0).

We derived the daily evapotranspiration data from the Global Land Evaporation Amsterdam Model (GLEAM v3.3a⁶⁸), which is a set of algorithms dedicated to developing terrestrial evaporation and root-zone soil moisture data. GLEAM fully assimilated the satellite-based soil moisture estimates from ESA CCI, microwave L-band vegetation optical depth (VOD), reanalysis-based temperature and radiation, and multi-source precipitation forcings. The direct assimilation of observed soil moisture allowed us to detect true soil moisture dynamic and its impacts on evapotranspiration. Besides, the incorporation of VOD, which is closely linked to vegetation water content^69,70, allowed us to detect the effect of water stress, heat stress, and vegetation phenological constraints on evaporation. Other observation-driven ET products from remote-sensing physical estimation and flux-tower are not included due to their low temporal resolution (i.e., monthly)⁷¹ or short duration^72,73. ET outputs of reanalysis products are not considered in our analysis, because the assimilation systems lack explicit representation of inter-annual variability of vegetation activities and thus may not fully capture hydrological response to vegetation changes^62,63,67.

We used land cover maps for the year 2001 from the Moderate-Resolution Imaging Spectroradiometer (MODIS, MCD12C1 C5⁷⁴) based on the IGBP classification scheme to exclude water-dominated and sparely-vegetated pixels (like Sahara, Arabian Peninsula). All climate and land cover datasets mentioned above were remapped to a common 0.25° × 0.25° grid and unified to daily resolution. The main characteristics of the datasets mentioned above are summarized in Supplementary Table 1.

Outputs of CMIP6 simulations

To understand how modeled dry season changes compare with observed changes, we analyzed outputs from the “historical” (1983-2014) runs of 34 coupled models participating in the 6th Coupled Model Inter-comparison Project⁷⁵ (CMIP6, Supplementary Table 3). We used these models because they offered daily outputs of all climatic variables needed for our analysis, including precipitation, latent heat (convert to E), and multiple meteorological variables for Ep (air temperature, surface specific humidity, wind speed, and net radiation). All outputs were remapped to a common 1.0° × 1.0° grid and unified to daily resolution.

Defining dry season length and timing

For each grid cell and each dry season definition (P < Ep, P < E and P < (bar{P})), we conducted a harmonic analysis to define the number of dry and wet seasons experienced per year, through a Fourier transform of the entire daily time series^15,16,76. We calculated the ratio of harmonic amplitudes at frequencies of one and two cycles per year to determine seasonality (Fig. 1a, Supplementary Fig. 2). A ratio greater than 0.75 indicates that the harmonic of two cycles per year (i.e., two dry and wet seasons, Fig. 1c) may fit the time series better, otherwise (ratio < 0.75) there is more likely a single dry/wet season (Fig. 1b).

The three definitions of “dry season” that we assessed were: (i) the period when daily precipitation (P) is persistently less than daily potential evapotranspiration (Ep), i.e., P < Ep, (ii) the period when daily precipitation is persistently less than daily actual evaporation (E) i.e., P < Ep, and (iii) the period when daily precipitation is persistently less than the multi-year average daily mean precipitation (left(bar{P}right.)), i.e., P < (bar{P}). Other definitions of “dry season” (e.g., based on a specific rainfall threshold^77,78) have been used in previous research. We chose these three definitions because they can be applied across the entire tropical land area (i.e., they are not locally determined by metrics such as a specific local rainfall threshold value). The dry season should be continuous, not the total number of intermittent dry days. We ensured the continuity for the definition with two rules: (1) regions with bimodal rainfall regime were identified through previous harmonic analysis and discussed separately, for which each single dry season should be continuous, (2) we adjusted the widely-used P < (bar{P}) dry season algorithm^15,16 to identify the arrival and end of the dry season for all the three definitions, which can avoid the influence of short-term climate anomalies. First, we calculated the mean P, Ep, and E for each day (j) of the calendar year (P_j, Ep_j, and E_j) and the daily mean rainfall (bar{P}) for all datasets for 1983-2016. To reduce the synoptic noise, we smoothed P, Ep, and E with a 30-day running window. Then, we calculated cumulative P − Ep, P − E, and P − (bar{P}) on day d, ranging from 1 Jan to 31 Dec, as:

A(d), B(d), and C(d) increase at day d when the daily precipitation is above the daily mean rainfall, daily potential evapotranspiration or actual evaporation, and decrease when the daily precipitation is below the corresponding diagnostic criterion. We defined the day of maximum A(d), B(d), or C(d) as the arrival of the climatological dry season (DSA) and the day of minimum cumulative value as the end of the climatological dry season (DSE). For regions with two or more dry seasons per year, we detected all days of maximum and minimum in the cumulative curve (Fig. 1c), but we used only the four days marking arrivals (dsa1, dsa2) and ends (dse1, dse2) of the two longest dry seasons for our analysis, usually a boreal summer (June–August) dry season and a boreal winter (December–February) dry season.

We calculated DSA, DSE, and DSL under each definition for each dataset (Supplementary Fig. 4), and we calculated the mean DSA, DSE, and DSL under each definition in Fig. 2 and Fig. 3. The mean annual precipitation values and the mean annual temperature values in Fig. 3 were derived from the ERA-5 datasets. We examined the uncertainty by calculating the standard deviation among all ensembles of P, Ep, and E under each definition.

To assess temporal changes, we calculated the arrival and end dates individually for each year from 1983 to 2016. We calculated the cumulative A(d), B(d), and C(d) for each day from DSA − 60 to DSE + 60 for each year instead of the entire calendar year from 1 Jan to 31 Dec, to ensure the correct season was captured. Since the dry season may potentially span multiple calendar years, the dry season arrival and end are not computed for the first and last year of each record. For regions with two dry seasons, the arrival and end dates were determined for the two dry seasons separately. For those regions, we calculated the cumulative function A(d), B(d), and C(d) for each day during DSA1 − 45 to DSE1 + 45 (for the first dry season detection) and DSA2 − 45 to DSE2 + 45 (for the second dry season detection). We used a shorter period (45 days, as opposed to 60 days used for regions with one dry season) in order to better capture the characteristics of the two dry seasons. Accordingly, DSL in days can be calculated as the difference between DSE and DSA, or between DSE1 + DSE2 and DSA1 + DSA2 for cases of two dry seasons. We calculated Water Deficit (WD) as the cumulative sum of P − (bar{P}), P−Ep, or P−E (dashed area in Fig. 1a), from the dates of DSA to DSE.

Long-term trend analysis

To assess temporal changes, we calculated annual dry season diagnostics (DSL, WD, DSA, DSE) individually for each year from 1983 to 2016. We estimated the trends of dry season diagnostics and climatic variables from the ordinary least squares linear regression. We defined each trend as the slope of this linear regression, and we determined statistical significance (P value) using two-tailed Student’s t tests. We used the nonparametric Mann–Kendall trend test to detect whether a significant monotonic increasing or decreasing trend exists, and to provide additional verification for the robustness of the linear regression trend analysis, as it is less sensitive to the beginning and end of the analysis period¹⁶. In addition, we calculated the time series of DSL, DSA, DSE, WD, and meteorological variables at the regional aggregated level using area-weighted averaging over the southern Amazonia, northern and southern central Africa, and southwestern Africa, to maximize large-scale features while minimizing local-scale variability and noise¹⁶. We estimated the linear trends at the regional level as at the grid level (Supplementary Table 2).

Considering the inconsistency of trends across precipitation and evapotranspiration datasets, we judged the level of consistency with the following criterion²⁹: “very likely” if the sign of the trend was the same and significant (P < 0.05) in six to eight precipitation datasets and no significant changes in the others, “likely” if the sign of the trend was supported by four or five precipitation datasets, “probably” if the sign of the trend was supported by one to three, “uncertain” when conflict trends (i.e., both significant increase and decrease trends existed) were found among different precipitation data sources, and “no change” when no significant changes for all of the six datasets were detected. Arid and humid regions (solid gray shaded area in Fig. 4) were excluded when calculating the percent area, since there is no climatology wet or dry season, thus no trends calculated under definitions of P < Ep or P < E.

Driving factors of Ep and E changes

To further illustrate the thermodynamic mechanism driving higher atmospheric water demand, we disaggregated the individual contributions of four meteorological variables (i.e., T, RH, u₂, and Rn) to the Ep trends. We derived the contribution of a certain meteorological variable I to Ep change (C_I) as the difference between the Ep calculated with all variables changing (i.e., ALL) and that calculated with I fixed at its daily climatological values (i.e., I_clim) (Eq. 6). I can be air temperature T, air humidity RH, surface wind speed u₂ or surface net radiation Rn. We calculated the linear trends of Ep and the respective contributions of meteorological variables for the period 1983–2016 (Fig. 5).

As for E, GLEAM estimated this flux through reducing Ep by an evaporative stress factor (S; E = Ep × S + Ei), based on satellite observations of Vegetation Optical Depth (VOD) and assimilated soil moisture⁶⁸. The latter are calculated using a multi-layer running-water balance. Interception loss (Ei) is calculated separately in GLEAM using a Gash analytical model, but its contribution to overall E changes was negligible. Hence, we analyzed changes of these parameters (S, VOD, Soil Moisture) representing the constraints of soil moisture and vegetation water content on evaporation. Daily VOD was derived from VODCA Ku-band⁷⁹, but only available for the period 1987–2016.

Data uncertainties

Due to the insufficient and unevenly distributed observation⁸⁰ in the rainfall data over tropics, we integrated daily meteorological station recode (Supplementary Fig. 13), gauge-based, satellite-combined, and reanalysis datasets to study the variations in precipitation and associated dry season change. Our analyses unravel an overall trend of tropic dry season lengthening and identifying some hotspot regions of changes. However, there are some discrepancies in regions like central Africa and Amazon Basin that may have resulted from data uncertainties and the different approaches used to generate homogeneous climate records. Gauge-based and satellite-combined datasets are quite sensitive to the number and density of observations used, but the observational station is sparse in these regions with large number of missing values in daily record (Supplementary Fig. 13). Different interpolation methods were adopted to fill data gaps and produce grid data, which might have generated errors in the rainfall products. For reanalysis, uncertainties are mainly caused by the biases in reanalyzing models, especially in regions with intricacy land surface process, such as Amazon rainforest.

For the Ep datasets, we used three independent sets of reanalysis data to verify the changes of atmospheric water demand. Our analyses indicate consistent rising in dry-season water demand, which exacerbated the lengthening of tropical dry seasons, from all datasets in southern Amazonia and southern central Africa. However, only a single dataset was used for E, due to the limited data availability at daily intervals, so uncertainty in evapotranspiration estimation has not been fully considered. Actual terrestrial evapotranspiration was modulated not only by surface meteorological conditions and soil moisture but also by the physiology and structures of plants. Changes in nonclimatic factors such as elevated atmospheric CO₂, nitrogen deposition, and land covers also serve as influential drivers. Uncertainties from those complex processes all contributed to the unclear uncertainty in E estimation. Therefore, more efforts should be made to identify and reduce these uncertainties.

Source: Ecology - nature.com