Rainfall intensification increases the contribution of rewetting pulses to soil heterotrophic respiration

Soil drying and wetting cycles promote carbon (C) release through large heterotrophic respiration pulses at rewetting, known as the “Birch” effect. Empirical evidence shows that drier conditions before rewetting and larger changes in soil moisture at rewetting cause larger respiration pulses. Because soil moisture varies in response to rainfall, these respiration pulses also depend on the random timing and intensity of precipitation. In addition to rewetting pulses, heterotrophic respiration continues during soil drying, eventually ceasing when soils are too dry to sustain microbial activity. The importance of respiration pulses in contributing to the overall soil heterotrophic respiration flux has been demonstrated empirically, but no theoretical investigation has so far evaluated how the relative contribution of these pulses may change along climatic gradients or as precipitation regimes shift in a given location. To fill this gap, we start by assuming that heterotrophic respiration rates during soil drying and pulses at rewetting can be treated as random variables dependent on soil moisture fluctuations, and we develop a stochastic model for soil heterotrophic respiration rates that analytically links the statistical properties of respiration to those of precipitation. Model results show that both the mean rewetting pulse respiration and the mean respiration during drying increase with increasing mean precipitation. However, the contribution of respiration pulses to the total heterotrophic respiration increases with decreasing precipitation frequency and to a lesser degree with decreasing precipitation depth, leading to an overall higher contribution of respiration pulses under future more intermittent and intense precipitation. Specifically, higher rainfall intermittency at constant total rainfall can increase the contribution of respiration pulses up to ∼10 % or 20 % of the total heterotrophic respiration in mineral and organic soils, respectively. Moreover, the variability of both components of soil heterotrophic respiration is also predicted to increase under these conditions. Therefore, with future more intermittent precipitation, respiration pulses and the associated nutrient release will intensify and become more variable, contributing more to soil biogeochemical cycling.

1 Introduction

Heterotrophic respiration pulses often occur after dry soils are wetted by rainfall or irrigation (Barnard et al., 2020; Borken and Matzner, 2009; Canarini et al., 2017; Jarvis et al., 2007; Kim et al., 2012). The respiration rates achieved at rewetting can be much higher than the rates maintained under permanently moist conditions, suggesting that the rewetting itself triggers a disproportionally high CO₂ production. Even if they are short-lived, these pulses can contribute a significant fraction of the annual CO₂ release (Kim et al., 2012; Li et al., 2004; Yan et al., 2014). Their occurrence had been documented as long ago as Birch (1958) – for which the phenomenon has been named the “Birch effect” – but they remain difficult to explain and predict.

Respiration pulses are larger when the change in soil moisture is larger and when the soil was drier before rewetting, as shown by observations under both laboratory (Birch, 1958; Fischer, 2009; Guo et al., 2014; Lado-Monserrat et al., 2014; Schaeffer et al., 2017; Williams and Xia, 2009) and field conditions (Cable et al., 2008; Carbone et al., 2011; Lopez-Ballesteros et al., 2016; Rubio and Detto, 2017; Unger et al., 2010; Yan et al., 2014). Besides CO₂ displacement at rewetting, several mechanisms linked to microbial processes have been postulated to explain these patterns (Barnard et al., 2020; Canarini et al., 2017; Kim et al., 2012; Schimel et al., 2007). It has been argued that cell lysis due to a rapid increase in water potential and subsequent consumption of the dead cells may cause the pulse (Bottner, 1985). Later measurements showed that little cell lysis occurs but that intracellular materials (osmolytes) can be released at rewetting, contributing to the respiration pulse (Fierer and Schimel, 2003). However, in some soils microbial cells become dormant during drying rather than accumulating osmolytes (Boot et al., 2013). It is thus possible that respiration pulses are triggered by a physical process associated with the rewetting event – possibly reestablishment of hydrologic connectivity between substrates and microorganisms (Manzoni et al., 2016), or physical disruption of soil aggregates releasing old organic matter (Homyak et al., 2018). Indeed, there is a strong correlation between the CO₂ production after rewetting and the amount of extractable organic C consumed, suggesting that extractable C accumulated during the previous dry period could fuel the respiration pulse (Canarini et al., 2017; Guo et al., 2014; Williams and Xia, 2009). It is likely that multiple mechanisms work in concert, shifting their relative importance under different conditions (Slessarev and Schimel, 2020).

The focus on the processes causing respiration pulses resulted in extensive work conducted under idealized laboratory conditions, in which soil moisture changes were controlled, typically following a regular pattern of drying and wetting (Fierer and Schimel, 2002; Miller et al., 2005; Shi and Marschner, 2014, 2015; Xiang et al., 2008). However, soil moisture varies randomly due to the stochastic nature of rainfall events (Katul et al., 2007; Rodriguez-Iturbe and Porporato, 2004), and this temporal variability can either promote or decrease soil organic C storage depending on its effects on soil microbes (Lehmann et al., 2020). Two features of soil moisture dynamics are particularly important because they directly affect the intensity of a respiration pulse – the duration of dry periods and the soil moisture increment at rewetting. Therefore, experimental designs based on regular cycles of drying and wetting do not allow exploration of how the stochastic nature of soil moisture fluctuations may affect respiration pulses. Capturing the effect of these stochastic fluctuations can be important as climatic changes are altering rainfall patterns – often lengthening the duration of droughts and increasing the intensity of the (less frequent) rainfall events (IPCC, 2012).

To quantify how the long-term mean heterotrophic respiration varies as a function of statistical rainfall properties (duration of dry periods and intensity), we developed a stochastic soil moisture and respiration model, parameterized using available respiration data. Specifically, we ask – how does variability in rainfall translate into variability in respiration pulses? How does the long-term mean contribution of respiration pulses vary along climatic gradients? These questions are motivated by the hypothesis that respiration pulses contribute a larger proportion of soil heterotrophic respiration under climates with more intermittent and intense rainfall events, compared to climates in which soil moisture variations are mild. If that is the case, future climatic conditions characterized by longer droughts and more intense rainfall events are expected to increase the overall role of respiration pulses in ecosystem C budgets.

2 Methods

2.1 Theory

The theoretical framework is illustrated in Fig. 1. We start from the premise that heterotrophic respiration follows changes in soil moisture during drying (R_d) and that respiration pulses occur immediately following rewetting. As such, respiration pulses depend on both the soil moisture at the end of the dry period and the soil moisture increase caused by rainfall (R_r). The stochasticity of rainfall timing and amount determines a range of possible durations of dry spells and soil moisture increments when rainfall occurs. As a result, respiration can be regarded as a stochastic process. To statistically characterize the two types of respiration, the statistical properties of both soil moisture and soil moisture changes at rewetting are needed. These statistical properties are included in the probability density function (PDF) of soil moisture and the joint PDF of soil moisture and its increase at rewetting. Both distributions are derived in Sect. 2.1.1. The PDF of respiration rates during drying and respiration pulses at rewetting are derived in Sect. 2.1.2 and 2.1.3, respectively. All symbols are defined in Table 1.

Figure 1 Schematic illustration of the theoretical framework developed to describe how the components of heterotrophic respiration change as a function of rainfall statistical properties. (a) Rainfall is treated as a stochastic process driving random fluctuations in soil moisture, which are captured by the probability density functions (PDF, indicated by p with a subscript for the variable of interest) of soil moisture (x) and soil moisture increments (y). (b) Respiration rate during drying (R_d) and respiration pulses at rewetting (R_r) respectively depend on soil moisture and on both soil moisture increments and soil moisture at the end of the dry period (x_d); based on the PDFs of x, y, and x_d, the PDFs of the two respiration components are obtained. (c) Using these PDF of respiration, long-term mean respiration rates during drying (〈R_d〉) and respiration pulses (〈R_r〉) are calculated, and their relations with the statistical properties of precipitation are analyzed.

Table 1 Symbol definitions and units. Symbol p_z represents the probability density function (PDF) of the stochastic variable z indicated in the subscript.

2.1.1 Soil moisture dynamics

Soil moisture varies in response to rainfall events and the subsequent loss of soil water by percolation below the rooting zone and evapotranspiration. The dynamics of soil moisture in the rooting zone (the most biogeochemically active soil layer) can be described by the mass balance equation (Laio et al., 2001; Rodriguez-Iturbe and Porporato, 2004),

