Interactive comment on “ Precipitation legacy effects on dryland ecosystem carbon fluxes : direction , magnitude and biogeochemical carryovers ”

Abstract. The precipitation legacy effect, defined as the impact of historical precipitation (PPT) on extant ecosystem dynamics, has been recognized as an important driver in shaping the temporal variability of dryland aboveground net primary production (ANPP) and soil respiration. How the PPT legacy influences whole ecosystem-level carbon (C) fluxes has rarely been quantitatively assessed, particularly at longer temporal scales. We parameterized a process-based ecosystem model to a semiarid savanna ecosystem in the southwestern USA, calibrated and evaluated the model performance based on 7 years of eddy-covariance measurements, and conducted two sets of simulation experiments to assess interdecadal and interannual PPT legacy effects over a 30-year simulation period. The results showed that decreasing the previous period/year PPT (dry legacy) always increased subsequent net ecosystem production (NEP) whereas increasing the previous period/year PPT (wet legacy) decreased NEP. The simulated dry-legacy impacts mostly increased subsequent gross ecosystem production (GEP) and reduced ecosystem respiration (Re), but the wet legacy mostly reduced GEP and increased Re. Although the direction and magnitude of GEP and Re responses to the simulated dry and wet legacies were influenced by both the previous and current PPT conditions, the NEP responses were predominantly determined by the previous PPT characteristics including rainfall amount, seasonality and event size distribution. Larger PPT difference between periods/years resulted in larger legacy impacts, with dry legacies fostering more C sequestration and wet legacies more C release. The carryover of soil N between periods/years was mainly responsible for the GEP responses, while the carryovers of plant biomass, litter and soil organic matter were mainly responsible for the Re responses. These simulation results suggest that previous PPT conditions can exert substantial legacy impacts on current ecosystem C balance, which should be taken into account while assessing the response of dryland ecosystem C dynamics to future PPT regime changes.


Introduction
Drylands play an important role in global carbon (C) cycle and future C sequestration (Houghton et al., 1999;Asner et al., 2003), as they cover 30-45 % of the earth's land surface (Asner et al., 2003;Reynolds et al., 2007), store about 15 % of the global soil organic carbon (Schlesinger, 1991), and represent 30-35 % of terrestrial net primary production (Field et al., 1998).Driven by sporadic precipitation (PPT) and nonlinear biological responses, dryland C fluxes are especially variable across time and space (Maestre et al., 2012;Collins et al., 2014), making the prediction of dryland C budgets a challenging task (Jenerette et al., 2012).Moreover, climate models predict that the intra-and interannual PPT variability may be further intensified in dryland regions with longer drought durations and more large-sized events (Solomon et al., 2007;Diffenbaugh et al., 2008;Cook and W. Shen et al.: Precipitation legacies on dryland C fluxes Seager, 2013).Further, sequences of wet years followed by sequences of dry years and vice versa are also increasingly likely (Peters et al., 2012;Sala et al., 2012).Understanding the response of dryland ecosystem C fluxes to PPT variation is, therefore, important to characterizing the global C cycle and predicting how future PPT regime changes will affect dryland C balance.
As a measure of ecosystem C balance, net ecosystem production (NEP) has a value that is positive when an ecosystem accumulates C and negative when an ecosystem loses C. Dryland NEP is closely tied to current-year PPT amount, with wetter-than-average years being a C sink, drier-than-average years being a C source, and years with average rainfall being C neutral (Flanagan et al., 2002;Hastings et al., 2005).Additionally, at seasonal scales, the distribution of PPT in addition to the total amount can have large influences on ecosystem production (Porporato et al., 2004;Katul et al., 2007).At interannual scales a PPT legacy effect, defined as the impact of past PPT conditions on the current structure and functioning of ecosystems (Lauenroth and Sala, 1992;Sala et al., 2012;Monger et al., 2015), has also been found to play an important role in shaping the temporal variability of dryland ecosystem C fluxes (Knapp et al., 2002;Huxman et al., 2004a, b;Heisler and Weltzin, 2006;Sala et al., 2012;Ogle et al., 2015).For example, Hasting et al. (2005) attributed the C sink status of a desert shrub ecosystem in the early spring of 2002 to the above-average rainfall in the late fall of 2001.Scott et al. (2009) and Hamerlynck et al. (2013) found that a cool-season (December-April) drought was followed by an unusually large net C loss during the following warm monsoon season (July-September) in a semiarid savanna and a semi-desert grassland.Moreover, the savanna ecosystem has recently been a net C source, and one hypothesized but untested explanation is due to an increase in current respiration of organic C that accumulated in the preceding wetter decade (Scott et al., 2009).While these studies reveal the existence of PPT legacy effects on NEP at the seasonal scale, only a few studies have quantitatively assessed the contribution of PPT legacy to the temporal variability of dryland NEP at interannual and interdecadal timescales (Williams and Albertson, 2006), mainly because it is methodologically difficult to separate the past and current PPT impacts on C fluxes with the limited observational data (Sala et al., 2012), and there is a general lack of field manipulative experiments to address the PPT legacies at these scales (Reichmann et al., 2013a).
Much of our current understanding of the PPT legacy effects on dryland C fluxes is based on aboveground net primary production (ANPP).A number of studies have documented that dryland ANPP is not only linearly related to current-year PPT but also closely related to the PPT amount and seasonality several months to years before (Lauenroth and Sala, 1992;Oesterheld et al., 2001;Huxman et al., 2004c).For example, field studies have found a positive legacy impact where ANPP is higher than expected if pre-ceded by a wetter year, or lower than expected if preceded by a drier year (Jobbagy and Sala, 2000;Oesterheld et al., 2001;Wiegand et al., 2004;Sherry et al., 2008;Sala et al., 2012).Proposed mechanisms explaining such observed positive PPT legacy effects on ANPP mainly involve the structural carryovers between years, which can be leaf and root biomass (Oesterheld et al., 2001); the composition of species differing in rooting depth and phenology (Paruelo et al., 1999;Jobbagy and Sala, 2000); or the density of seeds, tillers and plant individuals (Oesterheld et al., 2001;Yahdjian and Sala, 2006;Reichmann et al., 2013a).Alternatively, a negative legacy effect occurs when production is lower than expected if preceded by a wet period or higher than expected if preceded by a dry period (Jenerette et al., 2010).A negative PPT legacy effects may be influenced more by biogeochemical carryovers that influence the resource availability to respond to current PPT (Evans and Burke, 2013;Reichmann et al., 2013b), whereby increased growth in response to a higher PPT can reduce the available nutrients (e.g., nitrogen, N) for the following period and vice versa.Although various mechanisms have been proposed for the PPT legacy impacts on ANPP, few of them have been rigorously tested, and the key underlying mechanisms still remain poorly understood (Sherry et al., 2008;Williams et al., 2009;Sala et al., 2012;Monger et al., 2015).
Soil respiration (R s ), as a major component of ecosystem C efflux, has also been found to have lagged responses to PPT variations (Huxman et al., 2004b;Sponseller, 2007;Ma et al., 2012;Cable et al., 2013).This is particularly true at the event scale; after a period of drought, a rainfall event can result in a pulse of CO 2 efflux that may be orders of magnitude larger than that before the event and then decline exponentially for a few days to weeks (Xu et al., 2004;Jenerette et al., 2008;Borken and Matzner, 2009;Cable et al., 2013;Oikawa et al., 2014).At a seasonal scale, Vargas et al. (2010) found no lags between R s and soil moisture across 13 vegetation types, including four grasslands; however, Hamerlynck et al. (2013) presented longer-term ecosystem flux data that suggest seasonal drought legacies affect ecosystem respiration (R e ) in a semi-desert grassland in southeastern AZ, USA.They posited that the increased C substrate availability resulting from the previous cool-season drought-induced plant mortality was responsible for the higher R e in the following monsoon season.However, very few studies have been devoted to understanding the PPT legacy impacts on dryland respiration at greater than seasonal timescales.
In this study, we conducted simulation experiments with a widely used dryland ecosystem model, Patch Arid Land Simulator (PALS; Kemp et al., 1997Kemp et al., , 2003;;Reynolds et al., 2004;Shen et al., 2009), to analyze the PPT legacy effects on ecosystem-level C fluxes including NEP, gross ecosystem production (GEP), and R e .The PALS model was built on the pulse-reserve concept (Noy-Meir, 1973) and had been used to analyze the impacts of antecedent moisture conditions and the lagged responses of different plant functional types (FTs) in three North American deserts at the rainfall event scale (Reynolds et al., 2004).We parameterized, calibrated, and evaluated the model based on the long-term eddycovariance-measured fluxes in a semi-desert savanna ecosystem in the southwestern USA (Scott et al., 2009) to analyze the PPT legacy effects at interannual and interdecadal scales.Specifically, we addressed the following two questions.First, what are the direction and magnitude of ecosystem C flux responses to dry and wet legacies?We expected that the PPT legacy impacts would occur over annual and decadal scales in correspondence to PPT fluctuations at these scales, and the dry-and wet-legacy impacts would differ in direction and magnitude.Second, how are the direction and magnitude of PPT legacy effects related to the PPT characteristics of both the previous and the current year/period?We expected that greater variability in PPT would lead to corresponding increases in legacy effect.For PPT characteristics, we were interested not only in the annual and seasonal PPT amount but also in between-event interval and event size distribution, since these variables are widely recognized key PPT features to dryland ecosystems (Porporato et al., 2004;Katul et al., 2007;Shen et al., 2008a).

