Significance The significance of forest mortality on ecosystem services, and water, carbon, and nutrient cycling is indubitable. While there is a general agreement that climate change-induced heat and drought stress is expected to intensify forest mortality, the concurrent influence of changes in atmospheric humidity and CO 2 concentration remains unclear. Here, the response of mortality risk to projected climate change is evaluated in 13 biomes across the globe. Our results show that increasing specific humidity and CO 2 concentration partially offset the intensification of risk by changing precipitation and air temperature. The risk response is also mediated by plant hydraulic traits. The study provides a mechanistic foundation for estimating future responses of forest mortality risk, which can facilitate ecosystem management.

Abstract Climate-induced forest mortality is being increasingly observed throughout the globe. Alarmingly, it is expected to exacerbate under climate change due to shifting precipitation patterns and rising air temperature. However, the impact of concomitant changes in atmospheric humidity and CO 2 concentration through their influence on stomatal kinetics remains a subject of debate and inquiry. By using a dynamic soil–plant–atmosphere model, mortality risks associated with hydraulic failure and stomatal closure for 13 temperate and tropical forest biomes across the globe are analyzed. The mortality risk is evaluated in response to both individual and combined changes in precipitation amounts and their seasonal distribution, mean air temperature, specific humidity, and atmospheric CO 2 concentration. Model results show that the risk is predicted to significantly increase due to changes in precipitation and air temperature regime for the period 2050–2069. However, this increase may largely get alleviated by concurrent increases in atmospheric specific humidity and CO 2 concentration. The increase in mortality risk is expected to be higher for needleleaf forests than for broadleaf forests, as a result of disparity in hydraulic traits. These findings will facilitate decisions about intervention and management of different forest types under changing climate.

Forest mortality can lead to irreversible change in vegetation cover, thereby affecting many processes pertinent to water, carbon, and nutrient budgets (1). Multiple studies (2⇓⇓⇓⇓⇓⇓⇓–10) have noted close association between forest mortality and water and heat stress, owing to shifting precipitation patterns and rising air temperature. However, the influence of concurrent changes in specific humidity (SH) and CO 2 concentration, which affect plant response to stress by altering stomatal kinetics (11), have not received similar attention. Although elevated CO 2 concentration is expected to promote future forest productivity (12), the extent to which it affects forest mortality in the context of water and heat stress remains a subject of inquiry. Short-term records (3, 4) and long-term manipulative field studies in forests such as the Free Air CO 2 Enrichment experiments (13⇓–15) have tried to fill the knowledge gap; however, they do not cover the entire manifold of projected climate conditions. The goals of this study are to evaluate the individual and combined influence of projected changes in precipitation, temperature, SH, and CO 2 concentration on forest mortality risk and to investigate whether the response of mortality risk differs among plant functional types (PFTs).

Tree mortality may occur through several mechanisms, including hydraulic failure, carbon starvation, phloem transport limitation, and biotic attack (16, 17). Hydraulic failure is characterized as the malfunction of xylem water transport associated with cavitation, which is induced by low xylem water potential under limited soil water availability. Carbon starvation occurs when carbohydrate supply and storage cannot meet demand (17), which could result from low photosynthesis due to stomatal closure in response to low plant water potential and high atmospheric vapor pressure deficit (VPD). Reduced photosynthesis and plant water potential also pose limitations for phloem to maintain turgor pressure and may further impair phloem transport (18). Intense and prolonged stresses could weaken the defenses of forests to biotic attack (5) and may alter plant adaptation, seed production, and germination (2). Despite these mechanisms being far from thoroughly understood (17, 18), they primarily result from low plant water potential and restricted photosynthesis.