$$\begin{array}{}\text{(1)}& n{Z}_{\mathrm{r}}{\displaystyle \frac{\mathrm{d}s}{\mathrm{d}t}}=P\left(t\right)-E\left(s\left(t\right)\right)-L\left(s\left(t\right),t\right),\end{array}$$

where s is the saturation level (i.e., the relative volumetric soil moisture), n is the soil porosity, Z_r is the rooting depth, and P, E, and L represent precipitation inputs, evapotranspiration rate, and the combination of water losses due to percolation below the rooting zone and surface runoff. Equation (1) is interpreted at the daily timescale. Given our aim to describe the statistical properties of respiration rather than the details of soil moisture dynamics, we simplify the soil moisture mass balance equation to a form that is analytically tractable. Thus, we assume that evapotranspiration is the dominant water loss when soil moisture is lower than a threshold s_l (equivalent to the soil field capacity), whereas runoff and deep percolation dominate above this threshold. Also, runoff and percolation are assumed to occur rapidly compared to the timescales of the soil dry-down (free drainage conditions), so that, after a precipitation event that brings soil moisture above the level s_l, soil moisture decreases instantaneously to s_l. For simplicity, evapotranspiration is modeled as a linear function of soil moisture (Porporato et al., 2004),

$$\begin{array}{}\text{(2)}& E={\displaystyle \frac{s-s_{\mathrm{w}}}{s_{\mathrm{l}}-s_{\mathrm{w}}}}{E}_{max}=x{E}_{max},\end{array}$$

where E_max is the maximum rate of evapotranspiration, s_w is the plant wilting point (below which ET becomes negligible), and s_l is the threshold above which runoff and percolation are dominant. In the second equality, a normalized soil moisture denoted by x is introduced to further simplify the notation. With these assumptions and definitions, s ranges between s_w and s_l, while the normalized soil moisture varies between 0 and 1.

Precipitation is treated as a marked Poisson process with mean frequency λ and rain-event depths exponentially distributed with mean α. At each rain event, soil moisture increases by an amount corresponding to the rain event depth (normalized by nZ_r), unless the depth exceeds the available soil storage capacity (i.e., nZ_r(s_l−s)). Assuming that rainfall exceeding this capacity is routed to runoff, the PDF of soil moisture increments due to a rain event, y, for a given soil moisture at the end of the dry period, x_d, is given by (Laio et al., 2001)

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{p}_y\left(y\mathrm{|}{x}_{\mathrm{d}}\right)& =\mathit{\theta }\left[\left(\mathrm{1}-x_{\mathrm{d}}\right)-y\right]\mathit{\gamma }{e}^{-\mathit{\gamma }y}\\ & +\mathit{\delta }\left[y-\left(\mathrm{1}-x_{\mathrm{d}}\right)\right]e^{-\mathit{\gamma }\left(\mathrm{1}-x_{\mathrm{d}}\right)},\end{array}\end{array}$$

where p_y(y|x_d) is the PDF of y conditional on soil moisture at the end of the dry period, x_d; θ[⋅] is the Heaviside step function; δ[⋅] is the Dirac delta function; and γ is a parameter group defined as $\mathit{\gamma }=\frac{n{Z}_{\mathrm{r}}\left(s_{\mathrm{l}}-s_{\mathrm{w}}\right)}{\mathit{\alpha }}$ (γ can be interpreted as the number of average rainfall events needed to replenish the plant-available soil water). The first term on the right-hand side of Eq. (3) represents the probability density of a soil moisture increase y equal to the rainfall depth (θ[⋅] is equal to 1 for $y<\mathrm{1}-x_{\mathrm{d}}$; zero otherwise). The second term represents the probability of a soil moisture increase from the value x_d to the soil field capacity (x=1). This term is also referred to as an “atom of probability” because δ[⋅] is equal to zero for all soil moisture increments, except y=1-x_d, at which $\mathit{\delta }\left[\cdot \right]=\mathrm{\infty }$.

With this stochastic description of precipitation events and further assuming stochastic stationary conditions, the PDF of the normalized soil moisture driven by the dynamics in Eq. (1) can be obtained analytically and reads (Porporato et al., 2004)

$$\begin{array}{}\text{(4)}& p_x\left(x\right)=C_{\mathrm{l}}{\displaystyle \frac{e^{-x\mathit{\gamma }}{x}^{-\mathrm{1}+\frac{\mathit{\lambda }}{\mathit{\eta }}}}{\mathit{\eta }}},\end{array}$$

where η is a parameter group defined as $\mathit{\eta }=\frac{E_{max}}{n{Z}_{\mathrm{r}}\left(s_{\mathrm{l}}-s_{\mathrm{w}}\right)}$. C_l is a normalization constant that guarantees that the area under p_x(x) between x=0 and 1 is 1,

$$\begin{array}{}\text{(5)}& C_{\mathrm{l}}={\displaystyle \frac{\mathit{\eta }{\mathit{\gamma }}^{\frac{\mathit{\lambda }}{\mathit{\eta }}}}{\mathrm{\Gamma }\left[\frac{\mathit{\lambda }}{\mathit{\eta }}\right]-\mathrm{\Gamma }\left[\frac{\mathit{\lambda }}{\mathit{\eta }},\mathit{\gamma }\right]}},\end{array}$$

where Γ[⋅] and $\mathrm{\Gamma }\left[\cdot ,\cdot \right]$ are the complete and incomplete gamma functions (defined in Table 1). The PDF of soil moisture is the basis to obtain the PDF of respiration during soil drying (Sect. 2.1.2).

The last distribution needed to calculate the statistical properties of soil respiration pulses (Sect. 2.1.3) is the joint PDF of soil moisture at the end of a dry period and soil moisture increase due to precipitation events, denoted by $p_{y,x_{\mathrm{d}}}\left(y,x_{\mathrm{d}}\right)$ (note that both y and x_d are stochastic variables in this joint PDF). Thanks to the properties of the Poisson process, the PDF of soil moisture at the end of the dry period is equal to the PDF of soil moisture at a generic time (Cox and Miller, 2001), i.e., $p_{x_{\mathrm{d}}}\left(x_{\mathrm{d}}\right)=p_x\left(x\right)$. Because precipitation does not depend on antecedent soil moisture conditions in this model, the PDF of soil moisture at the end of a dry period is independent of the PDF of the subsequent precipitation event and soil moisture increase. Thus, the joint PDF of x_d and y is given by the product of the PDFs of x_d (Eq. 4) and of y conditional to x_d (Eq. 3),

$$\begin{array}{}\text{(6)}& p_{y,x_{\mathrm{d}}}\left(y,x_{\mathrm{d}}\right)=p_{x_{\mathrm{d}}}\left(x_{\mathrm{d}}\right)p_y\left(y\mathrm{|}{x}_{\mathrm{d}}\right).\end{array}$$

2.1.2 Heterotrophic respiration during soil drying

During a dry period, the heterotrophic respiration rate decreases in response to the gradual decrease in soil moisture, following a concave-downward trend (Manzoni et al., 2012; Moyano et al., 2012). Consistent with the hydrologic model setup, we assume that the soil drains rapidly and hence does not remain under saturated conditions long enough to develop anoxic conditions. It is thus reasonable to assume that respiration declines between the soil field capacity (equivalent to s_l in this model) and a lower soil moisture threshold for microbial activity. This lower threshold corresponds to water potential levels around −15 MPa in sieved soil samples (Manzoni and Katul, 2014), but here we assume that respiration becomes much smaller than rates under well-watered conditions already at the plant wilting point s_w, i.e., at a water potential of −1.5 MPa. This assumption is motivated by the observation that in intact soil cores and under field conditions respiration stops in wetter conditions than at −15 MPa (e.g., −2.7 MPa; Carbone et al., 2011). Moreover, this allows us to keep the parameter number to a minimum, consistent with the minimal soil moisture balance model of Eqs. (1) and (2) and the overall idealized representation of soil heterotrophic respiration. The respiration decrease with a lower threshold s_w (corresponding to x=0) can be captured by a parabolic relation,

