Stochastic Modeling of Soil Salinity

Large areas of cultivated land worldwide are affected by soil salinity. [Szabolcs (1989)] estimates that 10% of arable land in over 100 countries, and nine million km² are salt affected, especially in arid and semi-arid regions (Tanji, 1989). Salinity refers to large concentrations of easily soluble salts present in water and soil on a unit volume or weight basis (typically expressed as electrical conductivity (EC) of the soil moisture in dS/m, i.e. deciSiemens per meter at $25^{\circ}$ C; for $NaCl$ 1 mg/l $\sim$ $15\cdot 10^{-4}$ dS/m). High salinity causes both ion specific and osmotic stress effects, with important consequences for plant production and quality. Normally, yields of most crops are not significantly affected if EC ranges from $0$ to $2$ dS/m, while above levels of $8$ dS/m most crops show severe yield reductions (Ayars et al., 1993; Hillel, 2000). Prevention or remediation of soil salinity is usually done by leaching salts, and has resulted in the concept of leaching requirement (Richards, 1954; Hillel, 1998; Schleiff, 2008). Alternative amelioration strategies by harvesting salt-accumulating plants appear to be less effective (Qadir et al., 2000).

Salt accumulation in the root zone may be due to natural factors (primary salinization) or due to irrigation (secondary salinization). Several detailed numerical models have been developed to model soil salinization (e.g. Eldin et al., 1987; Schoups et al., 2006; Corwin et al., 2007). Generally, these models simulate unsaturated soil water flow via Richards and solute transport equations. These models are more suitable for local and short-term simulations, as they require precise soil characterization and are computationally demanding. Moreover, it is often difficult to identify cause-effect relationships or to synthetically compare the effects of different parameter scenarios from their numerical simulations.

Vertically-averaged soil moisture and salt balance equations have also been used (Allison et al., 1994; Hillel, 2000). Despite their simplicity, these models have the advantage of parsimony, thus allowing a direct analysis of the interplay of the main processes, and provide an ideal starting point to include external, random hydroclimatic fluctuations in the analysis of long-term salinization trends. The goal of this Letter is to offer a first step in this direction. With this purpose, here we present a minimalist model of soil primary salinization, describing analytically the long-term dynamics of salt in soils caused by wet (rain) and dry (aerosol) deposition. Our aim is to quantify the salt mass and concentration probability density functions (pdfs) in the root zone, and the probability of crossing the crops salt tolerance threshold as a function of the main hydro-climatic parameters. The model framework is potentially extendible to systems including salt input from groundwater and irrigation.

Our starting point is a spatially lumped model (Bras and Seo, 1987) for the vertically averaged dynamics of soil moisture and salt in the root zone. As a first step we will not consider input of salt due by irrigation or groundwater upflow. Following Rodriguez-Iturbe et al. (1999), Rodriguez-Iturbe and Porporato (2004) and Porporato et al. (2004), rainfall ($R_{t}$) is modeled as a marked point process with frequency $\lambda_{{}_{P}}$ and with daily rainfall depths exponentially distributed with mean $1/\gamma_{{}_{P}}$. The averaged soil moisture dynamics are modeled assuming constant (spatially and temporally averaged) soil and ecohydrological parameters, i.e., root depth, $Z_{r}$, porosity, $n$, and maximum evapotranspiration rate, $ET_{max}$. Assuming a rain salt concentration $C_{R}$ and a constant input $\mathcal{M}_{d}$ of salt mass per unit ground area and per unit time by dry deposition, the root-zone mass balance for soil moisture and salt mass $m$ is given by:

where $C$ is the salt concentration in the root zone; $L(s)$ represents deep percolation, while $ET(s)$ represents the losses resulting from plant transpiration and soil evaporation. As in Porporato et al. (2004), $ET(s)$ is assumed to be linear in the range of soil moisture comprised between the wilting point, $s_{w}$, and a suitable soil moisture threshold $s_{1}$ (an effective field-capacity threshold), at which $ET$ occurs at the maximum rate $ET_{max}$. All the rainfall input that cannot be accommodated is assumed to be lost as $L(s)$ at $s_{1}$. In this minimalist model the effect of salt-induced changes in osmotic potential may only indirectly be taken into account through an average reduction of $ET_{max}$. This is simply done here by keeping the same $ET_{max}$ for the two models (previous studies (Viola et al., 2008) have shown that, in the absence of osmotic effects, the minimalist model should have artificially higher $ET_{max}$ to account for percolation losses below $s_{1}$). A full account of how reduction in evapotranspiration affects salinization patterns (reduced evapotranspiration in turn increases the available soil moisture and thus reduces the concentration of salt in the soil and increases leaching frequencies) will be given elsewhere.