Model description
PALS is a process-based ecosystem model that consists of four modules: atmospheric forcing, a water cycling and energy budget, plant production and respiration, and soil organic matter (SOM) decomposition and heterotrophic respiration (R h ).The four modules are interactively linked by the cycling of C, N, and H 2 O through the atmosphere-plantsoil continuum.The PALS model explicitly considers seven plant FTs commonly found in the North American warm deserts: evergreen shrub, deciduous shrub, perennial forb, perennial C 3 and C 4 grasses, and native and exotic C 3 annual grasses (Reynolds et al., 1997;Shen et al., 2009).Since the detailed model structure and mechanistic relationships have been presented in several publications (Kemp et al., 1997(Kemp et al., , 2003;;Reynolds et al., 1997Reynolds et al., , 2000Reynolds et al., , 2004;;Gao and Reynolds, 2003;Shen et al., 2005Shen et al., , 2008aShen et al., , b, 2009)), here we briefly describe the four modules and refer to the specific literature for detailed description.
The atmospheric driving force module reads in data for atmospheric driving variables (e.g., atmospheric [CO 2 ], N deposition rate, daily maximum and minimum air temperatures, PPT, relative humidity, and solar radiation) and, based on these driving variables, calculates other important variables such as vapor pressure deficit (VPD), which directly influences stomatal conductance and indirectly influences soil temperature, SOM decomposition and soil respiration.Calculations of VPD and soil temperature can be found in Eqs. ( 2)-( 7) in Shen et al. (2005).
The water cycling and energy budget module mainly calculates soil water contents at six layers, the rates of water infiltration into and percolation out of a layer, and water losses via evaporation and transpiration from different layers.Water infiltration and percolation rates of a layer are determined by the effective PPT reaching the soil surface, previous water content, and the water holding capacity as a function of soil texture (Shen et al., 2005).Soil evaporation is determined by soil water availability and energy available in the two top soil layers (10 cm in depth).Water uptake by plants is partitioned among the soil layers according to the proportion of roots in each layer for all plant FTs (Kemp et al., 1997;Shen et al., 2008b).Canopy transpiration is calculated by using the energy budget and the canopy stomatal resistance (Reynolds et al., 2000;Gao and Reynolds, 2003).
The plant production and respiration module mainly simulates phenology, primary production, growth and maintenance respiration, photosynthate allocation, and litterfall of each plant FT.Three major phenophases (i.e., dates of germination, leafing, and dormancy) are determined in PALS based on the observed dates, air temperature, and PPT (Shen et al., 2009).Primary production for each FT is calculated based on the leaf area, potential net photosynthetic rate, stomatal conductance, leaf N content modifier, and the difference between intercellular and atmospheric [CO 2 ].The plant photosynthesis rate is estimated as a product of stomatal conductance and the partial pressure gradient between atmospheric and intercellular [CO 2 ].The stomatal conductance is calculated as an exponential function of leaf water potential that decreases linearly with atmospheric vapor deficit (see Eqs. ( 10)-( 14) in Shen et al., 2005).Photosynthate is allocated to different plant organs (leaf, stem, and root) using fixed allocation ratios after subtracting the maintenance respiration, which is estimated as a function of live biomass, basal respiration rate, and modifiers of temperature and plant water potential (Shen et al., 2008a).Growth respiration is calculated based on the growth yield coefficient and the net photosynthate used for growth (Shen et al., 2008a).Litterfall amount is mainly determined as a function of observed dormancy dates, maximum air temperature and drought conditions (Shen et al., 2008a(Shen et al., , 2009)).
The SOM decomposition and heterotrophic respiration module simulates the decomposition of metabolic and structural litter material; SOM in active, slow and passive pools; and CO 2 emissions associated with these decomposition processes (Kemp et al., 2003;Shen et al., 2009).The SOM decomposition rate or heterotrophic rate is calculated as a firstorder kinetic rate with a decomposition coefficient multiplied by the pool size and the temperature and moisture scalars (see Eqs. (A4)-(A11) in Shen et al., 2009).In addition, this module also simulates the dynamics of soil mineral N pool by using N mineralization and atmospheric deposition as the major inputs, and plant N uptake and leaching loss as the major outputs.Among these the N mineralization and plant uptake processes are modeled in more detail while the rates of the other processes are basically assigned with empirical constant values.The N mineralization processes are directly coupled to litter and SOM decomposition processes and are calculated as a product of the C flow rates and the C / N ratio of the corresponding litter or SOM pools (Parton et al., 1993;Kemp et al., 2003).The plant N uptake is a product of water transpiration and N concentration in soil solution (see Eq. 8 in Shen et al., 2008b).