$$\begin{array}{}\text{(7)}& R_{\mathrm{d}}=R_{\mathrm{d},\max }\left(\mathrm{2}x-x^{\mathrm{2}}\right),\end{array}$$

where R_d denotes the respiration rate during drying, and R_d,max is the maximum respiration rate in the absence of rapid rewetting (i.e., R_d at x=1 or s=s_l). Using other monotonic and concave-downward relations between respiration and soil moisture would not qualitatively alter the results.

In Eq. (1), soil moisture is a random variable, whose PDF follows Eq. (4). Therefore, R_d from Eq. (7) is also a random variable, which can be obtained from the PDF of soil moisture using the derived distribution approach, also referred to as the Jacobian rule (Kottegoda and Rosso, 1998),

$$\begin{array}{}\text{(8)}& p_{R_{\mathrm{d}}}\left(R_{\mathrm{d}}\right)=p_x\left(x\left(R_{\mathrm{d}}\right)\right)\left|{\displaystyle \frac{\mathrm{d}x}{\mathrm{d}{R}_{\mathrm{d}}}}\right|,\end{array}$$

where on the right-hand side the PDF of soil moisture is evaluated at moisture values corresponding to given respiration values. This is done by inverting Eq. (7) and expressing x as a function of R_d,

$$\begin{array}{}\text{(9)}& x\left(R_{\mathrm{d}}\right)=\mathrm{1}-\sqrt{\mathrm{1}-{\displaystyle \frac{R_{\mathrm{d}}}{R_{\mathrm{d},\max }}}}.\end{array}$$

We note that Eq. (7) is monotonic in the domain $\mathrm{0}\le x\le \mathrm{1}$, which allows unambiguous definition of the inverse of R_d(x). Had we used a nonmonotonic R_d(x) function (e.g., for applications of this approach to soils experiencing long saturation periods), the derived distribution approach would have required splitting the x domain into two – one for each monotonic branch of R_d(x). In turn, Eq. (9) allows the calculation of the slope of the x(R_d) relation, which is also needed in Eq. (8),

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{\mathrm{d}x}{\mathrm{d}{R}_{\mathrm{d}}}}={\left(\mathrm{2}{R}_{\mathrm{d},\max }\sqrt{\mathrm{1}-{\displaystyle \frac{R_{\mathrm{d}}}{R_{\mathrm{d},\max }}}}\right)}^{-\mathrm{1}}.\end{array}$$

The PDF of R_d is thus obtained from Eqs. (8)–(10) as

$$\begin{array}{}\text{(11)}& p_{R_{\mathrm{d}}}\left(R_{\mathrm{d}}\right)=C_{\mathrm{l}}{\displaystyle \frac{e^{-\mathit{\gamma }\left(\mathrm{1}-\sqrt{\mathrm{1}-r_{\mathrm{d}}}\right)}{\left(\mathrm{1}-\sqrt{\mathrm{1}-r_{\mathrm{d}}}\right)}^{-\mathrm{1}+\frac{\mathit{\lambda }}{\mathit{\eta }}}}{\mathrm{2}\mathit{\eta }{R}_{\mathrm{d},\max }\sqrt{\mathrm{1}-r_{\mathrm{d}}}}},\end{array}$$

where the normalized respiration $r_{\mathrm{d}}=\frac{R_{\mathrm{d}}}{R_{\mathrm{d},\max }}$ is introduced to simplify the notation. This PDF can now be used to analytically calculate the long-term mean of R_d, denoted by 〈R_d〉,

$$\begin{array}{}\text{(12)}& \begin{array}{rl}\langle R_{\mathrm{d}}\rangle & ={\displaystyle \frac{R_{\mathrm{d},\max }{C}_{\mathrm{2}}}{{\mathit{\gamma }}^{\mathrm{2}}}}\left\{\mathrm{\Gamma }\left[\mathrm{2}+{\displaystyle \frac{\mathit{\lambda }}{\mathit{\eta }}},\mathit{\gamma }\right]-\mathrm{2}\mathit{\gamma }\mathrm{\Gamma }\left[\mathrm{1}+{\displaystyle \frac{\mathit{\lambda }}{\mathit{\eta }}},\mathit{\gamma }\right]\right.\\ & \left.-{\displaystyle \frac{\mathit{\lambda }\left(\mathit{\eta }-\mathrm{2}\mathit{\gamma }\mathit{\eta }+\mathit{\lambda }\right)}{{\mathit{\eta }}^{\mathrm{2}}}}\mathrm{\Gamma }\left[{\displaystyle \frac{\mathit{\lambda }}{\mathit{\eta }}}\right]\right\},\end{array}\end{array}$$

where for convenience the parameter group $C_{\mathrm{2}}={\left(\mathrm{\Gamma }\left[\frac{\mathit{\lambda }}{\mathit{\eta }}\right]-\mathrm{\Gamma }\left[\frac{\mathit{\lambda }}{\mathit{\eta }},\mathit{\gamma }\right]\right)}^{-\mathrm{1}}$ is defined. The standard deviation of R_d, denoted by ${\mathit{\sigma }}_{R_{\mathrm{d}}}$, can not be obtained analytically, but it can be calculated through numerical integration of Eq. (11).

2.1.3 Heterotrophic respiration pulses at rewetting

Heterotrophic respiration pulses at rewetting are caused by mineralization of available C and microbial products at the end of the dry period, which in turn depend on how intense the rewetting event was. As a result of these processes, in a given soil, rewetting events depend on both soil moisture before the rewetting x_d and the change in soil moisture y (Birch, 1958; Lado-Monserrat et al., 2014). This relation can be captured by the empirical function (justified and parameterized in Sect. 2.2.1)

$$\begin{array}{}\text{(13)}& R_{\mathrm{r}}=R_{\mathrm{r},\max }{\displaystyle \frac{y}{\mathrm{1}+\frac{x_{\mathrm{d}}}{b}}}\mathit{\theta }\left[\left(\mathrm{1}-x_{\mathrm{d}}\right)-y\right],\end{array}$$

where R_r,max is the largest respiration pulse possible, which is achieved when an initially dry soil reaches saturation, i.e., y=1 and x_d=0. The parameter b accounts for the effect of antecedent soil moisture conditions – for a given value of x_d, the respiration pulse increases with increasing b. The last term in Eq. (13) is a Heaviside function limiting the relation between R_r and y to conditions in which soil moisture at most fills the available pore space (as in Eq. 3, θ[⋅] is equal to 1 only when $y<\mathrm{1}-x_{\mathrm{d}}$). If before the rain event soil moisture is at the plant wilting point (x_d=0) and the precipitation event is sufficient to reach s_l (i.e., $y=x-x_{\mathrm{d}}=\mathrm{1}$), the maximum respiration pulse is attained and $R_{\mathrm{r}}=R_{\mathrm{r},\max }$. Here, R_r represents an amount of C respired when the rewetting event occurs, so its dimensions differ from those of the respiration rate during drying, R_d; these two quantities are combined in the total heterotrophic respiration rate in Sect. 2.1.5.

Because both y and x_d are random variables that follow the PDF of Eq. (6), R_r should also be regarded as a random variable following its own PDF. Different from the PDF of R_d, which was obtained from the univariate PDF of soil moisture, the PDF of R_r has to be derived from the joint PDF of y and x_d. The derived distribution approach can still be used, but it requires the determinant of the Jacobian matrix of the transformation from y and x_d to R_r (Kottegoda and Rosso, 1998). To proceed, it is first convenient to introduce an auxiliary variable X=x_d, which is used together with Eq. (13) to find the transformation from the original variables y and x_d to R_r and X,