To quantify the risk of mortality induced by low plant water potential, previous studies (19, 20) used the safety margin [i.e., the difference between the minimum observed xylem water potential and the xylem water potential at 50% loss of conductivity ( ψ 50 )]. Plants with narrower or more negative safety margins are considered to be more susceptible to hydraulic failure. The safety margin provides a static assessment of plant susceptibility to hydraulic failure, although its representativeness may be undermined by limited field observations. It has also been suggested that, instead of the minimum water potential plants reach, the duration plants operate under high percentage loss of conductivity could more likely distinguish mortality (21, 22). Here, a duration-based hydraulic failure risk (HFR) is introduced, which quantifies the fraction of days when the daily minimum xylem water potential ( ψ x , min ) falls below ψ 50 . Because stomatal closure restricts photosynthesis (6, 17, 23), a stomatal closure risk (SCR) can also be formulated as the fraction of days on which stomata are completely closed (SI Appendix, section S1). The aggregated mortality risk is then defined as the fraction of days with occurrence of either hydraulic failure or stomatal closure, two physiological states contributing to dieback and eventual mortality. Alternative quantifications of risk that account for stress duration and severity are also considered to test the robustness of the analysis here (SI Appendix, section S5).

The mortality risk is evaluated by using a soil–plant–atmosphere continuum (SPAC) model, which computes hourly dynamics of xylem water potential and stomatal conductance (Materials and Methods and SI Appendix, section S1). By comparing against observed mortality at four sites, the modeled risk is shown to capture the temporal variation of mortality in response to climate stressors (SI Appendix, section S3). The risks under historical and future climate scenarios are then evaluated for 13 temperate and tropical forest biomes across the globe (Fig. 1). The selected biomes cover a variety of climates and PFTs, including evergreen needleleaf forest (ENF), deciduous broadleaf forest (DBF), and evergreen broadleaf forest (EBF). The response of mortality risk to changes in the following climate characteristics is analyzed: mean annual precipitation (MAP), precipitation seasonality (PS), mean annual air temperature (T), mean annual atmospheric SH, and atmospheric CO 2 concentration. PS is quantified as the fraction of MAP that falls within the growing season. Changes in these climate conditions are obtained from the projection of Coupled Model Intercomparison Project Phase 5 (CMIP5) (Table S2) from 1986–2005 to 2050–2069.

Fig. 1. Distribution of PFT and locations of the 13 investigated biomes. Biomes are the areas within the selected rectangular regions that are covered by a given PFT. PFTs shown in the map include evergreen broadleaf forest (EBF), deciduous broadleaf forest (DBF), evergreen needleleaf forest (ENF), and deciduous needleleaf forest (DNF).