Model parameterization
For this study, we modified and parameterized PALS to represent an upland mesquite savanna ecosystem in the Santa Rita Experimental Range (SRER; 31.8214• N, 110.8661 • W, elevation 1116 m), about 45 km south of Tucson, AZ, USA.Soils at this site are a deep sandy loam (Scott et al., 2009), and the mean groundwater depth likely exceeds 100 m ( Barron-Gafford et al., 2013).PPT was therefore considered as the only source of water input into the system.Based on the vegetation composition (Scott et al., 2009), there were five major plant FTs included in PALS: shrub (e.g., Prosopis velutina), subshrub (e.g., Isocoma tenuisecta), C 4 perennial grass (e.g., Digitaria californica), perennial forb (e.g., Ambrosia psilostachya), and C 3 annual grass, among which the velvet mesquite shrub with average height of ca.2.5 m accounted for ∼ 35 % of the total canopy cover and other FTs (mainly perennial grasses) accounted for ∼ 22 % (Scott et al., 2009).Therefore, we derived the site-characteristic parameters for the two major FTs (shrub and perennial grass) from previous studies carried out in SRER, with those for the other FTs being adopted from a generic parameter data set for the PALS model to be used in the North American warm deserts (Reynolds et al., 2004;Shen et al., 2005).These site-specific parameters mainly included plant-related parameters (e.g., canopy cover, C allocation ratio, rooting distribution ratio, and the initial values of living and dead plant biomass pools) and soil-related parameters (e.g., soil chemical and physical properties, C / N ratios, decomposition rates, and initial values of the litter and SOM pools).The values of these parameters are provided in Table S1 in the Supplement, with the cited literature also being listed below the table.
For the climatic variables used to drive the PALS model, we compiled a 30-year meteorological data set that included daily PPT, maximum and minimum air temperatures (T max and T min ), relative humidity (RH), and total solar radiation (S rad ) from 1981 to 2010.The T max , T min , RH, and S rad data from 1981 to 1990 were observations from the Tucson weather station (about 50 km north of the mesquite savanna site and lower elevation) and obtained by accessing Arizona Meteorological Network online data (AZMET: http:// ag.arizona.edu/azmet).The remaining 20 years (1991-2010) of T max , T min , RH and S rad data were observations from the Kendall Grassland meteorological site (about 85 km east of the mesquite savanna site and slightly higher elevation) and obtained by accessing Southwest Watershed Research Cen- ter (SWRC) online data (http://www.tucson.ars.ag.gov/dap/).The 30-year PPT data were observations from the Santa Rita watershed rain gage no. 5 (1.5 km from the site) and obtained also from the SWRC online data access.These different sources of meteorological data were adjusted based on the 7 years (2004)(2005)(2006)(2007)(2008)(2009)(2010) of meteorological data obtained from the AmeriFlux eddy-covariance flux tower at the mesquite savanna site (US-SRM; see Fig. S1 in the Supplement).Lastly, we used the AZMET and SWRC data from 1981 to 2003 plus the flux tower data from 2004 to 2010 to drive the model.
Since our simulation experiment was based on the manipulations of the 30-year (1981-2010) PPT data, we report the PPT characteristics here in more detail.In the past 30 years, the mean annual PPT (MAP) amount was 401 mm at the site, slightly greater than the long-term  mean of 377 mm (Scott et al., 2009).These 30 years were divided into two periods: a wet period from 1981 to 1994 with a MAP of 449 mm and a dry period from 1995 to 2010 with a MAP of 347 mm (Fig. 1a).For the analysis of PPT legacy effects at interdecadal scale, the wet period was treated as the previ-ous period and the dry period as the current period.For the analysis of PPT legacy effects at interannual scale, the annual scale was defined as being from July through June of the next year.To analyze the relationship between PPT legacy effects and seasonal rainfall characteristics, each year was further divided into four seasons (with their mean rainfall in parentheses): the main warm growing season from July to September (warm-GS, 224 mm), the cool dry season from October to November (cool-DS, 48 mm), the minor cool growing season from December to March (cool-GS, 104 mm), and the warm dry season from April to June (warm-DS, 26 mm).At the site, as in many other dryland regions (Sala et al., 1992;Heisler-White et al., 2008), most rainy days have only lightrainfall events.About 80 % of daily rainfall was < 10 mm, with medium-to large-sized events (10-50 mm) accounting for about 20 % and only 10 events larger than 50 mm in the 30 years (Fig. 1b).The no-rain-day duration between events (hereafter between-event interval or BEI) was ∼ 5 days on average in the warm-GS and ∼ 10 days in the cool-GS (Fig. 1c).The average BEI was ∼ 17 days in the cool-DS and 24 days in the warm-DS, but there could be no rain for 3 months in these dry seasons (Fig. 1c).

Model calibration and evaluation
After model parameterization, we calibrated the model based on 4 years (2004)(2005)(2006)(2007) of CO 2 and H 2 O flux data monitored using the eddy-covariance technique at the savanna site.Detailed descriptions of instrumentation, sensor heights and orientations, and data-processing procedures for the eddycovariance data can be found in Scott et al. (2009).During model calibration, we mainly adjusted the parameter values of photosynthate allocation ratios, live biomass death rates, and SOM decomposition rates to achieve a best fit between modeled and observed GEP and R e , since these parameters have been identified as the most sensitive and uncertain ones (e.g., photosynthate allocation ratios) in influencing the modeled ecosystem carbon fluxes (Shen et al., 2005).The model performed well in capturing the seasonal variation patterns of actual evapotranspiration (AET), GEP, R e , and NEP in the 4 calibration years (Fig. S2), with larger C fluxes during the warm-GS than in the other seasons.At the annual scale, simulated AET, GEP, and R e explained over 60 % of the variations in the observations (Fig. 2, left panels), but the correlation between the simulated and observed NEP was very weak (Fig. 2d).This was mainly because the model substantially overestimated GEP (120 g C m −2 simulated vs. 52 g C m −2 observed) in the cool-GS of 2006 (Fig. S3b).Further explanations on the possible causes of the GEP overestimation in 2006 shall be provided later in the Discussion section.If the data of this year were excluded, the explanatory power for annual NEP could reach 74 %.Since our goal was to use an empirically plausible model to understand the long-term temporal variations in ecosystem fluxes, we consider the calibration results acceptable.served and simulated fluxes with the corresponding currentyear PPT to see how the flux variations were explained by current-year PPT under baseline conditions (i.e., the PPT variations shown in Fig. 1).The explanatory power (R 2 ) for both the observed and simulated fluxes were mostly over 70 % (Fig. 2, right panels), which further indicates that the model is capable of capturing the impacts of PPT variability on ecosystem fluxes.The following simulation experiments were therefore designed to discriminate the contributions by previous-and current-year PPT impacts.

Simulation experiments
We designed two sets of simulation experiments to examine the interdecadal and interannual PPT legacy effects.To analyze the interdecadal legacy effects, we first changed the PPT of the 14-year previous period (1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)) by 0, ±10, ±30, ±50 and ±80 % (multipliers of existing daily PPT amounts in the record) while keeping the 16-year currentperiod (1995-2010) PPT unchanged.After these manipulations, the average PPT of the previous period ranged from 93 mm, corresponding to the 80 % decrease, to 837 mm, corresponding to the 80 % increase.This design detects how changes in previous-period PPT influence the current-period C fluxes and the associated C pool dynamics.On top of each previous period PPT manipulation level, we further changed the current-period PPT by 0, ±10, ±30, ±50, and ±80 %, which resulted in the average current-period PPT varying from 69 to 621 mm.This design detects how changes in the current-period PPT influence the legacies resulting from changes in the previous-period PPT.As a result, we conducted 73 simulation runs, corresponding to the 73 combinations of the above previous-and current-period PPT manipulations (9 previous PPT levels times 8 current PPT levels plus 1 baseline run).
To analyze the interannual legacy, we changed the PPT of each individual year by ±30 % while keeping the PPT of the subsequent years unchanged.This design resulted in 54 simulation runs (27 years from 1981 to 2007 times 2 PPT manipulation levels) and illustrates the effects of changes in the PPT of the previous 1 year on the C fluxes and resource pools of the current year(s).After a 30 % PPT change, annual PPT ranged from 162 to 925 mm in the 27 years, which was large enough to cover the PPT interannual variation at the study site.Another consideration of using 30 % as the PPT manipulation level was that future projected annual PPT variation in dryland regions will be −30 to +25 % (Bates et al., 2008;Maestre et al., 2012).