$$\begin{array}{}\text{(14)}& \begin{array}{rl}& \left\{\begin{array}{c}X=x_{\mathrm{d}}\\ R_{\mathrm{r}}=R_{\mathrm{r},\max }\frac{y}{\mathrm{1}+\frac{x_{\mathrm{d}}}{b}}\end{array}\right.\Rightarrow \left\{\begin{array}{c}{x}_{\mathrm{d}}=X\\ y=\frac{R_{\mathrm{r}}}{R_{\mathrm{r},\max }}\left(\mathrm{1}+\frac{X}{b}\right)\end{array}\right.\\ & \mathrm{for}y<\mathrm{1}-x_{\mathrm{d}},\end{array}\end{array}$$

where the inequality limits the soil moisture increments as the Heaviside function in Eq. (13). Second, the system on the left of Eq. (14) is inverted to express the original variables as a function of the transformed variables (reported on the right of Eq. 14), similar to the inversion done in Eq. (9). Third, we calculate the Jacobian matrix,

$$\begin{array}{}\text{(15)}& \mathbf{J}=\left[\begin{array}{cc}\frac{\partial x_{\mathrm{d}}}{\partial X}& \frac{\partial x_{\mathrm{d}}}{\partial R_{\mathrm{r}}}\\ \frac{\partial y}{\partial X}& \frac{\partial y}{\partial R_{\mathrm{r}}}\end{array}\right]=\left[\begin{array}{cc}\mathrm{1}& \mathrm{0}\\ \frac{R_{\mathrm{r}}}{R_{\mathrm{r},\max }}\frac{\mathrm{1}}{b}& \frac{\mathrm{1}}{R_{\mathrm{r},\max }}\left(\mathrm{1}+\frac{X}{b}\right)\end{array}\right],\end{array}$$

and the determinant of the Jacobian,

$$\begin{array}{}\text{(16)}& \left|\mathbf{J}\right|={\displaystyle \frac{\mathrm{1}}{R_{\mathrm{r},\max }}}\left(\mathrm{1}+{\displaystyle \frac{X}{b}}\right).\end{array}$$

Fourth, the joint PDF of the variables X and R_r is obtained using the derived distribution approach,

$$\begin{array}{}\text{(17)}& p_{X,R_{\mathrm{r}}}\left(X,R_{\mathrm{r}}\right)=p_{y,x_{\mathrm{d}}}\left(y\left(X,R_{\mathrm{r}}\right),x_{\mathrm{d}}\left(X,R_{\mathrm{r}}\right)\right)\left|\mathbf{J}\right|,\end{array}$$

where as in Sect. 2.1.2 all the terms on the right-hand side only depend on X and R_r, and $p_{y,x_{\mathrm{d}}}$ is given by Eq. (6). Finally, to obtain the (marginal) PDF of R_r, the joint PDF in Eq. (17) is integrated over all possible values of X,

$$\begin{array}{}\text{(18)}& \begin{array}{rl}{p}_{R_{\mathrm{r}}}\left(R_{\mathrm{r}}\right)& =\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{p}_{X,R_{\mathrm{r}}}\left(X,R_{\mathrm{r}}\right)\mathrm{d}X\\ & ={\displaystyle \frac{C^{\prime }}{\left(b+r_r\right)R_{\mathrm{r},\max }}}\mathit{\{}{e}^{-\mathit{\gamma }}{\left[{\displaystyle \frac{\mathit{\gamma }\left(\mathrm{1}-r_r\right)}{\mathrm{1}+\frac{r_r}{b}}}\right]}^{\frac{\mathit{\lambda }}{\mathit{\eta }}}{\displaystyle \frac{\mathrm{1}+b}{\mathrm{1}-r_r}}\\ & +e^{-\mathit{\gamma }{r}_r}{\left({\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\frac{r_r}{b}}}\right)}^{\frac{\mathit{\lambda }}{\mathit{\eta }}}\left[\mathit{\gamma }\left(b+r_r\right)\left(\mathrm{\Gamma }\left[{\displaystyle \frac{\mathit{\lambda }}{\mathit{\eta }}}\right]\right.\right.\\ & \left.-\mathrm{\Gamma }\left[{\displaystyle \frac{\mathit{\lambda }}{\mathit{\eta }}},\mathit{\gamma }\left(\mathrm{1}-r_r\right)\right]\right)+\mathrm{\Gamma }\left[\mathrm{1}+{\displaystyle \frac{\mathit{\lambda }}{\mathit{\eta }}}\right]\\ & \left.-\mathrm{\Gamma }\left[\mathrm{1}+{\displaystyle \frac{\mathit{\lambda }}{\mathit{\eta }}},\mathit{\gamma }\left(\mathrm{1}-r_r\right)\right]\right]\mathit{\}},\end{array}\end{array}$$

where on the right-hand side the normalized respiration pulse $r_{\mathrm{r}}=\frac{R_{\mathrm{r}}}{R_{\mathrm{r},\max }}$ is introduced to simplify the notation, and as before $C_{\mathrm{2}}={\left(\mathrm{\Gamma }\left[\frac{\mathit{\lambda }}{\mathit{\eta }}\right]-\mathrm{\Gamma }\left[\frac{\mathit{\lambda }}{\mathit{\eta }},\mathit{\gamma }\right]\right)}^{-\mathrm{1}}$. Due to the complexity of Eq. (18), the long-term mean and standard deviation of R_r, respectively denoted by 〈R_r〉 and ${\mathit{\sigma }}_{R_{\mathrm{r}}}$, need to be obtained via numerical integration.

2.1.4 Rewetting pulses only dependent on soil moisture change

It is useful to consider respiration pulses that only depend on the soil moisture increments; i.e., b≫x_d. In this case, Eq. (13) reduces to $R_{\mathrm{r}}=y{R}_{\mathrm{r},\max }$ (i.e., $y=R_{\mathrm{r}}/R_{\mathrm{r},\max })$ – equivalent to always having a completely dry soil before rewetting. Thanks to the simplicity of the respiration pulse equation, $p_{R_{\mathrm{r}}}\left(R_{\mathrm{r}}\right)$ can be obtained as a derived distribution from the marginal PDF of the soil moisture changes y,

$$\begin{array}{}\text{(19)}& p_y\left(y\right)=\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{p}_y\left(y\mathrm{|}{x}_{\mathrm{d}}\right)\mathrm{d}{x}_{\mathrm{d}}=\left[\mathrm{1}+\mathit{\gamma }\left(\mathrm{1}-y\right)\right]e^{-y\mathit{\gamma }},\end{array}$$

where p_y(y|x_d) is from Eq. (3). The $p_{R_{\mathrm{r}}}\left(R_{\mathrm{r}}\right)$ is then obtained as

$$\begin{array}{}\text{(20)}& \begin{array}{rl}{p}_{R_{\mathrm{r}}}\left(R_{\mathrm{r}}\right)& =p_y\left(y\left(R_{\mathrm{r}}\right)\right)\left|{\displaystyle \frac{\mathrm{d}y}{\mathrm{d}{R}_{\mathrm{r}}}}\right|\\ & ={\displaystyle \frac{\mathrm{1}+\mathit{\gamma }\left(\mathrm{1}-\frac{R_{\mathrm{r}}}{R_{\mathrm{r},\max }}\right)}{R_{\mathrm{r},\max }}}{e}^{-\frac{\mathit{\gamma }{R}_{\mathrm{r}}}{R_{\mathrm{r},\max }}}.\end{array}\end{array}$$

Thanks to the simplicity of Eq. (20), in this particular case the long-term mean and standard deviation of the respiration pulses are found analytically,

$$\begin{array}{}\text{(21)}& {\displaystyle }\langle R_{\mathrm{r}}\rangle ={\displaystyle \frac{R_{\mathrm{r},\max }}{{\mathit{\gamma }}^{\mathrm{2}}}}\left(e^{-\mathit{\gamma }}+\mathit{\gamma }-\mathrm{1}\right),\text{(22)}& {\displaystyle }{\mathit{\sigma }}_{R_{\mathrm{r}}}={\displaystyle \frac{R_{\mathrm{r},\max }}{{\mathit{\gamma }}^{\mathrm{2}}}}\sqrt{\left(\mathit{\gamma }-\mathrm{2}\right)\mathit{\gamma }-\mathrm{1}+\mathrm{2}{e}^{-\mathit{\gamma }}\left(\mathrm{1}+\mathit{\gamma }+{\mathit{\gamma }}^{\mathrm{2}}\right)-e^{-\mathrm{2}\mathit{\gamma }}}.\end{array}$$

Thus, when respiration pulses are simply proportional to the soil moisture change at rewetting, their mean only depends on the maximum pulse size R_r,max and the ratio of soil water storage capacity and mean precipitation depth (i.e., the parameter group $\mathit{\gamma }=\frac{n{Z}_{\mathrm{r}}\left(s_{\mathrm{l}}-s_{\mathrm{w}}\right)}{\mathit{\alpha }}$).

2.1.5 Combining heterotrophic respiration during soil drying and at rewetting

The total mean heterotrophic respiration rate is given by the sum of the mean respiration rate during soil drying 〈R_d〉 (Eq. 12; expressed in grams of carbon per square meter per day) and the mean rate of respiration resulting from the sequence of rewetting pulses over the study period (denoted by $\langle R_{\mathrm{r}}^{*}\rangle$ and also expressed in grams of carbon per square meter per day). The $\langle R_{\mathrm{r}}^{*}\rangle$ is calculated as the mean amount of respired carbon (〈R_r〉 from Eq. 18, expressed in grams of carbon per square meter) divided by the mean rainfall inter-arrival time, 1∕λ (expressed in days),

$$\begin{array}{}\text{(23)}& \langle R_{\mathrm{r}}^{*}\rangle =\mathit{\lambda }\langle R_{\mathrm{r}}\rangle .\end{array}$$

The mean total heterotrophic respiration rate is then obtained as

$$\begin{array}{}\text{(24)}& \langle R_{\mathrm{t}}\rangle =\langle R_{\mathrm{d}}\rangle +\langle R_{\mathrm{r}}^{*}\rangle .\end{array}$$

In what follows, the ratio of respiration pulse to total respiration (i.e., $\langle R_{\mathrm{r}}^{*}\rangle /\langle R_{\mathrm{t}}\rangle$) will also be considered, to evaluate the overall contribution of respiration pulses.

2.2 Data analysis

2.2.1 Laboratory incubation data for model calibration

The phenomenological respiration models in Eqs. (7) and (13) require knowledge of three parameters: the heterotrophic respiration rate at the soil field capacity (R_d,max), the maximum respiration pulse size (R_r,max), and the sensitivity of the respiration pulse to the initial soil moisture (b). To estimate these three parameters, we selected datasets where both the soil moisture before rewetting and the soil moisture increments were manipulated (Fischer, 2009; Guo et al., 2014; Lado-Monserrat et al., 2014). All data reported in these three publications were used, except data from the litter-amended soils in Lado-Monserrat et al. (2014) (we chose to focus on “natural” conditions) and data from small (y<0.3) rewetting events in Fischer (2009) (they exhibited small respiration peaks despite nearly stable soil moisture). The reported respiration amounts at rewetting were corrected to isolate the pulse size (R_r) from the respiration that would have occurred at constant soil moisture (R_d). This was done by calculating R_d,max from control soil samples kept constantly wet (Guo et al., 2014) or from the post-pulse respiration rate before soil moisture started to decline in experiments where drying was allowed in all samples (Lado-Monserrat et al., 2014). In contrast, respiration pulses had already been isolated by Fischer (2009). The last step of the parameter estimation involved fitting Eq. (13) to the data using a nonlinear least-square algorithm (fminunc function in MATLAB, R2018b, MathWorks, Inc.).

Because respiration amounts and rates in these laboratory incubations were expressed respectively in micrograms per gram and micrograms per gram per day (or on a per-unit soil organic C basis), units were converted to gram per square meter and gram per square meter per day using bulk densities and sampling depths reported in the original publications (results are shown in Table 2).

Table 2 Characteristics of the selected mineral and organic soil samples; estimates of the respiration model parameters in Eq. (13) (R_d,max: maximum respiration rate at the soil field capacity, R_r,max: maximum respired carbon at rewetting) and coefficients of determination (R²) for the least-square fit of the data (see also Fig. 2).

2.2.2 Field data for model validation

In addition to estimating the values of the three parameters in Eqs. (7) and (13), we validated the results from the whole stochastic model by comparing the predicted long-term mean heterotrophic respiration rates to observations along a rainfall manipulation gradient in a semiarid steppe (Zhang et al., 2017b, 2019). Briefly, the precipitation gradient was established by excluding 30 % and 60 % of precipitation with rain shelters and by increasing precipitation by 30 % and 60 % through irrigation. By design, only precipitation amounts (not timing) were altered, resulting in five mean rainfall depths α=2.6, 3.9, 5.1, 6.4, and 7.6 mm. Mean evapotranspiration rates, soil moisture, and heterotrophic respiration rates along the rainfall gradient were obtained from the published supplementary materials in Zhang et al. (2019) or from the Dryad dataset by Zhang et al. (2017a). Hydrologic parameters that were not provided were estimated as follows. The maximum evapotranspiration rate (assumed equal to the potential evapotranspiration) and the mean rainfall frequency were estimated from May–August CRU data at the rainfall manipulation site ($E_{max}=\mathrm{4.3}$ mm d⁻¹ and λ=0.41 d⁻¹). The soil at the site has a sandy loam texture (Bingwei Zhang, personal communication, 2019), and soil properties were obtained accordingly: n=0.42, s_w=0.11, s_l=0.52 (Table 2.1 in Rodriguez-Iturbe and Porporato, 2004). Finally, the rooting depth Z_r=0.2 m was estimated as the soil depth above which approximately 70 % of belowground productivity occurs, based on data from Zhang et al. (2020).

Regarding the parameters of the rewetting respiration function (Eq. 13), we assumed $R_{\mathrm{d},\max }=\mathrm{2}$ gC m⁻² d⁻¹ and b=0.1. These values are deemed reasonable for mineral soils based on Table 2 and accounting for a rooting depth about double the sampling depth of the incubation experiments (which doubles the R_d,max values in Table 2). Without specific information on respiration pulse sizes, we let R_r,max vary over a wide range. Additionally, we tested the simplified respiration model (Sect. 2.1.4), which does not require any assumption on b, against the same total heterotrophic respiration dataset.

3 Results

3.1 Dependence of heterotrophic respiration at rewetting on soil moisture

Laboratory incubation data were used to parameterize the functions linking heterotrophic respiration to soil moisture. As expected, the respiration pulses at rewetting depend on both rewetting intensity (y) and pre-wetting soil moisture (x_d), and this relation is well-characterized by Eq. (13) (Fig. 2). In Fig. 2, respiration pulses at rewetting are normalized by the amount of organic C in each soil to facilitate comparisons. However, after accounting for variations in organic C content, bulk density, and soil layer depth, the values of R_r,max and R_d,max per unit ground area are higher in the organic soils than in mineral soils (Table 2) and so is the ratio between R_r,max and R_d,max. The sensitivity parameter b shows milder variation across soils than the other parameters, with an average value b≈0.1. Based on this data analysis, in the following theoretical exploration we set parameter values intermediate between the extremes reported in Table 2 (i.e., $R_{\mathrm{r},\max }=\mathrm{5}$ gC m⁻², $R_{\mathrm{d},\max }=\mathrm{1}$ gC m⁻² d⁻¹, and b=0.1). In addition, we explore how the contribution of respiration pulses varies between mineral vs. organic soils, using the average parameter values reported in Table 2.

Figure 2 Relations between respiration pulse size (R_r, normalized by soil organic C content) and pre-wetting soil moisture (x_d) and soil moisture increment at rewetting (y), for six soils; top row: mineral soils; bottom row: organic soils. Symbols represent measured respiration pulses, and surfaces are fitted R_r functions from Eq. (13) (soil characteristics, fitting parameters, and data sources are reported in Table 2).

3.2 General model behavior

Figure 3 shows two examples of the simulated trajectories of soil moisture and heterotrophic respiration, for contrasting climatic conditions (more frequent precipitation in the left panels than in the right panels). It is important to note that in this comparison across climatic conditions (and in the comparisons that follow), the maximum respiration R_r,max and R_d,max are fixed, while in reality they are likely proportional to soil organic C availability, which in turn is the result of a long-term and soil-moisture-dependent balance between C inputs from vegetation and respiration (this limitation is discussed in Sect. 4.2). Respiration rates during dry periods follow soil moisture changes, declining as soil dries and returning to higher levels at rewetting (Fig. 3b, f). In addition to this rewetting-induced restoration of high respiration rates, rewetting causes CO₂ emission pulses, represented by vertical bars. Under the wetter climate (Fig. 3b), respiration pulses are more frequent than under the dry climate (Fig. 3f) because of the higher precipitation frequency. However, most of the respiration pulses are small because soil moisture increments at rewetting are often limited by the available soil pore space, and a relatively large fraction of precipitation is lost to runoff and deep percolation. In contrast, under dryer conditions, changes in soil moisture are large because on average soil moisture is low and the pore space is rarely filled up completely. As a result, the fewer respiration pulses can be larger under dry than under wet conditions.

Figure 3 Example of the dynamics of soil moisture and respiration for a wet (a–d) and a dry climate (e–h). Top panels (a) and (e) show the simulated trajectories of normalized soil moisture x; the middle panels (b) and (f) show the trajectories of respiration during dry periods (red solid curves, R_d) and the respiration pulse at rewetting (black vertical bars, R_r, on the same scale despite different units); the bottom panels (c), (d), (g), and (h) show the probability density functions of R_r and R_d ($p_{R_{\mathrm{r}}}\left(R_{\mathrm{r}}\right)$ and $p_{R_{\mathrm{d}}}\left(R_{\mathrm{d}}\right)$, respectively) overlapped to the histograms of the simulated data. In this figure, $R_{\mathrm{r},\max }=\mathrm{5}$ gC m⁻², $R_{\mathrm{d},\max }=\mathrm{1}$ gC m⁻² d⁻¹, b=0.1, γ=5, η=0.1 d⁻¹, and λ=0.3 and 0.1 d⁻¹ (a–d and e–h, respectively).

The bottom panels in Fig. 3 show the PDF of respiration for the same two climatic conditions analyzed in the upper panels. While the PDF of R_r is positively skewed regardless of climate (but with heavier tails under dry conditions, Fig. 3c, g), the PDF of R_d is strongly affected by climatic conditions – the probability of high values for R_d is higher under wet conditions (negatively skewed PDF) and lower under dry conditions (positively skewed PDF, Fig. 3d, h). This pattern is caused by the prevalence of high soil moisture values in the wet climate scenario, which maintain relatively high R_d. Figure 3c, d, g, h also show that the theoretical PDF (Eqs. 11 and 18) matches perfectly to the distribution of the numerically simulated data. The shape of the theoretical PDF of R_d in Fig. 3h might seem incorrect, as it increases sharply at high respiration values. This increase is due to the flat derivative of the R_d–soil moisture relation (Eq. 7), which causes an asymptote in the PDF at $R_{\mathrm{d}}=R_{\mathrm{d},\max }$ (Eq. 11). However, the area under this spike is vanishingly small when climatic conditions are dry as in the example of Fig. 3e–h, so that it is highly unlikely to have any respiration value around R_d,max.

3.3 Model test under field conditions

Field data were used to test whether the hydrologic and soil respiration models could capture trends in the mean evapotranspiration and heterotrophic respiration along a precipitation gradient (Fig. 4). The trend of the mean evapotranspiration rate with increasing mean rainfall depth was captured reasonably well (Fig. 4a), considering that no formal calibration was conducted, and all parameters were estimated based on independent information. Similarly, the model correctly predicts the trend in soil moisture (not shown), but with an overestimation bias around 0.05–0.1 (in terms of normalized soil moisture x). This overestimation is expected, because soil moisture had been measured in the drier top 0.1 m of soil, while the model considers average soil moisture over a 0.2 m depth. Also the trend in total heterotrophic respiration is predicted correctly by the full model, which explains 77 % of the variance in the respiration data (black curve in Fig. 4b). Calibrating the two parameters of Eq. (13) and R_d,max would allow a better fit, but since the goal here is to provide a qualitative model validation and not a quantitative performance assessment, we deem the model suitable for the following theoretical analyses.

Figure 4 Comparison of model results (curves) and observations (open circles) along an experimental rainfall gradient where the mean precipitation depth (α) was manipulated: (a) mean evapotranspiration rate 〈E〉, (b) mean total heterotrophic respiration rate 〈R_t〉, and (c) fraction of the total heterotrophic respiration rate due to rewetting pulses $\langle R_{\mathrm{r}}^{*}\rangle /\langle R_{\mathrm{t}}\rangle$. In panels (b) and (c), the dotted curves and shaded area indicate the variation caused by changes in R_r,max between 5 and 35 g m⁻² around the central value (solid curves) of 25 g m⁻²; the red curves indicate results using the simplified rewetting respiration model (Eq. 21). Parameter values are described in Sect. 2.2.2.

We also tested the simpler version of the model, in which respiration pulses only depend on the soil moisture increment. Without the effect of pre-wetting soil moisture, this version predicts higher mean respiration than the full model (red lines in Fig. 4b) and a higher contribution of rewetting respiration to the total heterotrophic respiration (red lines in Fig. 4c).

3.4 Dependence of heterotrophic respiration on rainfall statistical properties

Figure 5 shows the predicted effect of precipitation regimes on heterotrophic respiration during drying and at rewetting (Fig. 5a, b) on the total heterotrophic respiration rate (Fig. 5c) and on the fraction of respiration contributed by rewetting pulses (Fig. 5d). As in Fig. 3, R_r,max and R_d,max are fixed to focus on the role of climatic conditions, so the patterns shown in Fig. 5 should be interpreted as changes of mean respiration rates along gradients of precipitation frequency (λ) and mean depth (α) for given soil organic C stocks. Because in this minimal model the mean precipitation rate is given by 〈P〉=αλ, precipitation can be increased by assuming more frequent rain events (i.e., increasing λ), larger events (i.e., increasing α), or both. Any of these changes increase mean respiration during drying and at rewetting (Fig. 5a, b). As 〈R_d〉 increases with precipitation more than $\langle R_{\mathrm{r}}^{*}\rangle$, the relative contribution of respiration pulses to the total respiration rate, $\langle R_{\mathrm{r}}^{*}\rangle /\langle R_{\mathrm{t}}\rangle$, tends to decrease from drier to wetter conditions, especially when rain events become more frequent (as opposed to more intense) (Fig. 5d). This pattern is caused by the relatively larger respiration pulses occurring when soils are dry and rewetting causes large soil moisture increments (compare examples in Fig. 3b and f). Moreover, the relative change of $\langle R_{\mathrm{r}}^{*}\rangle /\langle R_{\mathrm{t}}\rangle$ is smaller than the change in 〈R_d〉 or $\langle R_{\mathrm{r}}^{*}\rangle$ as precipitation regimes are varied.

Figure 5 Effect of precipitation statistical properties (mean event frequency λ and depth α) on the mean heterotrophic respiration rates during dry periods 〈R_d〉 (a) and at rewetting $\langle R_{\mathrm{r}}^{*}\rangle$ (b), the mean total respiration rate 〈R_t〉 (c), and the fraction of the total heterotrophic respiration rate due to rewetting pulses $\langle R_{\mathrm{r}}^{*}\rangle /\langle R_{\mathrm{t}}\rangle$ (d). The white contour curves indicate combinations of λ and α that generate different annual precipitation rates ($\langle P\rangle =\mathit{\alpha }\mathit{\lambda }=\mathrm{0.25}$, 0.5, 1, and 2 m yr⁻¹ from dotted to solid lines). Other parameter values: $R_{\mathrm{r},\max }=\mathrm{5}$ gC m⁻², $R_{\mathrm{d},\max }=\mathrm{1}$ gC m⁻² d⁻¹, b=0.1, Z_r=0.3 m, n=0.5, s_w=0.2, s_l=0.7, $E_{max}=\mathrm{0.0037}$ m d⁻¹.

Not only the mean respiration rates, but also the variability of both respiration rates during drying and respiration pulses at rewetting vary with hydroclimatic conditions (Fig. 6). The standard deviation of R_d exhibits maxima at intermediate α when λ is fixed and at intermediate λ when α is fixed (Fig. 6a). This pattern is due to a shift in the shape of the PDF of R_d when moving from dry to wet conditions. Under dry conditions, the PDF of R_d has relatively low variance and is negatively skewed (Fig. 3h); as conditions become wetter the PDF flattens and the variance increases, and finally under wet conditions the PDF transitions again to a low-variance but positively skewed PDF (Fig. 3d). In contrast, the PDF of R_r is always positively skewed, with variance decreasing with increasing rainfall frequency (Fig. 6b; compare examples in Fig. 3c and g). However, increasing α for fixed λ is predicted to increase the variance of R_r. The coefficients of variation (CV) of R_d and R_r vary less than the corresponding standard deviations and tend to decrease as conditions move from dry to wet (Fig. 6c, d). Specifically, the CV of R_d decreases with both increasing λ and increasing α. In contrast, the CV of R_r is nearly independent of α but decreases with increasing λ.

Figure 6 Effect of precipitation statistical properties (mean event frequency λ and depth α) on the standard deviations of heterotrophic respiration rates during dry periods ${\mathit{\sigma }}_{R_{\mathrm{d}}}$ (a) and respiration pulses at rewetting ${\mathit{\sigma }}_{R_{\mathrm{r}}}$ (b), on the coefficients of variations of respiration rates during dry periods ${\mathrm{CV}}_{R_{\mathrm{d}}}$ (c), and on respiration pulses at rewetting ${\mathrm{CV}}_{R_{\mathrm{r}}}$ (d). The white contour curves indicate combinations of λ and α that generate different annual precipitation rates ($\langle P\rangle =\mathit{\alpha }\mathit{\lambda }=\mathrm{0.25}$, 0.5, 1, and 2 m yr⁻¹ from dotted to solid lines). Other parameter values are as in Fig. 5.

3.5 Effects of rainfall intensification and organic C availability on heterotrophic respiration

Results shown in Figs. 5 and 6 are based on average respiration model parameters; here, we explore how changing organic C content from mineral to organic soils affects the contribution of rewetting pulses to total soil heterotrophic respiration. We also focus on changes in respiration patterns along gradients of rainfall intensification, i.e., decreasing precipitation frequency λ while precipitation event depth α is increased and total precipitation is kept fixed (as along the white contour curves in Figs. 5 and 6). Figure 7 shows that rainfall intensification decreases 〈R_t〉 (Fig. 7a) but increases $\langle R_{\mathrm{r}}^{*}\rangle /\langle R_{\mathrm{t}}\rangle$ (Fig. 7b), regardless of soil organic C availability (black vs. gray curves) and total precipitation (dashed vs. solid curves). However, for a given total precipitation, organic soils (gray curves) exhibit both higher 〈R_t〉 and higher $\langle R_{\mathrm{r}}^{*}\rangle /\langle R_{\mathrm{t}}\rangle$ than mineral soils (black curves), due to their higher R_r,max (Table 2). As a result, in organic soils, the contribution of respiration pulses can be as high as 20 % of the total heterotrophic respiration, whereas in mineral soils it tends to be lower than 10 %. Moreover, in both soils, higher precipitation increases 〈R_t〉 while decreasing $\langle R_{\mathrm{r}}^{*}\rangle /\langle R_{\mathrm{t}}\rangle$ (compare solid vs. dashed curves).

Figure 7 (a) Total heterotrophic respiration 〈R_t〉 and (b) fraction of the total heterotrophic respiration rate due to rewetting pulses, $\langle R_{\mathrm{r}}^{*}\rangle /\langle R_{\mathrm{t}}\rangle$, as a function of mean precipitation event frequency λ, for given total precipitation (i.e., $\mathit{\alpha }=\langle P\rangle /\mathit{\lambda }$), and for both mineral and organic soils. The respiration model parameters are reported in Table 2, and other parameter values are as in Fig. 5.

4 Discussion

Heterotrophic respiration fluctuates at multiple temporal scales in response to hydroclimatic variability (Messori et al., 2019; Rubio and Detto, 2017) – from interannual variations due to climatic anomalies and extreme events (Reichstein et al., 2013), to seasonal variations partly linked to plant activity (Zhang et al., 2018), to short-term fluctuations induced by soil drying and rewetting (Daly et al., 2009). Here we focus on respiration fluctuations during drying–wetting cycles and how they are affected by precipitation regimes. Differently from most other modeling approaches to describe these dynamics, we develop a probabilistic model with analytical solutions for the probability density function of respiration rate (discussed in Sect. 4.1). For the sake of analytical tractability, this model rests on important assumptions (Sect. 4.2), but despite its simplicity it has the potential to assess the effect of precipitation variability (and its expected changes) on heterotrophic respiration (Sect. 4.3).

4.1 Comparison with previous stochastic approaches

Most biogeochemical models assume that heterotrophic respiration (and other processes) depend on a generic soil property φ following an empirical function f(φ) (Bauer et al., 2008; Moyano et al., 2013). As φ changes through time (e.g., soil moisture and temperature), the biogeochemical rate associated with φ also varies. Thus, the biogeochemical models use the function f(φ) to convert measured time series of soil moisture and other environmental variables into biogeochemical rates. The different approach we follow here consists in linking a known probability density of φ to the probability density of the function f(φ) to capture the propagation of the statistical properties of φ to f(φ). This can be done by the derived distribution approach, as in Eq. (8). This approach has been used to investigate gaseous nitrogen emissions in response to soil moisture fluctuations (Ridolfi et al., 2003), but the only example studying soil heterotrophic respiration rate we are aware of focused on respiration responses to temperature fluctuations (Sierra et al., 2011). These approaches provide simple and mathematically elegant solutions but have so far been limited to the effect of a single driver of the biogeochemical flux of interest. The responses of heterotrophic respiration to changes in soil moisture are more complex because rewetting pulses depend on both soil moisture increment and pre-wetting soil moisture (Fig. 2), requiring the solution of a bivariate stochastic process. Thus, our approach – by accounting for both these effects – is more general and applicable along gradients where the statistical properties describing the precipitation regime vary significantly (Fig. 4).

A previous stochastic approach focused on the CO₂ concentration in the pore space instead of respiration rates (Daly et al., 2008). Observations show that CO₂ concentration increases rapidly after rainfall and then decreases following a negative exponential function. This dynamic can be described as a stochastic process where CO₂ concentration is the random variable and precipitation represents the stochastic forcing (Daly et al., 2008). With this approach, the long-term mean CO₂ concentration was found to depend on the average rainfall rate (λα), while the standard deviation of CO₂ concentration depends on λα². This indicates that rainfall intensity (in terms of mean event depth α) plays a more important role than rainfall frequency in driving the variability of soil CO₂ concentration. Soil respiration was shown to be approximately proportional to CO₂ concentration in the pore space over a broad range of concentrations (Daly et al., 2008), so that respiration statistics are also expected to scale with rainfall statistics in the same way as soil CO₂ concentrations. This result is consistent with our finding that all components of heterotrophic respiration increase with both λ and α (Fig. 5).

Numerical process-based models have also been driven by randomly generated rainfall time series (e.g., Tang et al., 2019). These models do not allow analytical solutions for the respiration statistical properties to be found, but they offer insights into the individual processes affecting these properties. For scenarios of constant total rainfall and variable rain event frequency, Tang et al. (2019) found that rainfall intensification increased heterotrophic respiration in a semi-arid grassland, even though in their simulations soil organic C stocks also slightly increased due to higher plant productivity. This result differs from our finding that total heterotrophic respiration decreases with rainfall intensification (moving right to left along the curves in Fig. 7a) and was likely caused by how plant productivity and its feedback to soil organic C were modeled in their study.

4.2 Methodological limitations

Three model assumptions can alter the interpretation of our results: (i) that heterotrophic respiration pulses can be regarded as instantaneous, (ii) that the two parameters R_d,max and R_r,max are independent of climatic and vegetation conditions, and (iii) that hydroclimatic conditions are statistically stationary.

Respiration pulses are modeled as instantaneous events of CO₂ emission with a given size (Sect. 2.1.3). While mathematically convenient, rewetting respiration pulses are known to last for a few days after the rewetting has ended. Indeed, when analyzing laboratory incubation data, the pulse size is generally calculated by integrating through time the respiration rates above the rate occurring at stable soil moisture. The integration window ranges between 2 and 3 d (e.g., Fischer, 2009). This simplified approach to separate the actual rewetting pulse from the respiration rate at stable soil moisture requires some caution when rainfall events are frequent. In that case, pulses would overlap rather than being distinct. Moreover, with frequent rainfall, respiration could be inhibited due to water logging (Moyano et al., 2013; Rubio and Detto, 2017), and no respiration pulse might occur. Thus, to avoid these issues, our equations should not be used in wet environments with λ>0.3 d⁻¹.

We calculated the statistical properties of the heterotrophic respiration rate, but we did not consider the dynamics of the soil organic matter and plants that supply resources for microbial growth and respiration. Widely different precipitation amounts and distributions such as those depicted in Figs. 5 and 6 are associated with different plant communities, whose productivity increases along gradients of precipitation (Huxman et al., 2004; Luyssaert et al., 2007), providing litter and root exudates whose C is eventually stabilized into soil organic matter. Indeed, soil organic C stocks increase with increasing mean annual precipitation (Guo et al., 2006). Hence, soil organic matter probably varies along the axes of Figs. 5 and 6, which are instead interpreted here as purely climatic gradients. Such variations in organic matter content would affect the maximum respiration rate and pulse size, R_d,max and R_r,max (e.g., compare mineral and organic soils in Table 2). Because the mean respiration rates scale with the maximum rates (as apparent analytically from Eq. 21), it is reasonable to expect that higher organic matter content along precipitation gradients increases the sensitivity of respiration to changes in precipitation compared to predictions in Fig. 5. Indeed, even when keeping precipitation constant while varying the frequency and depth of precipitation events, the variations in total heterotrophic respiration are larger in organic soils than in mineral soils (Fig. 7a).

Moreover, soil C substrates might be depleted through multiple drying and rewetting events – a behavior we do not consider in the proposed statistically stationary model. While some experiments show sustained rewetting pulses (Miller et al., 2005; Xiang et al., 2008), others show reduced total heterotrophic respiration with increasing frequency of drying and rewetting, possibly due to substrate depletion (Shi and Marschner, 2014). To capture these dynamics, a more complex model describing the changes in substrate and microbial compartments would be needed (e.g., Brangarí et al., 2018; Lawrence et al., 2009; Tang et al., 2019) at the cost of losing the analytical tractability.

Our focus in this contribution is on heterotrophic respiration, but the data we used to parameterize the model are from laboratory studies without plants. Therefore, our heterotrophic respiration estimates neglect contributions from fresh C inputs from roots to the rhizosphere (Finzi et al., 2015; Kuzyakov and Gavrichkova, 2010). However, the timing of rhizodeposition depends on plant activity, which in turn depends on previous environmental conditions – differently from soil microbes that respond to soil moisture changes rapidly, plant responses integrate previous conditions, thereby partly decoupling root activity from current soil moisture. It is thus nontrivial to include rhizosphere processes in the current framework.

In addition to these limitations, our results should also be interpreted with caution when rainfall seasonality is important, because the assumption of stochastic stationarity (Sect. 2.1.1) may not be met, requiring the derivation of a different probability density function of soil moisture (e.g., Vico et al., 2017). Nevertheless, our results will still hold for parts of the year when the rainfall regime is relatively stable.

4.3 How are the statistical properties of heterotrophic respiration varying with changing precipitation regimes?

The axes of Figs. 5 and 6 can be interpreted in terms of changes in precipitation patterns caused by ongoing climatic changes. If rainfall in a semiarid or mesic environment increases (due to either more frequent or larger events), heterotrophic respiration also increases (Yan et al., 2014; Zhang et al., 2019) – this is not surprising as soils become on average wetter, removing water limitation and promoting microbial activity. These observations are consistent with our findings that the mean respiration pulse at rewetting and respiration during drying increases with increasing α, λ, or their product – i.e., total precipitation. However, the variability in respiration does not always change monotonically with increasing rainfall. Figure 6b shows that the standard deviation of the respiration pulses increases with more intense (higher α) and less frequent (lower λ) rainfall. In contrast, the standard deviation of the respiration rate during drying, R_d, peaks at intermediate α and λ and declines thereafter because the respiration response is flat and thus has higher variance at intermediate wetness (Eq. 7; Fig. 6a). Therefore, higher precipitation as driven by increasing α or λ is expected to increase the respiration pulses (Fig. 5b) and their variability (Fig. 6b), while decreasing their contribution to the total heterotrophic respiration (Fig. 5d).

It is perhaps more interesting to understand respiration responses to changes in rainfall patterns for given total rainfall amounts. When α and λ are changed simultaneously while keeping their product fixed (moving along the white curves in Figs. 5–6; or along the x axis in Fig. 7), the mean respiration pulse at rewetting and the standard deviations of both respiration components increase with more intermittent and intense rainfall events. In experimental rainfall manipulations that mimic the predicted climatic changes, increased variability in soil moisture associated with more intense but less frequent precipitation events decreases total soil respiration (Harper et al., 2005). This observation is consistent with our result that the mean total heterotrophic respiration decreases with rainfall intensification while maintaining a given mean precipitation rate (i.e., moving right to left along the curves in Fig. 7a). Our result is explained by the higher runoff and deep percolation losses predicted by the soil hydrologic model when precipitation events are large but rare (Rodriguez-Iturbe and Porporato, 2004). These water losses cause soil moisture to be on average lower as the precipitation regime becomes more intermittent – a pattern also confirmed empirically in rainfall manipulation experiments (Harper et al., 2005). Our approach neglects the lower plant C inputs and contributions to total soil respiration under a more intermittent precipitation regime (Harper et al., 2005), which further reduces the total (combined autotrophic and heterotrophic) soil respiration rate.

We also found that the contribution of rewetting pulses to the total heterotrophic respiration increases when rainfall becomes more intermittent and rainfall events larger (i.e., moving right to left along the curves in Fig. 7b). This result is consistent with observations in a temperate steppe (Yan et al., 2014). The rewetting pulse contribution is also larger in organic soils compared to mineral soils (gray vs. black curves in Fig. 7) – this effect is expected, because more C can be mobilized by drying and rewetting cycles in C-rich soils (Canarini et al., 2017). We can thus surmise that climatic changes causing longer dry period and more intense rainfall events (IPCC, 2012) will increase the role of pulse responses, including not only respiration but also nitrogen mineralization pulses that could release nitrogen at a time when plant uptake is low. In turn, this can cause a decoupling of nitrogen supply and demand, with possible negative consequences for ecosystem productivity (Augustine and McNaughton, 2004; Dijkstra et al., 2012).

Our findings are based on time-invariant relations between heterotrophic respiration and soil moisture, but temperature and other environmental conditions also affect microbial activity – in part directly and in part indirectly via rhizodeposition – raising the question of how our results could be impacted by other respiration-controlling factors. As a first approximation, temperature could be assumed to alter directly both respiration rates during drying and respiration pulses in a similar way. This implies that our results would hold even under fluctuating temperatures, at least during the growing season, when temperature variations are limited and precipitation can be described by a simple marked Poisson process (Sect. 2.1.1). However, a different modeling approach would be needed to quantify the mean heterotrophic respiration rate during seasons with frequent rainfall events, when respiration pulses are likely to be less important and anaerobic conditions (here neglected) could play a role. As the timescale expands from the growing season to the whole year, seasonal fluctuations in plant activity that delay the supply of C substrates to microbes will also play a role (Finzi et al., 2015; Kuzyakov and Gavrichkova, 2010), leading to a hierarchy of responses at multiple timescales – a more complex problem than the one addressed in this contribution.

5 Conclusions

Heterotrophic respiration depends nonlinearly on soil moisture – not only does it follow soil moisture during a dry period, but it also responds rapidly to rewetting. These rewetting responses occur in the form of pulses of CO₂ whose size increases with increasing soil moisture increment and decreasing pre-wetting soil moisture. We used this relation between respiration pulses and soil moisture to analytically characterize the statistical properties of respiration rates as a function of the statistical properties of the rainfall events that drive soil moisture changes. Consistent with empirical evidence, our model predicts that dryer climatic conditions (either lower rainfall depths or longer dry periods between two rain events) lower total heterotrophic respiration. More interestingly, we showed that the contribution of rewetting pulses to the total heterotrophic respiration increases in dryer climates, but also when the precipitation regimes shift towards more intermittent and intense events (even at constant total average rainfall). Therefore, our results suggest that the expected intensification of precipitation will increase the role of rewetting respiration pulses in the ecosystem C budgets.