Influence of Combined Changes in Climate Variables Based on the CMIP5 projections of four RCP scenarios in all 13 biomes (SI Appendix, Fig. S6), the response of mortality risk to changes in three combinations of climate conditions are examined (Fig. 3): (i) P+T; (ii) P+T+SH; and (iii) P+T+SH+CO 2 . Here, changes in P include combined changes in MAP and PS. For the 13 investigated biomes, on average, shifting precipitation patterns and rising temperature projected by RCP4.5 are found to intensify the risk by 158.8% for the period 2050–2069 relative to the historical risk. This increase in risk is consistent with previous studies highlighting the exacerbating effects of higher temperature (1, 3⇓⇓⇓–7). However, by incorporating increases in SH, the risk decreases by 46.6%. More remarkably, the risk drops an additional 91.2% under the added influence of elevated CO 2 concentration. In aggregate, changes in all four climate conditions increase the risk by 21.0% on average, which is much lower than the increase of 158.8% when only the changes in precipitation and air temperature are considered. Under high emission scenarios (RCP6.0 and RCP8.5), elevated humidity and CO 2 concentration might even overwhelm the effects of higher temperature, possibly resulting in a lower risk than the historical level (Fig. 3). These alleviating effects are robust across alternative risk measures (SI Appendix, Table S6 and Fig. S14). The alleviating effect of increasing atmospheric CO 2 concentration is in line with a reported decrease in stomatal conductance and increase in water use efficiency across various climate regions and species (2, 25⇓⇓⇓–29). Fig. 3. Response of mortality risk to combined changes in climate variables. P+T includes changes in MAP, PS, and air temperature; P+T+SH includes additional changes in SH; and P+T+SH+CO 2 includes additional changes in atmospheric CO 2 concentration. Gray dashed lines in the subplots show the risks under historical climates. Upper and lower boundaries of the boxes correspond to the 25th and 75th quantiles of the risk based on multimodel projections in each RCP scenario. Numbers in the subplots correspond to the biome locations in Fig. 1. On average, the combined changes of P+T+SH+CO 2 in RCP4.5 are found to increase the mortality risks by 101.1%, −18.3%, and 19.6% for ENF, DBF, and EBF biomes respectively (Fig. 4). The significantly higher increase for ENF compared with the other two PFTs results from their distinct risk sources (Fig. 4). SCR, the primary risk source for ENF, shows notably higher sensitivity to air temperature rise than HFR. Specifically, for a 1 °C increase in air temperature from historical climates, HFR and SCR are estimated to increase by 23.5% and 125.1%. Remarkably, the increase in SCR is close to the previously reported 116.3% °C−1 increase in die-off events of Pinus edulis (4), a conifer species likely threatened by SCR. From a mechanistic perspective, elevated temperature increases VPD and reduces stomatal conductance. This restricts carbon assimilation but promotes water loss, which results in a higher probability of full stomatal closure (i.e., higher SCR). The increase in HFR is smaller, as the reduction in stomatal conductance partly alleviates the increase in water loss due to increased VPD. Under projected changes of precipitation and VPD in RCP4.5, HFR increases by 135.5% on average, while SCR increases by 305.8%. When CO 2 is also considered, the aggregate changes are −8.6% and 83.7% for HFR and SCR, respectively. These findings imply a larger increase in mortality risk of ENFs, more specifically of isohydric species, than other PFTs under changes in the considered climate conditions. Fig. 4. Average risk under historical climates (open circles), changes in P+T (light filled circles), and changes in P+T+SH+CO 2 (dark filled circles) for each PFT based on RCP4.5. Gray filled triangles denote the relative contribution of hydraulic failure (HF) to the mortality risk (i.e., 1 denotes that all of the risk is due to HF, and 0 denotes all of the risk is from stomatal closure).