Data analysis
The legacy effect was quantified as the C flux (or resource pool size) of the current period/year after PPT changes in the previous period/year minus that without PPT changes in the previous period/year.As an example, the following equation calculates the legacy effect of increasing the previous-period PPT by 30 % on the current-period NEP: where NEP CP PPT+30 % is the cumulative NEP throughout the current period (1995-2010) under a 30 % previous-period (1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994) PPT increase; NEP CP PPT+0 % is the cumulative NEP throughout the current period with no previous-period PPT change (i.e., the baseline PPT conditions shown in Fig. 1).This method directly quantifies whether changes in PPT of the previous period will impose a positive, a negative, or no legacy effect on the C fluxes (or resource pools) of the current period.For simplicity, hereafter we refer to the legacy effect resulting from the decreased previousperiod/year PPT as the dry legacy and that resulting from the increased previous-period/year PPT as the wet legacy.Spearman correlation analysis was used to detect the relationships between legacy effects and PPT characteristics, including PPT amount, BEI, and the number of large (≥ 10 mm) vs. small (< 10 mm) events at yearly and seasonal scales.The correlation analysis was performed in SPSS 16.0 (Chicago, IL, USA).
The simulated absolute magnitude of the PPT legacy influence on ecosystem C fluxes (i.e., GEP, R e , and NEP) generally increased with the absolute magnitude of changes in the previous-period PPT (Figs. 3,4).Increasing the currentperiod PPT generally amplified the legacy effects compared to decreasing the current-period PPT (comparing the left to the right panels of Fig. 3).The magnitude of the PPT legacies was also significantly correlated with the PPT difference between the current and previous period ( PPT, equals to the current-period PPT minus the previous-period PPT; Fig. 4).If the previous period was wetter than the current period (i.e., PPT < 0 or a wet-to-dry period transition), the legacy effect on R e was negatively correlated with PPT (Fig. 4c) but that on NEP was positively correlated with PPT (Fig. 4e), indicating more current-period C release after a wetter previous period.In contrast, if the previous period was drier than the current period (i.e., PPT >0 or a dry-to-wet period transition), the correlations were all positive for GEP, R e and NEP (Fig. 4, right panels), indicating more current-period C sequestration after a drier previous period.
The resource pool dynamics were also shaped by the alterations in the previous-and current-period PPTs.We only showed the 30 % decrease and increase in the previous-and current-period PPT (i.e., 4 out of 72 pairs of PPT change combinations) as representative examples in Fig. 5, because the major response patterns for the other paired combinations were similar.The PPT legacy impacts generally lasted for about 6-8 years for plant biomass, litter mass and soil water content (SWC), and much longer for soil organic matter (SOM) and soil mineral N (N soil ; Fig. 5).Based on the resource pool responses in the early 1-2 years (i.e., 1995 and 1996) of the current period, the dry legacies decreased biomass, litter and SOM (Fig. 5a-f), but positively impacted N soil (Fig. 5g-h).Contrastingly, the wet legacies increased biomass, litter and SOM (Fig. 5a-f) but negatively impacted  (1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994).Interdecadal legacy effects on carbon fluxes (e.g., NEP) are calculated as the difference between the current-period flux with previous-period PPT changes and that without previousperiod PPT changes.Sample size is 41 for the wet-to-dry period transition (left panels) and 23 for the dry-to-wet period transition (right panels).GEP represents gross ecosystem production, R e ecosystem respiration, and NEP net ecosystem production.R 2 is the coefficient of determination, and P is probability.
N soil (Fig. 5g-h).Similar to the influences on C fluxes, increasing the current-period PPT (Fig. 5, right panels) amplified the PPT legacy impacts on biomass and litter (Fig. 5ad), and hastened the recovery rates of SOM and N soil to their baseline levels (Fig. 5e-h).