A complete numerical model, in which the impact of osmotic stress in reducing $ET$ is explicitly included (Bras and Seo, 1987), has been also studied. Moreover, in the detailed model runoff takes place at saturation ($s=1$), while percolation occurs for $s>s_{fc}$ (the soil moisture field capacity), and it is proportional to the soil hydraulic conductivity $K_{sat}s^{c}$, where $c$ is a soil-pore connectivity index and $K_{sat}$ is the saturated hydraulic conductivity (Rodriguez-Iturbe and Porporato, 2004). A comparison between the results of the two soil moisture models, presented in Figure 1a, suggests the viability of the simplified model. Simulations for wetter climates confirm this result.

The system (1) and (2) can be further simplified if one considers that the typical timescales for salt mass dynamics in the root zone are orders of magnitude larger than the ones characterizing rainfall (and thus wet deposition). Moreover, soil moisture typically reaches steady-state conditions within a growing season (e.g., $<5-7$ months), while the salt mass balance only does so on much longer times scales (e.g., $>$ decades). Accordingly, at those long timescales, say $T$, the salt mass input flux can be assumed to take place at a constant rate, $\Upsilon$, that is $\int_{t}^{t+T}(\mathcal{M}_{d}+C_{R}R_{t})dt^{\prime}\sim\mathcal{M}_{d}T+T\,C_{R}\lambda_{{}_{P}}/\gamma_{{}_{P}}=\Upsilon T$, and be interrupted by instantaneous and unfrequent leaching events induced by percolation. As a result, (2) can be rewritten as

Leakage may then be modelled as a marked point process, with percolation depths exponentially distributed with parameter $\gamma_{{}_{P}}$ (Botter et al., 2007). For reasons of analytical tractability, the percolation events are assumed to occur according to a Poisson process with frequency $\lambda$ given by the frequency of soil moisture crossing the threshold $s=s_{1}$. This can be expressed in terms of the soil moisture pdf as $\lambda=\rho(s_{1})p(s_{1})$, where the term $\rho(s)=(ET(s)+L(s))/nZ_{r}$ represents the normalized catchment-scale loss function (i.e. the total losses from the system due to evapotranspiration and leakage as a function of the soil moisture) (Rodriguez-Iturbe and Porporato, 2004). Adopting the soil moisture minimalist model (for which the pdf is a truncated gamma distribution, e.g., Porporato et al. (2004)), the leaching frequency is $\lambda=\eta\exp(-\gamma)\gamma^{\lambda_{{}_{P}}/\eta}/\Gamma(\lambda_{{}_{P}}/\eta,\gamma)$ (Botter et al., 2007), where $\Gamma(x,y)$ is the lower incomplete gamma function, $\eta=ET_{max}/(nZ_{r}(s_{1}-s_{w}))$ and $\gamma=\gamma_{{}_{P}}nZ_{r}(s_{1}-s_{w})$.

A leaching-efficiency parameter $b$ is used to account for incomplete salt dissolution, further assuming that the typical value of soil moisture during leaching events can be approximated by the value $s_{1}$. With the above assumptions, the dynamics of the salt mass in the root zone can be described by a single equation

where $L^{\prime}_{t}$ is a marked Poisson noise (Van Den Broeck, 1983) with frequency $\lambda$, and (dimensionless) exponential marks with mean

Figures 1b and 1c compare the results of both salinity models. The free parameters $s_{1}$ and $b$ are fitted with respect to the complete model of salt mass and concentration, respectively.

From a mathematical viewpoint, equation (4) is a stochastic differential equation with multiplicative white (jump) noise. In our case, since the soil solution can be considered in equilibrium during leaching events, one has to interpret (4) in the Stratonovich sense (Van Den Broeck, 1983). Accordingly, the normal rules of calculus are preserved, and equation (4) can be transformed into

where $y(t)=\ln[m(t)]$.

The stationary solution of (6) can be obtain as in Rodriguez-Iturbe et al. (1999). Then using the derived distribution for $m$, i.e., $p(m)=p(y)dy/dm$, we obtain the probability distribution for the salt mass in the root zone

where $\mathcal{N}=(\lambda/\Upsilon)^{\frac{1+\mu}{\mu}}/\Gamma(\frac{1+\mu}{\mu})$ and $\Gamma(x)$ is the Gamma function. Equation (7) summarizes the soil salinity statistics as a function of climate, soil and vegetation parameters.

Figure 2 is a graphical representation of the dependence of the mean salt concentration $\langle C\rangle=\langle m\rangle/nZ_{r}\langle s\rangle$on the yearly rainfall and $\lambda_{{}_{P}}$. The contour-lines connect equal values of the mean salt concentration in the soil, for a given input of salt $\Upsilon$. The latter has been calculated for two different geographic regions. Typical salt inputs in coastal areas are $100-200$ kg/(ha yr) of salt, while values drop of an order of magnitude in continental regions (Hillel, 2000).