Discussion and Implications The study evaluates how projected climate change will affect mortality risks and how the risks may be mediated by different PFTs across the globe. In this regard, the study introduces a measure of mortality risk that accounts for the duration that plants operate under high percentage loss of conductivity or stomatal closure. Although large uncertainty exists in the exact physiological mechanisms that cause mortality (17, 18), the proposed mortality risk measure captures two of the fundamental causes (i.e., low water potential and severely restricted carbon assimilation) which contribute to the downstream mortality mechanisms. Notably, the quantification of mortality is made possible by synergistic coupling of multiple prior submodels connecting plant physiological status to hydrological and meteorological conditions. This coupling allows the SPAC model to resolve hourly dynamics of xylem water potential and stomatal conductance, the variables required to evaluate mortality risk (SI Appendix, sections S1 and S5). For example, by accounting for the feedback between evapotranspiration and atmospheric boundary layer (ABL) development, the SPAC model is able to simulate a physically consistent hourly dynamic of air temperature and SH during droughts. During drought, when evapotranspiration is restricted by low soil moisture, the model partitions a larger fraction of incoming energy into sensible heat, thus enhancing the ABL and raising the temperature during daytime. This evapotranspiration–ABL coupling allows SPAC to consider the cooccurrence of extreme drought and heat stress, which has been pointed out as the main environmental trigger of tree mortality (3⇓–5, 7). The SPAC model also uses an optimization-based stomatal conductance representation that accounts for the effect of plant hydraulic limitation (SI Appendix, section S1). The representation is an advantage over several widely used dynamic global vegetation models, where the stomatal regulation is disconnected or empirically connected with soil water stress (6, 21). Many of these semiempirical models are derived from observations under ambient CO 2 concentration, and their parameter values are subject to change in an elevated CO 2 environment (11), thus undermining their efficacy under future climate. In contrast, the optimization-based stomatal regulation model used here has been demonstrated to predict stomatal response to stress under both historical and elevated CO 2 concentration (30). Although the modeled risk shows strong correspondence with observed mortality (SI Appendix, section S3), the estimated risk should be interpreted with care. Given the uncertainties inherent in model structure and parameters, and the complexity of the forest ecosystem, it is unrealistic to accurately assess the exact magnitude of mortality risk. Large variations in plant hydraulic traits, tree height, diameter at breast height, and stand density (31, 32) may impact the actual mortality risk. Sensitivity of mortality risk to aforementioned factors and to uncertainties in model structure are examined (SI Appendix, section S4). Results indicate that, despite their influence on the magnitude of mortality risk, the alleviating effect of increasing SH and CO 2 concentration is still found to be robust. Notably, actual mortality risk may also be altered by forest fire frequency and insect outbreak, rooting profile, seed production, community-level competition, local acclimation to drought, and adaption to long-term climate change (33⇓⇓–36), factors whose characterization is still fraught with uncertainties (37). Their impacts in relation to the direct influence of climate conditions discussed here deserve further investigation. However, independent of these indirect influences, results reported here demonstrate a ubiquitous and robust alleviating effect of elevated atmospheric humidity and CO 2 concentration, which is comparable in magnitude to the intensifying effect of changes in precipitation patterns and air temperature. The combined influence of changes in these climate variables on mortality risk is also strongly mediated by plant hydraulic traits. These results highlight that ignoring the influence of elevated atmospheric humidity and CO 2 concentration may lead to overestimation of future forest mortality risk.