Interannual legacy
At the interannual scale, a 30 % decrease or increase in PPT could have legacy impacts on ecosystem C cycling lasting for 2-12 years (Fig. 6a-b).Notably, the direction of GEP and R e responses to decreasing or increasing previous-year PPT could be positive or negative (Fig. 6c-f).The dry-or wet-legacy effects on these two fluxes were variable; idiosyncratic; and, in some cases, large at this timescale.However, the simulated dry legacies mostly increased NEP (Fig. 6g), whereas the simulated wet legacies mostly decreased NEP (Fig. 6h), which was similar to legacy responses at the interdecadal scale (Fig. 3e-f).The correlation analysis showed that not only rainfall amount but also BEI and event size distribution could influence the magnitude of the simulated dry and wet legacies (Table 1).The warm-GS PPT of a previous year was significantly correlated with the dry legacies for NEP and the wet legacies for GEP and NEP (Table 1).On the other hand, the cool-GS PPT of a current year influenced the dry and wet legacies for C fluxes, but not all of them were statistically significant (Table 1).These results indicate that the legacies were mainly generated in the warm-GS of a previous year, but the current-year cool-GS PPT conditions could influence the C flux responses to the previous-year legacies.Unlike at the interdecadal scale (Fig. 4), our correlation analysis showed that only the dry legacies for NEP had significant correlations with the PPT difference ( PPT) between 2 consecutive years or cool-GSs (Table 1), indicating that the larger the PPT difference between a previous dry year and a current wet year, the greater the legacy impacts on NEP imposed by the previous dry year.
To analyze the interannual PPT legacy impacts on the dynamics of resource pools (i.e., biomass, litter, SOM, N soil , and SWC), 2 wet years (1983 and 1994) and 2 dry years (1986 and 1995) were chosen as examples (see Fig. 1a).The simulated dry legacies reduced biomass, litter and SOM but increased N soil and SWC in the first current year (Fig. 7).In contrast, the simulated wet legacies imposed the opposite direction of impacts on the five resource pools (Fig. 7).The simulated PPT legacy impacts on the resource pools could also last for several years, and the direction and magnitude of the legacy impacts in the following years could differ from those in the first year as described above.For example, increasing the PPT of 1995 by 30 % caused a positive legacy impact on the biomass of the first following year (i.e., 1996; Fig. 7b), but it became negative in the later following years (e.g., in 1998; Fig. 7b), further indicating that current-year PPT conditions can influence the direction and magnitude of previous-year PPT legacies.

