Greening-induced increase in evapotranspiration over Eurasia offset by CO₂-induced vegetational stomatal closure

The Penman–Monteith–Leuning model version 2 (PML-V2) was developed by coupling the widely-used photosynthesis model (Farquhar et al 1980) and a canopy stomatal conductance model (Yu et al 2004) with the Penman–Monteith energy balance equation (Monteith 1965) to jointly estimate GPP and ET (Leuning et al 2008, Zhang et al 2010b, 2016, 2019, Gan et al 2018). We used the PML-V2 for this study is due to that the PML-V2 model shows good performances for considering both effects from rising CO₂, and increasing LAI and climate forcings on ET estimation, compared to the older version PML-V1 or other ET models, e.g. Global Land Evaporation Amsterdam Model (GLEAM) that did not consider rising CO₂'s impact on canopy stomatal conductance. The PML-V2 model has been applied to successfully produce the MODIS LAI-based global GPP and ET products at a 500 m and 8 d resolution since 2002, which were noticeably better than most widely used diagnostic GPP and ET products (Zhang et al 2019).

The PML-V2 model consists of three components for separately calculating the ET's three components (equations (1)–(3)): canopy transpiration (${E_c}$), soil evaporation (${E_s}$) and rainfall interception (${E_i}$),

$$\begin{equation}{\text{ET}} = {E_c} + {E_s} + {E_i}.\end{equation} \tag{ 1 }$$

The PML-V2 model simulates ${E_i}$ (see supplementary text) based on a revised Gash-model scheme (van Dijk and Bruijnzeel 2001) with key parameters (see table S1 (available online at stacks.iop.org/ERL/16/124008/mmedia)) calibrated using filed measurements. The canopy transpiration (${E_c}$) is estimated using the Penman–Monteith energy balance equation but the canopy stomatal conductance (${G_c}$) is not a constant and changes with vegetational dynamics (equations (2) and (4)) through coupling with the photosynthesis process (${A_{gs}}$). The photosynthesis process (${A_{gs}}$) changes with increasing CO₂ and LAI across different vegetation types, and then produces ${G_c}$ (equations (4)–(6)) and GPP (see supplementary text). The soil evaporation (${E_s}$) depends on absorbed energy flux and soil water deficit,

$$\begin{equation}{E_c} = \frac{{\varepsilon {A_c} + \;\left( {\rho {c_p}/\gamma } \right)D{G_a}}}{{\varepsilon + 1 + {G_a}/{G_c}}}\end{equation} \tag{ 2 }$$

$$\begin{equation}{E_s} = \frac{{f\varepsilon {A_s}}}{{\varepsilon + 1}}\end{equation} \tag{ 3 }$$

where the surface available energy ($A = {R_n} - G$) is divided into canopy absorbed energy (${A_c}$) and soil absorbed energy (${A_s}$), ${A_c} = \left( {1 - \tau } \right)A$ and ${A_s} = \tau A$, $\tau = {\text{exp}}\left( { - {k_A}{\text{LAI}}} \right)$, ${k_A} = 0.6$; ${R_n}$ is the net radiation and $G$ is the ground heat flux (W m⁻²). $\varepsilon = \frac{{{\Delta }}}{\gamma }$, and ${{\Delta }}$ is the slope of the curve relating saturation water vapor pressure to temperature (kPa °C⁻¹). $\rho$ is the density of air (g m⁻³); ${c_p}$ is the specific heat of air at constant pressure (J g⁻¹ °C⁻¹); $D$ is vapor pressure deficit (kPa). ${G_a}$ is the aerodynamic conductance (m s⁻¹); ${G_c}$ is the canopy conductance (m s⁻¹); $f$ is a dimensionless variable that determines the water availability for soil evaporation, which is a fraction function of the accumulated precipitation and soil equilibrium evaporation for each grid cell (Zhang et al 2019).

In PML-V2 model, the canopy stomatal conductance (${G_c}$) are calculated by the photosynthesis process (${A_{gs}}$) for each plant functional type (PFT) with the constraint of vapor pressure deficit ($D$) at surface,

$$\begin{equation}{G_c} = \frac{{m{A_{gs}}{ }}}{{{C_a}\left( {1 + D/{D_0}} \right)}}\end{equation} \tag{ 4 }$$

$$\begin{align}{A_{gs}} = & \frac{{{P_1}{C_a}}}{{k\left( {{P_2} + {P_4}} \right)}} \nonumber \\ & \times\left\{ {k{\text{LAI}} + {\text{ln}}\frac{{{P_2} + {P_3} + {P_4}}}{{{P_2} + {P_3}\exp \left( {k{\text{LAI}}} \right) + {P_4}}}} \right\}\end{align} \tag{ 5 }$$