Materials and Methods SPAC Model. The SPAC model consists of three process components: a soil–water balance; a plant water transport that is based on cohesion–tension theory and associated hydraulic properties; and an ABL development model that permits evapotranspiration to alter the height, temperature, and SH of the boundary layer (SI Appendix, section S1 and Fig. S1). Soil is characterized as a two-layer bucket, where the soil moisture is controlled by precipitation, interception, vertical flux between the two soil layers, deep percolation, soil evaporation, and plant root extraction. Water transport within plants is modeled as a resistance system with no capacitance. Water vapor and CO 2 exchange at the leaf level are modeled by combining Fickian diffusion of gases and the Farquahar photosynthesis model (38, 39), where the stomatal kinetics are determined by optimizing carbon gain while minimizing water losses (11). The stomatal conductance is affected by both atmospheric conditions and plant water status. The water flux through the soil–plant–atmosphere system is solved given soil and atmospheric conditions. Subdaily temperature and SH are obtained from the ABL development model. The energy and mass components in ABL development are affected by feedback of total water flux from the ground surface (40), including interception, soil evaporation, and plant transpiration. The coupled SPAC model simulates ecohydrologic states at an hourly interval. Soil and Vegetation Properties. Based on the global land cover type from the moderate imaging spectroradiometer (MODIS; MCD12C1) (41), 13 forest biomes were selected across the globe. The biomes cover three PFTs and a variety of climate types (Fig. 1 and SI Appendix, Table S1) with MAP ranging from 500 to 3,000 mm. Regions with shallow groundwater (42) and snow-dominated climate (43) were avoided, as the influence of groundwater and snow is not considered in the model. The SPAC model in each biome was parameterized with local soil and representative plant properties. Soil texture compositions were obtained from the Harmonized World Soils Database (44). Soil hydraulic properties were calculated based on the generalized statistical relations (45). The annual cycle of the Leaf Area Index (LAI) was extracted from the level-4 MODIS global LAI and Fraction of Photosynthetically Active Radiation product (MCD15A2) (46). Plant hydraulic traits were obtained from a global database containing hydraulic traits of 866 species (47). Photosynthetic parameters were derived from a cross-species study (48). Stomatal optimization parameters were obtained based on a metaanalysis study across PFTs and climates reported in a previous study (49). These plant properties were obtained at a biome level by averaging the properties of trees belonging to the same PFT and climate type (47) as found in the given biome. Historical and Projected Climates. The SPAC model in each biome is forced by local daily climate, including stochastic precipitation, net shortwave radiation, and initial and boundary conditions of potential temperature and SH of ABL. At the beginning of each day, the ABL is reset with the corresponding initial and boundary conditions. The stochastic precipitation is represented as a marked Poisson process characterized by frequency and mean rainfall depth statistics (50). Daily historical climates are calculated based on the National Centers for Environmental Prediction (NCEP)/National Center for Atmospheric Research (NCAR) reanalysis data (51) from 1986 to 2005. Projected climate changes are obtained from multimodel outputs of CMIP5 experiments under four RCP scenarios (Table S2). For each model under each RCP scenario, changes in climate variables are quantified as the difference between the averages for 1986–2005 and 2050–2069. Future precipitation statistics and other climate forcings for the model are generated by incorporating these changes into the historical climates from NCEP/NCAR (SI Appendix, section S2) to eliminate the influence of biases in climate model outputs (52). Historical and future atmospheric CO 2 concentrations under the four RCP scenarios are provided in ref. 53. Experimental Design and Statistics. The mortality risk under a given climate was quantified based on plant dynamics by running the SPAC model at hourly resolution for 30 annual ensembles after a 5-y warm-up period. Influence of changes in each individual climate variable was analyzed by keeping the others the same, while only changing the target climate variable. Influence of combined changes in climate variables as projected by multimodels were grouped together to evaluate the overall response of mortality risk under each RCP scenario. Each reported change in risk is the average of changes of all of the biomes, unless stated otherwise. Change in each biome was calculated as the difference between the historical risk and the future risk (i.e., the average risk based on multimodel projections), in proportion to the historical risk. Biomes with historical risks <0.01% were excluded from the statistics. Note. Details on model formulation, data processing, model validation, sensitivity analyses, and alternative quantifications of risks are provided in SI Appendix.

Acknowledgments We thank the National Aeronautics and Space Administration Earth Observing System Data and Information System Land Processes Distributed Active Archive Center for maintaining and providing the MODIS data; the National Oceanic and Atmospheric Administration/Office of Oceanic and Atmospheric Research/Earth System Research Laboratory Physical Sciences Division for providing NCEP data; the World Climate Research Program’s Working Group on Coupled Modeling, which is responsible for CMIP; the climate modeling groups (listed in SI Appendix, Table S2) for providing their model outputs; and the Numerical Terradynamic Simulation Group at University of Montana and Max Planck Institute for Biogeochemistry for sharing their data. M.K. was supported by National Science Foundation (NSF) Grant EAR-1454983. M.K. and A.P. were supported by NSF Grant EAR-1331846. A.P. and G.G.K. were supported by NSF Grant DGE-1068871. A.P. was additionally supported by NSF Grants EAR-1316258 and FESD EAR-1338694. G.G.K. was additionally supported by NSF Grant EAR-134470 and US Department of Energy Grant DE-SC0011461. C.-W.H. was supported by NSF Grant DEB-1557176.

Footnotes Author contributions: Y.L., A.J.P., M.K., and A.P. designed research; Y.L. performed research; C.-W.H. and G.G.K. further developed model assumptions and scenarios; Y.L. analyzed data; and Y.L., A.J.P., M.K., C.-W.H., G.G.K., and A.P. wrote the paper.

Conflict of interest statement: Ignacio Rodriguez-Iturbe and A.P. were coauthors of a 2015 paper [Feng X, Porporato A, Rodriguez-Iturbe I (2015) Stochastic soil water balance under seasonal climates. Proc R Soc A 471: 20140623].

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1704811114/-/DCSupplemental.