Direction and magnitude of the simulated PPT legacies
Through this simulation analysis we demonstrated that previous PPT could impose substantial legacy impacts on current ecosystem C fluxes at interannual and interdecadal timescales.Notably, our simulation results support the hypothesis proposed for our study site (Scott et al., 2009) that the accumulated SOM during the previous wet period contributed to the net C release from the ecosystem during the current dry period.This specific test illustrates a major finding from our simulation study of a negative PPT legacy effect on NEP; i.e., decreasing previous PPT increased current NEP, whereas increasing previous PPT decreased current NEP (Figs. 3,6).Increasing prior PPT (wet legacy) led to limited changes in GEP but consistently increased R e .Decreasing prior PPT (dry legacy) led to more variable effects for both GEP and R e that were strongly conditioned on current-period PPT such that increasing current PPT was associated with increases in the dry-legacy effect.Overall, the effects on GEP were larger than R e for reduced prior PPT and smaller for increased prior PPT, which resulted in a consistent negative PPT legacy on NEP regardless of current PPT.
The complexity in the legacy effects on ecosystem C cycling we show here are in part influenced by the contrasting PPT legacy responses of C uptake and emission and their distinct interactions with current PPT distributions.
In projecting future dryland C dynamics, the effects of PPT legacies increase the complexity of ecosystem responses to PPT variability.One consistent interaction between legacy and current PPT effects was that larger between-period PPT differences could result in larger legacy effects (Fig. 4), which is in agreement with what has been found in some field studies.For example, the magnitude of drought legacy on ANPP is proportional to the severity of the drought (Yahdjian and Sala, 2006;Swemmer et al., 2007), and dry-or wet-year legacies on ANPP are linearly related to the PPT difference between years (Sala et al., 2012;Reichmann et al., 2013a).Our simulation analysis detected that not only annual PPT amount but also finer-scale PPT characteristics such as GS rainfall, BEI, and event size could be important in determining the interannual-scale PPT legacy effects (Table 1).These simulation results suggest that PPT legacy effects may play a more important role in shaping the temporal variability of dryland ecosystem C fluxes under the projected increase in future PPT variability (Solomon et al., 2007;Cook and Seager, 2013) but that their characterization remains a challenge.
The influence of PPT legacies on dryland ecosystem C balance may strongly interact with other sources of variability in dryland C balance, including current-year PPT (Flanagan et al., 2002;Hastings et al., 2005), growing-season length (Xu and Baldocchi, 2004;Ma et al., 2007), seasonal drought (Scott et al., 2009(Scott et al., , 2010;;Hamerlynck et al., 2013), and other factors such as temperature and vegetation composition (Hui et al., 2003;Hamerlynck et al., 2010;Barron-Gafford et al., 2012;Scott et al., 2014).These interactions are shown by several examples from our simulations.While PPT was wetter than normal in 1987 (537 mm), the NEP was −85 g C m −2 yr −1 (a C source), due to the negative wetlegacy impacts on NEP from several previous wet years before (1982-1985; see Fig. 6h).PPT was nearly normal in 2008 (402 mm), but the simulated NEP was 80 g C m −2 yr −1 (a C sink), due to the positive dry-legacy impacts on NEP from several previous dry years (2002)(2003)(2004)(2005)(2006)(2007); see Fig. 6g).Our findings of substantial PPT legacy effects are consis-tent with a recent analysis of 14 years (1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011) of eddycovariance measurements, where Zielis et al. (2014) reported that inclusion of previous year's weather (PPT and temperature) into the linear predicting models for NEP increased the explained variance to 53 % compared to 20 % without accounting for previous year's weather, indicating that previous year's weather also played an important role in determining the C balance of the subalpine spruce forest in Switzerland.Although response patterns generated from this simulation study compared well with previous field observations, there exists no field study that, to our knowledge, provides a similarly comprehensive analysis of PPT legacies.The simulation experimental design of this study provides helpful insights into designing field manipulative experiments to further test the modeled patterns by focusing on contrasting wet and dry legacies, separating ecosystem production and decomposition, and exploring the difference in prior and current PPT on the magnitude of the PPT legacy effect.