$$\begin{equation}{V_m} = \frac{{{V_{m,25}}\exp \left[ {a\left( {T - 25} \right)} \right]}}{{1 + {\text{exp}}\left[ {b\left( {T - 41} \right)} \right]}}\end{equation} \tag{ 6 }$$

where ${P_1} = {A_m}\beta {I_0}\eta$, ${P_2} = {A_m}\beta {I_0}$, ${P_3} = {A_m}\eta {C_a}$, ${P_4} = \beta {I_0}\eta {C_a}$, and ${A_m} = 0.5{V_m}$. ${I_0}$ is the photosynthetically active radiation (PAR, in mol) from shortwave downward radiation. $\,{C_a}$ is the atmospheric CO₂ concentration (in ppm or mol/mol). ${V_{m,25}}$, $\beta$, $\eta$, $m$ are key parameters for photosynthesis process, i.e. maximum photosynthetic carboxylation efficiency (${V_{{\text{cmax}}}}$), initial photochemical efficiency ($\beta$), initial value of the slope of CO₂ response curve ($\eta$), and stomatal conductance coefficient ($m$), and D_min, D_max, ${D_0}$ are parameters to identify the constraint of atmospheric vapor pressure deficit in PML-V2 model (see table S1). All values of these key parameters for global ten major PFTs in table S1 were calibrated using half-monthly GPP, ET and climate forcing over 2002–2014 derived from 95 eddy covariance (EC) flux sites of the FLUXNET2015 dataset and AVHRR GIMMS3g-based LAI. The PFTs are defined from the NASA MCD12Q1 IGBP land cover classes (Friedl et al 2010). The detailed calibration method has been shown in Zhang et al (2019), which optimize the key parameters for each PFT by minimizing the objective function (F): $F = 2 - \left( {{\text{NS}}{{\text{E}}_{{\text{ET}}}} + {\text{NS}}{{\text{E}}_{{\text{GPP}}}}} \right)$, where ${\text{NS}}{{\text{E}}_{{\text{ET}}}}$ and ${\text{NS}}{{\text{E}}_{{\text{GPP}}}}$ are the Nash–Sutcliffe Efficiency (NSE) for ET and GPP, respectively, which are obtained by comparing the PML-V2 simulations and FLUXNET2015-based observations.

Using the parameter-calibrated PML-V2 model, we estimated grid cell-scale half-month ET and GPP using observational climate forcing, AVHRR GIMMS3g-based LAI and CO₂ data, and compare them to observations at 95 eddy flux sites from FLUXNET2015-based dataset (figure S1). We show that compared to the FLUXNET2015-based observations, the PML-V2 model presents a good model performance with e.g. a NSE of 0.70, a root mean square error (RMSE) of 0.7, and a bias of −0.21% for half-month ET; and a NSE of 0.74, a RMSE of 2.05, and a bias of 3.8% for half-month GPP (figures S1(a) and (b)). Even under the cross-validation mode of the PML-V2 model, the comparing results are only slightly degraded, e.g. the NSE becomes 0.6 for ET and 0.61 for GPP (figures S1(c) and (d)).

