Lateral subsurface flow modulates forest mortality risk to future climate and elevated CO2

Forest mortality has been widely observed across the globe during recent episodes of drought and extreme heat events. But the future of forest mortality remains poorly understood. While the direct effects of future climate and elevated CO2 on forest mortality risk have been studied, the role of lateral subsurface water flow has rarely been considered. Here we demonstrated the fingerprint of lateral flow on the forest mortality risk of a riparian ecosystem using a coupled plant hydraulics-hydrology model prescribed with multiple Earth System Model projections of future hydroclimate. We showed that the anticipated water-saving and drought ameliorating effects of elevated CO2 on mortality risk were largely compromised when lateral hydrological processes were considered. Further, we found lateral flow reduce ecosystem sensitivity to climate variations, by removing soil water excess during wet periods and providing additional water from groundwater storage during dry periods. These findings challenge the prevailing expectation of elevated CO2 to reduce mortality risk and highlight the need to assess the effects of lateral flow exchange more explicitly moving forward with forest mortality projections.


Introduction
Drought and heat-induced forest mortality has been increasingly observed throughout the globe [1]. However, forest mortality risk remains challenging to predict [2], despite its important implications for ecosystem services and sustainability [3,4]. Accurate predictions of forest mortality require not only incorporating the physiological mechanisms underlying plant sensitivity to drought, but also characterizing the various hydrological processes mediating water availability during drought [5]. While hydraulic transport through the soil-plant-atmospheric continuum has been well accepted to mechanistically capture plant sensitivity to environmental cues and droughtinduced mortality [6][7][8][9][10], the role of subsurface flow in modulating plant water availability has only recently received attention [5,[11][12][13]. Existing evaluations of forest mortality risk only considered climatic factors and water exchanges in the vertical dimension, ignoring the potential effects of lateral flows happening belowground. Here we posit the lateral redistribution of groundwater could be as important as climate drivers for mediating mortality risk and thus should be incorporated into future model developments.
Forest response to future climate and CO 2 concentrations has been widely studied. Elevated CO 2 is generally expected to increase photosynthesis and reduce stomatal conductance [14][15][16]. Its impact on soil water dynamics appears much more complex and direct observational evidence remains lacking [17][18][19][20][21]. Modern Earth System Models (ESMs) suggested stomatal closure at elevated CO 2 will reduce evapotranspiration (ET), increase soil moisture, and further increase streamflow [22][23][24][25][26]. Further, recent work of plant hydraulics modeling suggested that elevated CO 2 may largely alleviate mortality risk increase from global warming and increased evaporative demand through stomatal closure and water savings [27][28][29], although counter examples also exist [30]. However, these previous predictions are based on soil moisture schemes that solve water fluxes in the vertical dimension only (figure 1(a)), ignoring water fluxes that move laterally through the landscape.
Landscapes are comprised of soil columns that are interconnected rather than isolated. Various water exchanges in the lateral direction could occur to influence plant water availability and consequently mortality risk (figures 1(b)-(d)). For instance, water excess from CO 2 -induced stomatal closure could be consumed by neighboring plants if there is a soil matric potential gradient (figure 1(b)) [21,31], or by downslope plants due to topography-induced redistribution (figure 1(c)) [12,32,33]. In addition, groundwater dependent ecosystems may rely on water supply from groundwater or streams that have traveled a long distance from places with rather distinctive climate regimes (figure 1(d)) [11,12,34,35]. These lateral hydrological processes might create responses that are decoupled from local climate, or spatially heterogenous even for the same species in the same climate [34,36], both of which are currently missing from mortality risk evaluations. Importantly, lateral flow can increase or decrease local soil moisture depending on gradients of topography and soil matric potential, and can only be quantified by mechanistic models that solve transient water fluxes through the soil-plant-atmosphere continuum vertically and across the landscape laterally [12,19,21,35,37].
To quantify the role of lateral flow in mediating ecosystem responses to projected future climate scenarios, we updated a state-of-the-art, coupled plant hydraulics-distributed hydrology model, Parallel Flow-Terrestrial Regional Ecosystem Exchange Simulator (ParFlow-TREES) [12], with a hydraulic-based stomatal optimization routine [7]. Plant hydraulic transport and the impairment of plant hydraulic systems have been shown to mechanistically represent the physiological stress that trees experience during drought [10,38,39], and provide a useful metric to quantify drought-induced mortality risk [40][41][42]. We forced ParFlow-TREES with multiple ESM projections from the Coupled Model Intercomparison Project Phase 5 (CMIP5) experiments and evaluated forest function and mortality risks at whole ecosystem and across the landscape. To separate the effects of lateral flow, we further performed simulations disabling lateral flow to mimic previous model assumptions that only accounted for water fluxes in the vertical dimension. We focused on riparian ecosystems, because they are simultaneously influenced by various lateral subsurface flow processes (figure 1) and climate [43]. Although focused on one site, the implications of our results are broadly applicable for other riparian as well as upland ecosystems. Specifically, we asked: (a) How does lateral subsurface flow mediate plant responses to future climate factors at whole ecosystem scale? (b) How does lateral subsurface flow create spatially heterogeneous responses within an ecosystem?

Model and evaluation
We used ParFlow-TREES [12] to solve the transient, three-dimensional water fluxes in the subsurface and along the soil-plant-atmospheric continuum. ParFlow-TREES integrates a plant physiology model, TREES [44] to a variably saturated groundwater model, ParFlow [45]. Details on the model were reported previously [12,46]. Briefly, TREES estimates photosynthesis and transpiration at 30 min time steps using the bigleaf approach with sun and shade partition. ParFlow solves the variably saturated flow (i.e. saturated and unsaturated groundwater and surface water) in a single matrix based on hydrodynamic principles, without a priori specification of the types of flow in certain portion of the domain [45,47,48]. The integrated model, ParFlow-TREES, explicitly simulates the landscape represented by grid cells, with each cell covered by one plant species or bare soil and thus effectively captures transient water fluxes through the soil-plantatmosphere continuum vertically and across the landscape laterally. In this study, we further modified ParFlow-TREES [12] by incorporating a stomatal optimization routine based on photosynthetic gain and hydraulic cost [7,41,49]. Plant hydraulic transport and hydraulic-based stomatal optimization have been shown to reliably model whole plant responses to environmental fluctuations of soil water potential, evaporative demand and CO 2 [28,41,[50][51][52], although variable performances have been reported among different hydraulic modeling schemes [53,54]. ParFlow-TREES is very suitable for investigating lateral water exchanges among neighboring soil patches, between upland and lowland, as well as between stream channel and floodplain, based on gradients of soil water, gravity, and subsurface hydraulic properties.
Our study site is a natural riparian cottonwood forest along a 3 km river corridor within the Oldman River valley near Lethbridge, Alberta, Canada (49.702 • N, 112.863 • W). The mean annual temperature was 5.9 • C and average annual precipitation was 380.2 mm [55]. Model domain setup and major parameters (table S1 (available online at stacks.iop.org/ ERL/16/084015/mmedia)) were taken from a previous study [12]. Observations of photosynthetically active radiation, precipitation, temperature, wind speed, vapor pressure deficit (VPD), leaf area index, ET, and net ecosystem exchange (NEE) were made during the growing season of 2014 and 2015 [55]. Gross primary productivity (GPP) was estimated from Eddy covariance measurements of NEE, along with meteorological measurements of light intensity and temperature [55][56][57]. Stream flow data was gathered from the closest gauge station (water ID station 05AD007, data archived by the Water Survey of Environment Canada). These observations were used to drive model simulation and evaluate model predictions of carbon and water fluxes (figure S1). The model performed well at predicting GPP (R 2 = 0.81 in 2014 and R 2 = 0.78 in 2015) and ET (R 2 = 0.82 in 2014 and R 2 = 0.85 in 2015) at daily time steps.

Historic and future hydroclimate
We used a high quality, global dataset of retrospective meteorological forcing [58] and observed streamflow records from the nearest gauge station, for the historic period of 1997-2016. Ambient CO 2 concentration was set at 380 ppm. We used projections of future meteorological and runoff fields from selected ESM predictions under the CMIP5 experiments. We first identified ESMs that (a) reported 3 hourly fields of air temperature, precipitation, humidity, shortwave radiation, air pressure, wind speed, and land surface runoff; and (b) included historic, RCP4.5 and RCP8.5 experiments. These criteria yielded 14 ESMs (table S2). Second, we extracted monthly temperature and precipitation from all 14 ESMs at the grid containing the focal study watershed for the historical period and evaluated their consistency with the historic, meteorological database [58]. For each ESM, we calculated the normalized root mean squared error (RMSE) of precipitation and temperature. We used the sum of normalized RMSE to rank the 14 ESMs (the lower normalized RMSE score means better agreement with historic meteorological database). We chose the four models with highest ranking (table S2) to drive future simulations. We downloaded the 3 hourly meteorological forcing variables for the future period (2081-2100) from the four selected ESMs and two Representative Concentration Pathway (RCP s) (Ca = 531.9 ppm for RCP4.5 and 804.0 ppm for RCP8.5) and linearly interpolated into half-hourly time steps. In most ESMs, runoff is generated as excess of precipitation that does not evaporate and needs horizontal integration along prescribed global river network using river routing model to convert to stream flow [59]. However, this approach was not appropriate to estimate streamflow at our local site. We developed an alternative approach that translates runoff predicted at ESM grid level to local streamflow (figure S2, Methods S1).

Numerical simulations
For both historic and future simulation, uniform forcing was imposed across model domain given its small area. We used the first year of forcing for model spin-up and restarted from the first year to run continuous simulations while re-initializing plant hydraulic states in the beginning of each year. Growing season length was determined based on annual plots of cumulative thermal degree days following a previously described approach [28]. Soil water status was also solved during non-growing season with plant functions turned off. For historic simulation, we used the maximum percentage loss of whole-plant hydraulic conductance (PLKmax) as a metric for mortality risk, adopting a previously published approach [28]. We also calculated the average ET and GPP for each growing season. We used the 20 year average of PLKmax, ET and GPP, either aggregated to the entire ecosystem or spatially distributed across the landscape, to provide the baseline to which future simulations can be compared against. For each ESM and RCP scenario, we followed the same procedure as the historic simulation, and quantified the relative changes in GPP (∆GPP %) and ET (∆ET %) compared to historic average. We quantified changes in mortality risk (∆PLKmax) as the difference between future and historic PLKmax.
To evaluate the potential effects of lateral flow on ecosystem function, we repeated the historic and future simulation with the model parameter of lateral soil permeability set to zero. This manipulation effectively mimicked current ESM assumptions of ignoring lateral flow, and plant water supply is solely driven by vertical water fluxes (climate-only prediction). Further, we performed sensitivity analyses to assess the individual effects of future atmospheric CO 2 , precipitation, and VPD with and without considering lateral flow. We used the historic forcing as the baseline and replaced one variable (CO 2 , precipitation, or VPD) at a time with its future projections by different RCPs and ESMs. Although this approach may neglect the inherent co-variation among hydroclimate variables, it has been a useful approach to disentangle the relative importance of different climate factors on ecosystem functions [60,61]. We quantified the effect size of VPD, precipitation and CO 2 on ecosystem function by calculating the absolute values of ∆GPP%, ∆ET% and ∆PLKmax from each simulation compared to the historic average.

Climate and lateral flow effects on ecosystem responses
At the focal site, a wide range of hydroclimate conditions was projected by selected ESMs for the period 2081-2100 ( figure 2, table S2). Compared to historic conditions, mean annual temperature and precipitation were projected to increase. Future VPD was estimated to increase by most models although ESM inmcm suggested otherwise. Streamflow was predicted to have earlier peaks and lower mean annual flows. Forced by these hydroclimate inputs, predictions of mortality risk, ET and GPP varied across climate conditions projected by different ESMs and RCPs, as expected ( figure S3). Compared to historic conditions, changes in future mortality risk ranged from −30% to 40%, ET changes ranged from −20% to 40%, and GPP was predicted to increase by 10%∼70%. Compared to RCP4.5, the ecosystem was likely to experience slightly higher mortality risk (p < 0.01 for one ESM, and p < 0.1 for two ESMs), no significant differences in ET (p > 0.1), and significantly higher GPP (p < 0.01) under RCP8.5 scenarios.
Predictions of ecosystem functions also varied by whether or not considering lateral flow ( figure 3, open vs solid dots). Averaging across all ESM projections, when climate and lateral flow were used together the decline in PLKmax was negligible (∼1.3%), compared to a decrease of 13% when only climate was considered. Similarly, increases in ET and GPP were also much smaller than the anticipated increase predicted from only considering climate (figures 3(b) and (c)). In particular, climate-only predictions suggested decreasing PLKmax in three out of four ESMs (figure S4), while incorporating lateral flow resulted in little change or a slight increase ( figure S3). This discrepancy can be attributed to the fate of water surplus during periods when water-savings from elevated CO 2 exceeds plant water demand. Without lateral flow, water surplus stays within the local soil column and can be used to sustain plant water uptake during future drought, which is the main mechanism that elevated CO 2 buffers plant mortality risk. However, with lateral flow, a water surplus flowed to neighboring or downstream soil patches instead of staying within the local soil column for future use. Thus, the beneficial effects of elevated CO 2 on buffering mortality risk were largely compromised when lateral flow was accounted for. On the other hand, when CO 2 -induced water savings were not sufficient to counteract climatic drought, as projected by model CNRM (figure 2), predictions with lateral flow suggested lower PLKmax, and greater increases in ET and GPP than climate-only predictions (figure 3). This is because the lateral flow from streams into the floodplain provided additional water supply to buffer plants from drought stress.
The effect size of changing climate on ecosystem functions were much smaller when lateral flow was considered ( figure 4, solid vs open bars). This was consistent across different climate variables and different ecosystem functions. In particular, elevated CO 2 had smaller impacts on ET and PLKmax, compared to its role in boosting GPP ( figure 4(a)). Its impact was further reduced when lateral flow was considered. Finally, lateral flow had stronger influence over PLKmax, compared to its impacts on ET and GPP, suggesting the critical role of lateral flow in estimating mortality risk.

Landscape response to future hydroclimates
Within the same climatic forcing, lateral flow created spatially heterogeneous responses across the landscape, with two sets of emerging spatial patterns (figures 5, S5 and S6). On the one hand, in a wetter future such as the one suggested by inmcm RCP4.5, mortality risk, ET and GPP either monotonically decreased or increased in a direction moving away from the stream channel (figures 5(a)-(c)), or equivalently, with increasing depth to water table ( figure 6). This result appeared in line with previous suggestions that elevated CO 2 would benefit dry sites the most [62][63][64][65]. However, sensitivity analysis suggested the greater benefits towards upland areas (with deep water table) was primarily due to higher precipitation that counteracted topography-driven lateral flow (figures 6(a)-(c), green triangles). CO 2 boosted GPP across the landscape, but minimally influenced ET and PLKmax (figures 6(a)-(c), grey dots). On the other hand, in a drier future predicted by MIROC RCP8.5 (even though with increased precipitation,

Discussion
Improved understanding of the various mechanisms underlying forest mortality is critical for future planning and management. While the impacts of climate factors on mortality risk have been examined, the effects of lateral flow have been neglected. Here, we employed a state-of-the-art, integrated plant hydraulics-hydrology model to quantify the impacts of lateral subsurface flows on ecosystem response to future hydroclimate scenarios. We demonstrated the critical role that lateral flow could have in modulating future forest mortality risk (figure 3), buffering ecosystem sensitivity to individual climatic factors (figure 4), and creating heterogeneous responses across the landscape ( figure 5). Importantly, we found the water-saving and drought  ameliorating effects of elevated CO 2 were much smaller than previously anticipated after accounting for lateral flow ( figure 4(a)). Taken together, our results provided supporting evidence that lateral subsurface flow mediates ecosystem water dynamics and thus should be incorporated into future ESM developments [66] and evaluations of forest mortality risk [28].
Elevated CO 2 strongly influenced the GPP boost (figures 3 and 4), due to the well-known, carbon fertilization effects [14][15][16] in the absence of nutrient limitation [67,68]. But its impact on ET and mortality risk was limited, in contrast to previous studies suggesting elevated CO 2 will reduce ET, increase soil moisture, and alleviate plant water stress [22][23][24][25][26]. This is because water fluxes in the lateral direction removed extra water saved from stomatal closure associated with elevated CO 2 , compromising its potential benefit for local plants. Moreover, concurrent changes in VPD, precipitation and streamflow could overwhelm the water-saving effects of high CO 2 , resulting in increases of both ET and PLKmax ( figure 3). This emphasized that future projections of mortality risk should cover the entire manifold of climate and hydrological conditions as they will likely occur in concert [28]. Taken together, the role of elevated CO 2 in alleviating forest mortality risk is likely to be smaller than previously anticipated from solely from leaf-level stomatal closure.
Lateral flow from streams, or more broadly groundwater storage, provided critical water sources for plants when water from precipitation is not sufficient to meet the water demand [12,13]. But this process has been ignored in previous investigations of climate change droughts. Current ESMs collected water surplus from the bottom of land cells, but do not allow this water to go back to land and support plant use as it travels across the landscape [34,66]. However, the lateral redistribution of surface and subsurface water has been shown to create systematic variations in soil moisture, ET, productivity and drought-induced tree mortality not only along river corridors, but also in upland ecosystems [5,[69][70][71]. Neglecting lateral flow might lead to overestimation of drought stress and mortality risk, such as the climate-only predictions forced by CNRM ( figure 3). While our focal site was likely to receive higher precipitation, globally most land areas are likely to experience a warmer future with either less precipitation or little change in precipitation [72,73]. Thus, subsurface flow will be critical in modulating the future fates of forests that depend on water supply from groundwater in addition to precipitation [74,75]. Further, climate-driven changes in human water demand and associated management practices may intensify the effects of changing stream flow [76]. Moreover, previous studies emphasized the impacts of growing season climate on plant function [27,28], Figure 6. Relative changes in gross primary productivity (∆GPP%; (a), (d)) and evapotranspiration (∆ET%; (b), (e)), and changes in mortality risk (∆PLKmax; (c), (f)) versus averaged historic water table depth (WTD), based on simulations forced by inmcm RCP4.5 (a)-(c) and MIROC RCP8.5 (d)-(f), respectively. Inmcm RCP4.5 projects a future with higher precipitation, temperature, streamflow and lower VPD, while MIROC RCP8.5 projected a future with higher precipitation, temperature, VPD and lower streamflow. Symbol and colors represented different simulations using the full suite of future hydroclimate projections (black dots) or using one future variable projection at a time (colored symbols). but our results indicated non-growing season climate could be equally important through processes that recharge groundwater and contribute to the generation of streamflow, such as snowpack dynamics in snow-dominated environment [77,78].
Lateral subsurface flow created spatial heterogeneities of plant responses within the same ecosystem (figures 5 and S5), which cannot be captured by models that only consider vertical water fluxes [11,74]. Trees distant from streams are minimally influenced by stream flow and their response to future climate were essentially the same as upland trees [76]. Trees growing next to the stream had little change in mortality risk, because they are completely buffered from water stress due to persistent water supply from the stream [12], regardless of changes in VPD and precipitation. Trees in the intermediate zone demonstrated the most variable responses, depending on the interactions among stream flow, climate, tree water uptake, gradients of soil water and gravity. This integrated, dynamic response cannot be predetermined, and can only be captured by models solving transient water fluxes through the soil-plantatmosphere continuum and across the landscape, such as ParFlow-TREES. Further, the more complete physics-based representations of water fluxes is expected to be more robust when applied to novel environmental conditions [79].
As in any modeling study, there were a number of assumptions that warrant further investigations. First, we used the maximum percentage loss of hydraulic conductivity as a proxy of droughtinduced mortality risk. While hydraulic damage has been successful in predicting several episodes of mortality [5,41,46,80], the exact cause and timing of tree death remains complicated to pinpoint and additional processes might be critical [30,[81][82][83][84][85]. Second, plant traits such as rooting depth were assumed uniform across the landscape, whereas plants acclimate and adjust traits to their local environment [86][87][88][89].
This lack of consideration might be responsible for the model estimation of rather high PLKmax (∼60%) for some trees under historic conditions (figure 6(i)), which was likely to be an overestimation. Finally, our current study focused on the riparian corridor and relied on ESM projected runoff to represent future changes in streamflow. Future work should extend to the whole watershed and explicitly investigate the influence of changing climate and management practices on streamflow generation and plant water stress.
While our results should be regarded as potential responses in the absence of these limitations, they illustrated how ecosystem response to elevated CO 2 might deviate from prevailing expectations as lateral subsurface flow is incorporated (figures 3 and 4). The lateral transport of groundwater or stream water could be critical mechanism for maintaining hotspots of ecological refugia to changing climate (figure 5) [90][91][92]. Since the future of groundwater and streamflow is highly uncertain [93] due to both direct and indirect effects of climate change and management decisions, particular focus should be devoted to improve the representation of ecosystem water dynamics when evaluating the response of ecosystem to future climate and elevated CO 2 .

Conclusion
Using a state-of-the-art model that solves transient, three-dimensional water fluxes, we quantified the effects of lateral subsurface flow on ecosystem function and mortality risk to future climate projections. We found lateral subsurface flow could strongly modulate future ecosystem functions and their sensitivity to climatic factors both at the entire ecosystem scale and across the landscape. These findings challenge the prevailing expectation of elevated CO 2 to reduce mortality risk and call for re-evaluation of vegetation persistence in altered climatic and hydrological conditions. Further, our results provide supporting evidences of lateral flow as an imperative to incorporate into ESMs.