Potential mechanisms of the modeled PPT legacies
There are three basic mechanisms explaining why PPT legacy impacts can occur in a model system like PALS.First, the rate of C fluxes is a function of not only various environmental factors (e.g., PPT and temperature) but also the pool size itself.For example, soil heterotrophic CO 2 efflux (R h ) rate is a product of the decomposition coefficient, two scalar functions accounting for temperature and moisture influences, and also the size of the SOM pool (Kemp et al., 2003;Shen et al., 2009).Change in the SOM pool size from previous PPT thereby affects current R h .Second, different C pools have different turnover rates that determine whether biogeochemical materials (e.g., biomass or SOM) can be carried over.If the material produced in a previous year has a turnover rate of less than 1 year, it will not be carried over to the next year to form a legacy impact as explained in the first mechanism.In addition, the turnover rates of different C pools also determine legacy duration.For example, SOM pools in the model have relatively slower turnover rates than biomass pools (Shen et al., 2005(Shen et al., , 2008b)), thus resulting in the longer-lasting legacy impacts on SOM than on biomass or litter pools (Figs. 5 and 7).Third, the interactions between C fluxes and resource (e.g., N or water) availability also determine the direction and magnitude of legacy effects.For example, N carryover as a legacy of a prior dry period (Fig. 5g,  h) can impose impacts on the current-period GEP only when the current-period PPT is not so limiting (Fig. 3b).These are the general mechanisms explaining the occurrence of the modeled PPT legacies from a systems perspective.Below we discuss more specifically the major patterns and the responsible biogeochemical carryovers found in this study.
An intuitive first explanation of the simulated wet legacies would be the carryover of water.However, in most cases soil water carryover did not occur because the wet legacies on SWC were mostly negative or close to zero at the beginning of the current period/year (Figs.5i-j; 7i-j).Soil water carryover was therefore not the major contributor to the modeled PPT legacy effects at interdecadal and interannual scales.This simulation result corroborates field studies that have shown that carryover of water across long temporal scales is rare in dryland ecosystems, because the rainy growing seasons or wet years are often separated by dry dormant seasons or dry years resulting in short residence times (Oesterheld et al., 2001;Reichmann et al., 2013a;Scott et al., 2014).
The carryover of soil N (N soil ) is mainly responsible for the modeled GEP responses.In the PALS model, the photosynthetic rate is linearly related to N availability if plant N demand is not fulfilled (Reynolds et al., 2004;Shen et al., 2005).Therefore, the enhanced N soil from dry legacies (Figs.5g, h and 7g, h) generated mostly positive responses of GEP (Figs. 3a, b and 6c).The simulated dry legacies increased N soil mainly through suppressed plant growth (e.g., the reduced biomass and litter production shown in Figs. 5  and 7) that limited N uptake, which is consistent with the results of many field measurements that N soil accumulates under drought conditions (Reynolds et al., 1999;Yahdjian et al., 2006;Yahdjian and Sala, 2010;de Vries et al., 2012;Evans and Burke, 2013;Reichmann et al., 2013b).Although diverse mechanisms of inorganic N accumulation during dry periods have been proposed in field studies -such as the diffusion restriction of N ions in thin water films of dry soil, the reduced N immobilization by microbial growth and plant uptake, and the reduced N loss from the soil via leaching (Yahdjian et al., 2006) -our simulation results suggest that reduced plant uptake may be the main contributor to the N soil accumulation during dry periods.Given the accumulated N soil as a dry legacy, how ecosystem C fluxes such as GEP respond to this dry legacy may be influenced by current PPT conditions.When current PPT conditions were favorable (e.g., the increasing current-period PPT treatment shown in Fig. 3b and the relatively wet years shown in Fig. 6c), GEP mostly increased with a dry legacy (or the accumulated N) because both N and H 2 O availabilities were favorable for plant growth (or GEP).Contrastingly, when current PPT conditions were unfavorable (e.g., the decreasing current-period PPT treatment shown in Fig. 3a and the relatively dry years shown in Fig. 6c), the GEP responses could be reduced because of the constrained plant growth and the reduced biomass in previous dry years (see Figs. 5c and 7b).
Similarly, the mostly negative responses of GEP to wet legacies (see Figs. 3a, b and 6d) can be explained by the reduced N soil (Figs. 5g,h and 7g,h).The decrease of N soil with increasing PPT in the PALS model is mainly attributed to the increases in plant N uptake and the N leaching loss that is calculated as a linear function of PPT amount (Shen et al., 2005).Similar to our simulation results, several field studies found that N uptake increases and N soil decreases under wet conditions in dryland ecosystems (McCulley et al., 2009;McCalley and Sparks, 2009;Yahdjian and Sala, 2010;Reichmann et al., 2013b).However, contrary to our model assump-tion that N leaching loss is greater in wet than in dry years, some recent field studies have reported that the N leaching loss actually is higher in dry than in wet years or at wet sites (McCulley et al., 2009;Evans et al., 2013;Reichmann et al., 2013b;Homyak et al., 2014), resulting in a more "open" N cycle under drier conditions.If these recent field study results are also true for our semi-desert savanna ecosystem, the model assumption could potentially cause an overestimation of N soil carryover effects as shown in Figs. 3 and 6.Further studies are needed to discriminate the relative contributions of different N processes (e.g., plant uptake, microbial immobilization and mineralization, denitrification, ammonia volatilization, and leaching) to the dynamics of soil inorganic N pools.Nevertheless, this simulation analysis highlights the importance of interactions between N and H 2 O availabilities in creating the legacy impacts of PPT and in shaping the temporal variability of dryland ecosystem C fluxes.
The carryover of organic material (biomass, litter and SOM) is mainly responsible for the modeled R e responses.In the PALS model, the autotrophic (R a ) and heterotrophic (R h ) respiration rates are linearly related to the size of biomass, litter and SOM pools (Kemp et al., 2003;Shen et al., 2008aShen et al., , 2009)).The previous wet condition benefited biomass, litter and SOM accumulation (Figs. 5 and 7), which resulted in the mostly positive wet-legacy impacts on R e (Figs. 3c,d and 6f).Conversely, the dry legacy decreased these pools (Figs. 5 and 7) and therefore resulted in the mostly negative dry-legacy impacts on R e (Figs. 3c,d and 6e).Contrary to our simulation results that dry legacies are mostly negative on SOM and R h , some field studies suggest that the labile C resulting from litter decomposition in a dry season may stimulate R h in the following wet season (Jenerette et al., 2008;Scott et al., 2009;Ma et al., 2012).This is likely because the labile soil C pool in the PALS model only accounts for ∼ 3 % of the total SOM and has a very short residence time (1.7 year; see Table S1); small amounts of seasonal labile C carryover therefore may not exert obvious legacy impacts on the total SOM pool size and R h across interannual and interdecadal scales.These results imply that the PPT legacy effects differ in direction and magnitude, depending on the type of C fluxes under consideration, the type of legacies (i.e., dry vs. wet), and the temporal scale of analysis.
Several lines of future research will likely be needed to continue improving the understanding of ecosystem legacy dynamics.Structural shifts in vegetation composition such as woody-plant encroachment (Potts et al., 2008;Scott et al., 2014), exotic-species invasion (Hamerlynck et al., 2010;Scott et al., 2010), and changes in microbial communities (de Vries et al., 2012;Evans and Wallenstein, 2012;Collins et al., 2014) may also interact with the biogeochemical processes to shape the PPT legacy effects on the temporal variability of dryland C fluxes.Furthermore, we need to better understand the legacy effects of extreme events such as the cool-GS drought in 2006 (see Fig. 1a) so that these important events can be adequately simulated.This cool-GS drought may have caused increased plant mortality as reported for a semi-desert grassland near our study site (Scott et al., 2010;Hamerlynck et al., 2013), but that is poorly represented in the model and may have caused the overestimation of the modeled GEP in comparison with the observation (see Fig. S3b).Finally, our approach that uses a highly resolved process model provides information complementary to contrasting analytical approaches that evaluate ecosystem responses to statistical rainfall regimes (Rodrigo-Iturbe et al., 2006;Katul et al., 2007;Porporato and Rodríguez-Iturbe, 2013).Improvement of these alternative modeling approaches is needed to understand both general and specific ecosystem responses to changing PPT regimes at temporal scales from events to decades.