In this study, we applied the PML-V2 model to produce historical coupled estimates of GPP, ET and its components ${E_c}$, ${E_s}$ and ${E_i}$ over 1982–2014 at 0.05° and half-monthly resolutions. The forcing inputs for the PML-V2 model are land surface climate variables (i.e. precipitation, near-surface air temperature, pressure and relative humidity, wind speed at a 10 m height, downward longwave and shortwave radiations) from the GLDAS V2.0 observational 3 h 0.25° climate dataset (Rodell et al 2004, Beaudoing and Rodell 2019), the global satellite-based half-monthly LAI from the satellite AVHRR GIMMS3g V3 product (Zhu et al 2013), and the global satellite-based 8 d surface albedo for shortwave and surface emissivity for longwave radiation from the GLASS V4 product (Liu et al 2013, Cheng and Liang 2014). Atmospheric CO₂ concentration data are provided by NOAA global monthly mean atmospheric CO₂ records (ftp://aftp. cmdl.noaa.gov/products/trends/co2/co2_mm_gl.txt).

To isolate individual environmental drivers' contribution to historical changes in ET, ${E_s}$, ${E_c}$ and ${E_i}$, we designed a group of modeling sensitivity experiments (see table 1) that are similar to those used in the Multi-Scale synthesis and Terrestrial Model Intercomparison Project (Huntzinger et al 2013). Using the PML-V2 model driven by observational climate forcing, observed atmospheric CO₂, satellite-derived LAI, and albedo and emissivity datasets, we conducted five scenarios of sensitivity simulations (S1–S5, see table 1) by setting all drivers change with time (S1), or one or more drivers using fixed climate at 1982 (constant-1982) with others dynamic (S2–S4), or all constant-1982 drivers (S5 as the control run) as inputs over 1982–2014. Therefore, differences between pairs of these sensitivity simulations allow for a robust assessment of model-derived carbon and water fluxes (e.g. ET) sensitivity (including variability and trend) to individual drivers (S1 and S2 for albedo and emissivity's contribution (equation (7)), S2 and S3 for LAI's contribution (equation (8)), S3 and S4 for CO₂'s contribution (equation (9)) and S4 and S5 for climate's contribution (equation (10))),

$$\begin{align}\Delta {\text{E}}{{\text{T}}_{{\text{Ald}}\, \&\, {\text{Ems}}}} = {\text{E}}{{\text{T}}_{{\text{S1}}}} - {\text{E}}{{\text{T}}_{{\text{S2}}}},\end{align} \tag{ 7 }$$

$$\begin{align}\Delta {\text{E}}{{\text{T}}_{{\text{C}}{{\text{O}}_{\text{2}}}}} = {\text{E}}{{\text{T}}_{{\text{S2}}}} - {\text{E}}{{\text{T}}_{{\text{S3}}}},\end{align} \tag{ 8 }$$

$$\begin{align}\Delta {\text{E}}{{\text{T}}_{{\text{LAI}}}} = {\text{E}}{{\text{T}}_{{\text{S3}}}} - {\text{E}}{{\text{T}}_{{\text{S4}}}},\end{align} \tag{ 9 }$$

$$\begin{align}\Delta {\text{E}}{{\text{T}}_{{\text{CLM}}}} = {\text{E}}{{\text{T}}_{{\text{S4}}}} - {\text{E}}{{\text{T}}_{{\text{S5}}}}.\end{align} \tag{ 10 }$$

Interactions between two or more drivers (S2−S4 for LAI and CO₂'s contribution, S3−S5 for climate and CO₂'s contribution, S2−S5 for LAI and CO₂ and climate's contribution) can also be estimated. In this study, we estimated the overall interaction that is calculated by difference between the historical change (S1−S5) and the sum of individual driver's direct contribution (equations (7)–(10)). It is noted that we did not separate variables of climate forcing and therefore make our modeling sensitivity experiments smaller and results easy to interpret.

We calculated the relative contributions of individual drivers (table 1) to the global or regional long-term changes or trends in the studying variables (e.g. annual ET) by

$$\begin{align}& {\text{Relative contribution of }}{X_i}{ } \nonumber \\ & = \frac{{{X_i}\,\,{\text{induced change or trend in }}Y}}{{\mathop \sum \nolimits_{i = 1,4} \left| {{X_i}\,\,{\text{induced change or trend in }}Y} \right|{ }}} \times 100\% \end{align} \tag{ 11 }$$

where ${X_i}$ is one of the four major model drivers (climate change, CO₂, LAI, surface albedo and emissivity), and $Y$ is a response variable (e.g. ET, ${E_c}$, ${E_s}$ or ${E_i}$). Equation (11) indicates that a negative (positive) value of the relative contribution of ${X_i}$ means that driver ${X_i}$ has a negative (positive) contribution to the overall change or trend in the response variable $Y$. Therefore, the overall contribution should be summed from all absolute values of individual relative contributions (i.e. $\mathop \sum \limits_{i = 1,4} \left| {{\text{Relative contribution of }}{X_i}} \right|$) to be sure that the summed result is 100%.

2.4.1. Observational annual ET derived from basin-scale water balance

At a large scale (e.g. basin-scale in this study), the ET can be easily estimated by the water-balance method, considered as the difference between precipitation (P) and basin streamflow (Q) and change in water storage (ΔS), that is, ET_wb = P − Q − ΔS. For a given basin, P and Q are observed from meteorological or hydrological stations, and ΔS has been widely estimated from the gravity recovery and climate experiment (GRACE) satellite-based gravity field. Therefore, the derived ET_wb is often regarded as 'observational' to evaluate model performance (Zeng et al 2014, Han et al 2015, Ma et al 2019). In this study, we selected 15 large basins over Eurasia where annual ET_wb over 2003–2014 is available and each has a contribution area >40 000 km² to minimize uncertainties from the valid coarse spatial resolution of GRACE-based ΔS) (figure 1(a)). Specifically, annual ET_wb, which was from a derived dataset (Ma et al 2021), was calculated by the GPCC 0.25° × 0.25° observational P (Becker et al 2013), gauged Q and the GRACE‐observed annual ΔS (available at http://grace.jpl.nasa.gov) (Tapley et al 2004, 2019). For each basin, the annual ET datasets described in section 2.4.2 were calculated from an area-weighted averaging method over all grid cells within the basin.

Zoom In Zoom Out Reset image size

2.4.2. Three types of long-term ET products

Firstly, we compared the historical annual ET estimated by PML-V2 with the observational annual ET_wb derived over 2003–2014 from basin-scale water-balance method, and the result indicates a good performance of PML-V2 ET at basin scales (figure 1). Across 15 river basins in Eurasia, the PML-estimated annual ET are strongly correlated with the annual ET_wb with a coefficient of determination (R²) of 0.70 (figure 1(b)). Due to different climate forcing used in the ET datasets, significant positive bias (∼11%) for all basins between the annual ET_wb and the PML ET can be found (figure 1(b)), especially for the Yangtze River basin (∼30%) and the basins with small annual ET (e.g. the Kolyma, the Lena and the Yenisey rivers with biases of 50%–100%).

Secondly, we evaluated trends in annual ET over 2003–2014 (figure 1(c)). There are 12 out of 15 basins selected for trend analysis, since three basins (i.e. the Lena, the Songhua, and the Volga rivers) do not have complete annual streamflow data over 2003–2014. For the ET_wb trends across the 12 basins, 8 out of the 12 have strong increasing trends, ranging from 0.8 to 6.0 mm year⁻², while the basins with strong negative trends mainly locates at the Europe (e.g. the Danube and the Rhine rivers). In comparison, the PML-estimated ET in 10 out of the 12 basins have the same directions in trends of ET_wb, with an R² of 0.41 between trends in the two ET datasets (figure 1(c)). Relatively large differences (>50%) in trends between PML-V2 ET and ETwb can been seen for the basins of the Danube, Elbe, Ob, PRYPYAT, Yangtze and Zhujiang rivers (figure 1(c)). With area-weighted averaging applied over the 12 basins, the PML-V2 ET shows a trend of 1.83 mm year⁻² over 2003–2014, which is close to the ET_wb trend (1.57 mm year⁻²). These results indicate that the PML-estimated ET is adequate for further analysis of historical changes in ET and its driving factors' contributions over Eurasia.

With continually increased anthropogenic carbon emissions, global atmospheric CO₂ has increased by 15.5% (∼56 ppm) over 1982–2014. In response to the CO₂-induced climate change and vegetational fertilization, the annual LAI (or NDVI) over Eurasia has increased by 7.36% (or 5.79%) during this period (figure S2). The most significant increasing trends in LAI (>15%) or NDVI (>10%) can be found over the Europe, the Indian Peninsula and southeast China, occupying about 35% of the Eurasian area (figures 2(a), (b) and S3(a), (b)). The climate forcing over Eurasia has also experienced remarkable change (figures 2(c)–(f)), with surface air temperature (in °C) increased by 13.63%, and precipitation increased by 5.42% and other variables less than 5% (figure S2). The temperature increases exceed 20% over about 30% of the Eurasian land (figure S3(c)).

Zoom In Zoom Out Reset image size

In response to these increases in surface driving variables, the annual ET over Eurasia also has largely increased. By comparing the historical simulated ET trends from PML-V2 with the ERA5, GLEAMv.3.3a, and CRv1.0 ET products (figure 3), we find that all four ET datasets show significant increasing trends in ET over Eurasia for 1982–2014 (figure 3(a)). On area-weighted average over Eurasia, the GLEAM ET shows the largest trend (10.89 ± 0.90 mm year⁻¹ decade⁻¹, p < 0.01), followed by the PML ET trend (6.20 ± 1.13 mm year⁻¹ decade⁻¹, p < 0.01) and the CR ET trend (5.21 ± 1.15 mm year⁻¹ decade⁻¹, p < 0.01). The ERA5 ET trend is the smallest, but still strong (3.70 ± 1.05 mm year⁻¹ decade⁻¹, p < 0.01). It is worthy noted that uncertainty in trends between the four ET products may be due to different climate forcing used, but other important drivers (e.g. LAI, CO₂) with large long-term trends whether being considered or not as input in the ET algorithms could strongly alter long-term change/trend of regional ET. In a word, the PML ET trend is close to the ensemble mean (6.51 ± 3.1 mm year⁻¹ decade⁻¹) of the other three. Spatial patterns of trends in ET are generally consistent within the four datasets (figures 3(c)–(f)). Most of high increasing trend values are from highly vegetated regions (about 60% of Eurasian area) over central and eastern Europe, the Indian Peninsula and southeast China where vegetation greening trends are largest (figures 2(a) and (b)). The ET trends over these high greening regions are larger than 20 mm year⁻¹ decade⁻¹ for both the PML-V2 and GLEAM datasets.

Zoom In Zoom Out Reset image size

Using the historical estimations from PML-V2 driven by all observational forcings, we quantified the relative contributions of the three ET components, i.e. ${E_c}$, ${E_s}$, ${E_i}$ to the long-term changes of annual ET during 1982–2014 over Eurasia (figures 4 and S4). On area-weighted average over the whole Eurasia, we estimated that the long-term increase in annual ET comes from the increases in ${E_c}$ and ${E_i}$, and the decrease in ${E_s}$ (figure 4(a)). The long-term trend in ET over Eurasia (6.20 mm year⁻¹ decade⁻¹) is primarily dominated by the increases in ${E_c}$ (4.18 mm year⁻¹ decade⁻¹) with a relative contribution of 61.1%. Although ${E_s}$ has a larger interannual variation over 1982–2014 than ${E_i}$, ${E_s}$ has a much smaller contribution (−4.47%) to the trend of annual ET, while ${E_i}$ has a 34.46% of trend contribution (figure 4(b)). Regionally, the strong increase of ET over Eurasia is mainly from the highly vegetated regions where ${E_c}$ dominantly contributes 60%–80% of the ET trend (figures 4(c) and (d)). Over most of these regions, ${E_s}$ shows a negative contribution of 0%–20%, and ${E_i}$ contributes 0%–40% of reginal ET trends (figures 4(e) and (f)). For the desert or sparse vegetated regions such as the Arabian Peninsula, the Middle East, and the Qinghai–Tibet Plateau, ${E_s}$ dominates 80%–100% of the regional ET (positive or negative) trends (figures 4(c) and (e)).

Zoom In Zoom Out Reset image size

With the modeling sensitivity experiments using PML-V2 (see table 1), we further isolate the relative contributions of individual environmental drivers, i.e. climate forcing, CO₂, LAI, albedo and emissivity, and their interactions (figures 5 and S5). We find that on area-weighted average, the summed effect of all individual drivers' contribution explains ∼87% of the trend in annual ET during 1982–2014 (figure 5). The nonlinear effect of driver interactions presents a negative contribution (∼13%) to the overall ET trend (figures 5(a) and (b)). Noticeably, we find that the timeseries of climate forcing-induced change in annual ET almost overlaps the historical change in annual ET driven by all drivers (figure 5(a)), indicating the dominant role of climate on both interannual variation and trend of ET over Eurasia. Due to the strong greening (an increase of 7.36% over Eurasia) over 1982–2014 as observed by satellites, the LAI alone has induced a strong increasing trend (3.88 mm year⁻¹ decade⁻¹) in annual ET. In contrast, the increased CO₂ has resulted in a substantial decreasing trend (−4.67 mm year⁻¹ decade⁻¹) in annual ET (figures 5(a) and (b)), due to the CO₂-induced suppression effect on ${E_c}$ through reducing vegetational stomatal conductance (equation (4)). Consequently, the historical increasing trend in ET over Eurasia is mainly due to the positive relative contributions of changes in climate forcing (39.55%), LAI (21.94%) and surface albedo and emissivity (6.06%), and the negative relative contribution of increased CO₂ (26.40%). Regionally, the dominance of climate forcing on ET trend covers 68.26% of Eurasian land area except for the South Asia and a part of western Europe where CO₂ dominates (9.38%), and the India Peninsula and a part of northern Siberia where the LAI controls (13.06%, see figure 5(f)). Over most vegetation greening regions (figures 5(c)–(f)), such as central and eastern Europe and southeast China, the LAI-induced ET increasing trends (5–10 mm year⁻¹ decade⁻¹) are almost offset by CO₂-induced ET decreasing trends (−10 to −5 mm year⁻¹ decade⁻¹). Even over the India Peninsula where LAI dominates the ET trends (15–20 mm year⁻¹ decade⁻¹), the CO₂ has induced ET strong decreasing trends (−15 to −10 mm year⁻¹ decade⁻¹).

Zoom In Zoom Out Reset image size