Between the black region and the light gray ones in Figure 2a, the behavior of $\langle C\rangle$ changes substantially. Above a certain total rainfall per year, the input of salt related to rainfall frequency becomes immaterial as leaching effectively washes out the salt mass from the root zone. For lower total rainfall values, however, the salt in the soil increases with increasing $\lambda_{{}_{P}}$. For a given annual precipitation depth, with low rainfall frequencies, rainfall events carry enough water to trigger leaching. Conversely, if $\lambda_{{}_{P}}$ is high, evapotranspiration dominates, leaching is largely reduced, thereby causing salt accumulation in the root zone. Therefore, $\langle m\rangle$ strongly increases with $\lambda_{{}_{P}}$. Relatively small reductions of rainfall at the transition between these two regimes may entail a dramatic increase in long-term soil salinization. Figure 2 also shows the threshold of soil salinity below which vegetation is practically unaffected (e.g., $\langle C\rangle<2$ dS/m) and the thresholds above which regular (e.g. non-halophytic) vegetation is damaged (e.g., $\langle C\rangle>2$ dS/m). For coastal areas soil salinization may occur even in relatively more humid regions, especially when rainfall events are not very intense. On the contrary, in continental regions only arid climates may begin to develop soil salinization (in the absence of irrigation and groundwater input). Indeed, through our model one can evaluate the risk of soil salinization in rain-fed agriculture just by estimating the typical salt inputs, total rainfall per year and the rainfall frequency. For example, a rain-fed crop in a semi-arid climate (e.g., rainfall depth of $70$ cm/yr) in a continental region risks salinization only when rainfall events are not very intense (e.g., $\gamma_{{}_{P}}^{-1}\leq 0.4$ cm or $\lambda_{{}_{P}}\geq 0.48$ d^-1 ). If the same crop is located in a coastal area, salinization occurs for a wider range of rainfall parameters (e.g., $\gamma_{{}_{P}}^{-1}\leq 1$ cm or $\lambda_{{}_{P}}\geq 0.18$ d^-1 ).

The solution (7) may be used in conjunction with soil moisture statistics to obtain a full characterization of the salt concentration in the root zone. Because one may safely assume that equations (1) and (3) are decoupled over short time scales, the soil moisture $s(t)$ and the salt mass $m(t)$ may be treated as statistically independent random variables. By observing that the salt concentration in the root zone is equal to $C(t)=m(t)/nZ_{r}s(t)$ and assuming $s_{w}\sim 0$, we find the stationary probability distribution of the salt concentration $p(C)$ as the quotient distribution of two independent random variables (Curtiss, 1941),

The comparison between analytical solutions and numerical simulations (Figure 3) shows that the analytical solution reproduces reasonably well the pdf of the complete model.

By integrating equation (8) from a given concentration value $C^{*}$ to infinity, one obtains the cumulative pdf of $C$, $P(C^{*})$, which is the probability of having a salt concentration greater than a certain critical concentration value, $C^{*}$, as a function of the soil-plant-atmosphere parameters. The inset of Figure 3 confirms the impact that climate change may have on soil salinity. Note, in particular, that such an impact is marked only for semi-arid or drier climates (see Figure 2). For example with a reduction from $\lambda_{{}_{P}}=0.2$ to $\lambda_{{}_{P}}=0.15\,$d^-1, the probability of crossing $C^{*}=6$ dS/m is more than tripled. When coupled to a crossing analysis of concentration levels, the previous results may be used to evaluate the risk of plant salt stress. The analytical form of the results makes it suitable for computations of salinity risk at the global scale as a function of few measurable parameters, and facilitates their coupling with other models of long-term soil-plant biogeochemistry.

In this Letter we have presented an analytical approach to stochastic modeling of soil salinity, where the complexity of the problem is reduced by employing simplifying assumptions that permit us to describe high-dimensional, unpredictable components via suitable random terms. By assuming time-averaged inputs of salt and instantaneous percolation processes, a decoupling from soil moisture equation results in a simplified stochastic mass balance equation for the soil salt mass amenable to exact solution.

Soil salinity statistics are obtained as a function of climate, soil and vegetation parameters. These can be combined with soil moisture statistics to obtain a full characterization of soil salt concentrations and the ensuing risk of primary salinization.

This modeling framework can be extended to investigate additional salt inputs from irrigation and groundwater by modifying accordingly the average salt input parameter $\Upsilon$ and calculating the corresponding soil moisture pdfs (e.g. see Vervoort and Van Der Zee (2008) for groundwater inputs and Vico and Porporato (2009) for irrigation).