Conclusions
We learned through this simulation analysis that (1) previous PPT conditions can impose substantial legacy impacts on the C balance of dryland ecosystems, with dry legacies fostering more current C sequestration and wet legacies causing more current C release; (2) the responses of ecosystem C fluxes to the simulated dry and wet legacies are mostly opposite in direction and asymmetrical in magnitude, with dry legacies being greater for GEP than for R e and wet legacies being greater for R e than for GEP; (3) the carryover of N soil is mainly responsible for the GEP responses, and the carryovers of biomass, litter and SOM are mainly responsible for the R e responses; and (4) the simulated PPT legacy effects can last for several years even with a 1-year PPT change, and therefore the direction and magnitude of interannual PPT legacy effects are less predictable than interdecadal ones.These simulation results suggest that dryland ecosystems such as these in the southwestern USA may emit more C that was sequestered in the past into the atmosphere with the predicted drying trend in the region (Seager et al., 2007;Solomon et al., 2007).The temporal variability of ecosystem C fluxes may be further intensified in the region due to the increasing PPT variability and the associated legacy impacts.
The Supplement related to this article is available online at doi:10.5194/bg-13-425-2016-supplement.

Figure 1 .
Figure 1.Precipitation characteristics in the 30 years (1981-2010) at the Santa Rita mesquite savanna site.(a) Annual and seasonal precipitation amount; (b) frequency distribution of daily rainfall; (c) mean and maximum between-event interval (BEI).Horizontal lines within (a) indicate mean annual and seasonal precipitation.The warm growing season (warm-GS) is from July through September, the cool dry season (cool-DS) from October to November, the cool growing season (cool-GS) from December through March, and the warm dry season (warm-DS) from April through June.Error bars in panel (c) represent standard deviations, and n is the number of rain event pairs used to calculate the between-event interval in the 30 years.

Figure 2 .
Figure 2. Comparison of the model-simulated water and carbon fluxes with the eddy-covariance observations at the mesquite savanna site.Left panels show the comparison between the modeled and observed fluxes in 4 calibration (2007-2007; solid dots) and 3 validation years (2008-2010; open dots).Right panels show the relationships of the simulated (solid dots) and observed (open dots) fluxes with precipitation in the 7 years (2004-2010).R 2 is the coefficient of determination describing the proportion of the variance in measured fluxes explained by the model for the left panels or that explained by precipitation for the right panels.AET represents actual evapotranspiration, GEP gross ecosystem production, R e total ecosystem respiration, and NEP net ecosystem production.

Figure 3 .
Figure 3. Interdecadal legacy effects of changing the previousperiod (1981-1994) precipitation on the cumulative carbon fluxes of the current period (1995-2010).Interdecadal legacy effects on carbon fluxes (e.g., NEP) are calculated as the difference between the current-period flux with previous-period PPT changes and that without previous-period PPT changes.Dashed lines with open symbols represent different levels of decreasing the current-period precipitation (left panels).Solid lines with filled symbols represent increasing the current-period precipitation (right panels).

Figure 4 .
Figure 4. Spearman correlations of interdecadal precipitation legacy effects with the precipitation difference between periods ( PPT).Interdecadal PPT is calculated as the mean PPT of the current period (1995-2010) minus that of the previous period(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994).Interdecadal legacy effects on carbon fluxes (e.g., NEP) are calculated as the difference between the current-period flux with previous-period PPT changes and that without previousperiod PPT changes.Sample size is 41 for the wet-to-dry period transition (left panels) and 23 for the dry-to-wet period transition (right panels).GEP represents gross ecosystem production, R e ecosystem respiration, and NEP net ecosystem production.R 2 is the coefficient of determination, and P is probability.

Figure 5 .
Figure 5. Interdecadal precipitation legacy effects on the resource pool dynamics.Left panels show the resource pool responses under a 30 % decrease, while right panels show those under a 30 % increase in the precipitation (PPT) of the current period from 1995 to 2010.Legacy effects on pool size (e.g., Biomass) are quantified as the difference between the current-period pool size with previousperiod PPT change and that without previous-period PPT change.Dashed lines represent a 30 % decrease, while solid lines represent a 30 % increase in the PPT of the previous period from 1981 to 1994.SOM represents soil organic matter, N soil soil mineral nitrogen, and SWC soil water content.

Figure 6 .
Figure 6.Interannual precipitation legacy effects on the ecosystem carbon fluxes.(a) and (b) show the lasting duration of dry (left panels) and wet (right panels) legacies, respectively.The legacy lasting duration is quantified as the number of years during which the legacy impacts on ecosystem fluxes exist after a previous-year PPT change.(c) through (h) show the responses of gross ecosystem production (GEP), ecosystem respiration (R e ) and net ecosystem production (NEP) to dry (left panels) and wet (right panels) legacies.Bars in the background of (a) and (b) represent the previous-year PPT amount after a 30 % decrease and increase, respectively.

Figure 7 .
Figure 7. Interannual precipitation legacy effects on resource pool dynamics.Left panels show the legacy effects on pool dynamics in 2 representative wet years, while right panels for 2 representative dry years.Legacy effects on pool size (e.g., Biomass) are quantified as the difference between the current-year pool size with previousyear PPT change and that without previous-year PPT change.Solid lines represent a 30 % decrease, while dashed lines represent a 30 % increase in the previous-year precipitation (PPT).SOM represents soil organic matter, N soil soil mineral nitrogen, and SWC soil water content.