Hysteresis in soil hydraulic conductivity as driven by salinity and sodicity – a modeling framework
Declining soil-saturated hydraulic conductivity (Ks) as a result of saline and sodic irrigation water is a major cause of soil degradation. While it is understood that the mechanisms that lead to degradation can cause irreversible changes in Ks, existing models do not account for hysteresis between the degradation and rehabilitation processes. We develop the first model for the effect of saline and sodic water on Ks that explicitly includes hysteresis. As such, the idea that a soil's history of degradation and rehabilitation determines its future Ks lies at the center of this model. By means of a “weight” function, the model accounts for soil-specific differences, such as clay content. The weight function also determines the form of the hysteresis curves, which are not restricted to a single shape, as in some existing models for irreversible soil processes. The concept of the weight function is used to develop a reversibility index, which allows for the quantitative comparison of different soils and their susceptibility to irreversible degradation. We discuss the experimental setup required to find a soil's weight function and show how the weight function determines the degree to which Ks is reversible for a given soil. We demonstrate the feasibility of this procedure by presenting experimental results showcasing the presence of hysteresis in soil Ks and using these results to calculate a weight function. Past experiments and models on the decline of Ks due to salinity and sodicity focus on degradation alone, ignoring any characterization of the degree to which declines in Ks are reversible. Our model and experimental results emphasize the need to measure “reversal curves”, which are obtained from rehabilitation measurements following mild declines in Ks. The developed model has the potential to significantly improve our ability to assess the risk of soil degradation by allowing for the consideration of how the accumulation of small degradation events can cause significant land degradation.
Soil degradation caused by the use of saline and sodic irrigation water is a major threat to agricultural production and food security. Each year, as much as 1.5 million ha of farmland are irretrievably lost due to salinity- and sodicity-induced degradation, while crop output is decreased on an additional 20 to 46 million ha (FAO and ITPS, 2015). Together, this results in widespread economic losses and declines in production that elevate stress on food supplies (Qadir et al., 2014). Growing strain on freshwater resources, rising human populations, and climate change are all likely to increase reliance on saline and sodic irrigation water to support food production, further raising the risk of debilitating soil degradation in the near future.
In this paper, we focus on changes to saturated hydraulic conductivity, Ks, as the major indicator of soil degradation (McNeal and Coleman, 1966; Yaron and Thomas, 1968; Chen and Banin, 1975; Shainberg and Singer, 2011; Hillel, 2000). Saline and sodic conditions cause declines in Ks through a combination of different soil processes, including clay swelling, clay dispersion, and slaking. While swelling is usually considered to be a reversible process, clay dispersion is thought to cause irreversible changes in the soil matrix and, thus, Ks (Quirk and Schofield, 1955; McNeal and Coleman, 1966; Chen and Banin, 1975; Ezlit et al., 2013; Dang et al., 2018a, b). No clear dividing line exists, however, between the onset of swelling and dispersion, and experimental studies have shown that declines in Ks are characterized by smooth curves as water quality is changed, rather than by abrupt jumps (Quirk and Schofield, 1955; McNeal and Coleman, 1966; Chen and Banin, 1975). We hypothesize that this smooth behavior is due to soil heterogeneity, with different aggregates exhibiting different thresholds between swelling and dispersion, and the overall behavior of the soil corresponding to the integrated response of all aggregates. This point constitutes a cornerstone of the framework upon which our model lies, and we will return to it during the model development.
While numerous experimental studies have led to models describing declines in Ks due to salinity and sodicity, we are unaware of any model that considers differences between degradation and rehabilitation (McNeal and Coleman, 1966; McNeal, 1968; Yaron and Thomas, 1968; Lagerwerff et al., 1969; Russo, 1988; Ezlit et al., 2013; Ali et al., 2019). That is, all existing models assume that these two processes are strictly reversible, meaning decreases in Ks can be restored simply by changing the chemical properties of the soil water (McNeal and Coleman, 1966; McNeal, 1968; Yaron and Thomas, 1968; Lagerwerff et al., 1969; Russo, 1988; Ezlit et al., 2013). There is strong reason to believe, however, that this is not the case. Clay dispersion, for instance, is widely understood to be irreversible (Goldenberg et al., 1983; Yaron et al., 2008). For this reason, a significant reduction in Ks, i.e., a decrease of more than 15 %–25 % compared to its initial level, is considered a threshold beyond which a soil will suffer irreparable damage (Quirk and Schofield, 1955; McNeal and Coleman, 1966; Rengasamy et al., 1984; Cook et al., 2006; Bennett et al., 2019). Likewise, the experimental evidence that has considered degradation and rehabilitation together suggests that changes in Ks are defined by hysteresis (Dane and Klute, 1977), the exact nature of which can be expected to vary on a soil-specific basis (Levy et al., 2005).
Failure to account for hysteresis limits the ability to forecast whether a particular climate and irrigation regime will lead to long-term land degradation (Kramer and Mau, 2020). When even a crude hysteresis mechanism is considered, it has been shown that changes in Ks may be compounded, such that initially small declines in Ks may lead to serious degradation (Kramer and Mau, 2020). Considering the limited nature of land and fresh water resources, maximizing land production through the use of saline and sodic water is likely to prove essential for maintaining food security. It is similarly critical, however, to realize this goal without causing long-term damage to soils. The long-term viability of irrigation practices with poor-quality water therefore hinges on improving our understanding of how hysteresis affects the processes of degradation and rehabilitation.
In this paper, we present a model for the effect of salinity and sodicity on Ks that considers hysteresis. That is, in contrast to existing models, the presented model accounts for how a soil's history of degradation and rehabilitation will affect Ks in the future. Our goal is to provide a modeling framework that can be incorporated into numerical models of water and solute dynamics to tackle questions of sustainable water resources use. As such, the presented model represents an important step in bringing forward the ability to assess the risk of soil degradation and to identify optimal rehabilitation strategies.
In constructing our model, we rely on the Preisach framework, which was originally developed in magnetism and has been applied to a wide variety of fields, including the modeling of soil water retention (Flynn and Rasskazov, 2005; Flynn et al., 2006; O'Kane and Flynn, 2007; McNamara, 2014) and the effect of stress and strains on soil and rock structures (Guyer, 2006). A comprehensive overview of the Preisach model can be found in Mayergoyz (1986, 2003). We begin this paper by introducing the Preisach framework and how it can be applied to model hysteresis in soil Ks. We then discuss the experiments necessary to parameterize the model and to account for soil-specific properties. We demonstrate a limited version of this experimental procedure and show how the parameterized model can be used to infer future behavior. Finally, we incorporate data from past degradation experiments to demonstrate how our model can be used to improve future assessments of the risk of soil degradation.
Our modeling of hysteresis in Ks is based on Preisach's framework. We will briefly present this framework and then apply it to model Ks. We highly encourage the reader to use the interactive widgets found in the GitHub repository for this paper (Kramer et al., 2020), as they greatly enhance the explanatory power of the figures in the following sections.
2.1 The Preisach framework
At the foundation of the Preisach framework is the idea of the hysteron, which functions as a system's basic memory unit. The hysteron admits two output states, “on” (f=1) and “off” (f=0), that are determined by an input, u(t), and its history. Figure 1a presents a single hysteron and the hysteresis loop between two threshold values, α, β, as follows:
where k is the previous state of f and, by convention, α≥β.
The conceptual leap offered by the Preisach framework is to imagine a myriad of hysterons in all possible threshold configurations, each infinitesimally contributing to the system's overall behavior. This step recognizes the heterogeneity within a system, where its different constituents may not respond in tandem to external forcing, which we hypothesize is the case in soils. To this end, the Preisach model is often interpreted geometrically over a limiting triangle (α≥β), in which each unique pair of points (αi, βi) corresponds to a hysteron, as depicted in Fig. 1b. The system's output can then be determined by integrating over the limiting triangle, considering the fraction of hysterons that are on. Increases in input are modeled by turning on hysterons for which u(t)≥α. Decreases in input are modeled by turning off hysterons for which u(t)≤β (see Widget 1 at https://github.com/yairmau/hysteresis-python, last access: 8 April 2021 on GitHub). For simplicity's sake, we will use the more natural ordering (α, β), although in all the graphs in this paper, α is the ordinate (vertical) axis and β is the abscissa (horizontal) axis.
The geometric interpretation of the Preisach framework is demonstrated in Fig. 2, where the panels represent the following steps:
We begin at t=0 with input u=u0. All of the hysterons within the limiting triangle are on, and output, f, is at its maximum level. This state is known as positive saturation.
Input u is decreased to u1. We swipe left on the limiting triangle, turning off all hysterons for which β>u1, causing a decrease in f.
Input is increased to u2. We swipe up on the limiting triangle, turning on all hysterons for which α<u2, leading to an increase in f.
A further decrease in input to u3 will again cause a swipe left, turning off all hysterons for which β>u3, leading to a decline in f.
2.2 Application of the Preisach framework to soil hydraulic conductivity
The Preisach model can be adapted to study the effects of salinity and sodicity on soils by setting Ks as the output, while the inputs are given by the chemical composition of the soil water. In our model, the two inputs are electrolyte concentration (C; mmolc L−1) and sodium adsorption ratio (SAR; mmol L). In doing so, the framework presented in Sect. 2.1 must be expanded to accommodate two input variables, but this is easily accomplished (see Sect. A). In addition to the two input variables (C and SAR), many covariates can have an important role in the Ks response, such as clay content, clay mineralogy, soil pH, organic matter, etc. The Preisach framework is able to account for soil-specific characteristics by using a weighting function, μ.
The weighting function μ(α,β) assigns a relative weight to each hysteron within the limiting triangle. When integrating over the limiting triangle to compute the output, the value of each hysteron at (αi, βi) is then multiplied by its respective weight μ(αi,βi). An interactive example of the role of the weighting function can be found in Widget 2 at https://github.com/yairmau/hysteresis-python (last access: 8 April 2021) in our GitHub repository.
The distribution of weights significantly affects the shape of the hysteresis curves and, therefore, the output value, given an input history. Assuming a constant SAR so that salinity, C, is the only input variable, Fig. 3 presents three weight functions together with their corresponding hysteresis curves. Every weight function (left-hand column) is associated with a unique set of hysteresis curves showing how Ks responds to changes in C (right-hand column). The hysteresis graphs depict major loop curves (solid lines) and first-order reversal curves (hereafter FORCs, indicated by dashed lines). The major hysteresis loop determines the limits of the system. The upper half of the major loop defines the resilience of an initially stable soil to degradation, while the lower half indicates the path to rehabilitation for a completely degraded soil. The smaller the gap between the upper and lower curves, the greater the degree of reversibility if Ks is fully degraded. That is, when the gap is thin, as in Fig. 3a, degradation in Ks is highly reversible. When the gap between the two curves of the major hysteresis loop is wider, as in Fig. 3c, degradation in Ks is less reversible. The FORCs correspond to the system's behavior within the limits of the major loop, showing the degree of hysteresis when degradation in Ks is limited. The closer the FORCs are to the top curve of the major loop, the higher the degree of reversibility. The FORCs in Fig. 3 show high, medium, and low degrees of reversibility, for Fig. 3a–c, respectively. Because the weight function uniquely defines a hysteresis graph, characterizing the degree to which a given soil is susceptible to irreversible changes in Ks can be done by finding that soil's weight function μ(α,β).
2.3 Reversibility index
To facilitate quantitative comparison of the degree to which different soils are reversible, we use the weight function to define a reversibility index, Ri. The full development of the reversibility index is contained in Sect. C, but the fundamental idea behind it is that the closer the weight distribution is to the diagonal of the limiting triangle (i.e., the line α=β), the more reversible the soil, as in Fig. 3a. Conversely, as the weight distribution becomes closer to the vertex of the triangle in the top left (αmax, 0), the more irreversible the soil, as in Fig. 3c. The reversibility index is defined to range from 0 to 1, with a value of one indicating a completely reversible soil (no hysteresis) and a value of zero corresponding to a completely irreversible soil.
As a result of the one-to-one correspondence between the weight function and the hysteresis graph, a soil's weight function can be determined by measuring its FORCs. In this section, we outline the experimental framework needed to measure FORCs and demonstrate how the data collected from these experiments can be used to determine a soil's weight function. While the experimental work necessary to compute weight functions for multiple soils is outside the scope of this paper, in what follows we present a “proof of concept” wherein we show limited FORC measurements and use these to construct a partial weight function. This proof of concept relies on soil column experiments, but we also note that the same framework could be used to measure FORCs in the field.
3.1 Experimental setup
Experiments to measure FORCs should begin with the same setup commonly used to study degradation of Ks due to salinity and sodicity (McNeal and Coleman, 1966; Ezlit et al., 2013; De Menezes et al., 2014). In these experiments, a fixed water head (H) is imposed on the soil column, so that Ks can be found using Darcy's law, , where q is the discharge flux, and ∇H is the hydraulic head gradient. The SAR of the input water is typically held constant, as salinity is gradually reduced, with the procedure repeated for a range of SAR values. Conversely, the role of the variables can be flipped, with C held constant and the SAR of the applied water gradually increased. In either case, changes in Ks are recorded as the quality of the input water is changed. In the context of the Preisach framework, these experiments begin from “positive saturation”, i.e., the state when all hysterons are on. This corresponds to the point where both Ks and the input variable are maximal, as seen in the annotation in Fig. 3a. As the quality of the applied water is changed (C decreased or SAR increased), the experiments effectively sample the upper curve of the major hysteresis loop.
Measurement of FORCs can be achieved by extending the procedure described above so that, following the initial decline in Ks, we “reverse” the trend in water quality (now either towards higher C or lower SAR) and measure the response of Ks to conditions that favor improvements in soil structure. Several FORCs can be measured by repeating this procedure and varying the reversal point. For example, for the case when SAR is fixed, one should alter the point of minimum salinity at which the trend in C is reversed. In Fig. 3, two FORCs (the dashed lines) are shown for each hysteresis graph. Because the presented model considers changes in Ks as a result of both C and SAR, calculation of a particular soil's complete weight function requires the measurement of two sets of FORCs, one in which changes in Ks are measured as SAR is held constant, and one in which C is held constant.
After measuring several FORCs, the weight function can be calculated by taking the mixed partial derivative of the output (i.e., Ks) at points along these curves (Mayergoyz, 2003, see Eq. A4). An annotated example of the code necessary to compute the weight function from experimental data is included at https://github.com/yairmau/hysteresis-python (last access: 8 April 2021). We note that in the case where SAR is varied while C is kept constant, positive saturation occurs for low input values and not high values, as assumed in the presented graphs. In order to keep the mathematical framework the same for all cases, we define the two inputs in the code as C and −SAR, thus maintaining uniformity; Ks declines as the inputs are decreased.
3.2 Initial results
To showcase how FORCs can be used to calculate a soil's weight function, we conducted a limited experiment to measure FORCs on a test soil. In doing this, our goal is to demonstrate both the feasibility of such experiments and the suitability of the tools we developed to analyze the results. As such, this experiment does not constitute the whole set of procedures delineated above, which are necessary for the full characterization of a soil, but must be seen as a proof of concept. To the best of our knowledge, the results shown in this subsection are the first experimental account of the partial reversibility in Ks as influenced by salinity and sodicity.
In this experiment, we used a brown steppe soil (Ravikovitch, 1992), from the Kiryat Gat region of Israel, with 25 % clay content. We held SAR equal to 20 (mmolc L−1) constant and measured Ks, as the value of C was repeatedly decreased and then increased. Measurements of Ks were recorded following the equilibration of the flow, the electrical conductivity, and the sodium concentration between the soil column's inflow and outflow solutions (in three replicates). The typical time for equilibration of a single measurement point was 3 h, after which determination of Ks was conducted, and the solution was changed for the next step of the procedure. For more details on the soil properties and experimental setup, see Sect. B.
Figure 4a shows the major loop of degradation and first-order reversal curves, starting from concentrations 5, 10, 15, 20, and 50 mmolc L−1. In what follows, we will refer to each FORC by its lowest concentration value. The degradation in Ks experienced by the soil before salinity was decreased below C equal to 20 mmolc L−1 was found to be fully reversible, as FORCs 20 and 50 depict (they overlap with the major loop). FORCs 5, 10, and 15 show partial reversibility only, with Ks not returning to its original value, even for very concentrated solutions of 400 mmolc L−1 (not shown in the graph).
We used the data from the FORCs in Fig. 4a to determine a limited weight function, shown in Fig. 4b. This is a limited weight function because it contains information on one SAR value only. Because the decay in Ks, shown by the major loop curve, is most pronounced for low C values, most of the weight is concentrated on the left of the triangle. Furthermore, the moderate rehabilitation in Ks shown by the FORCs also takes place in low C values; therefore, we see a higher density on the bottom of the triangle. Note that the triangle's upper input bounds are 300 mmolc L−1, while the experiment's upper bound was 200 mmolc L−1. This expansion of the triangle is necessary because the reversal curves do not go back to full rehabilitation (relative Ks=1.0), and according to Preisach's framework, the system always reaches full saturation upon arrival at the upper input bound. The expansion of the triangle leaves some weight beyond the experimental bounds (see the dotted line in Fig. 4b), thus allowing for incomplete rehabilitation.
In order to make sure that the weight triangle contains the information of the hysteresis graph, we modeled the soil's response using this weight function. The modeled results, shown in Fig. 4c, are quite similar to the experimental data of Fig. 4a. It is worth noting that the modeled FORCs return to higher Ks values, i.e., there is a slight overestimation of rehabilitation, and several FORCs converge as C increases. This is due to the effort to build a high-resolution picture of the weight function from very few data points. It is possible to eliminate this issue by increasing the number of FORCs (as shown in Fig. B1). Of course, this comes at the cost of additional time and resources, and we believe that our results are quite satisfactory, given the high level of upscaling involved in the process. The algorithm devised to reconstitute the weight function from experimental data (see the GitHub repository) incorporates both the qualitative and quantitative behavior of the soil under salinity and sodicity changes.
3.3 Insights for future experiments
This limited experiment yields a number of important insights, and we hope that these will prove useful to future efforts to characterize trends of partial reversibility in soil Ks.
The experimental setup is relatively simple, and the whole procedure is straightforward, although full characterization of a soil's weight function has the potential to be long and laborious. Determination of a soil's complete weight function requires measuring hysteresis graphs for at least five fixed SAR values (i.e., when C is varied, as seen in Fig. 4a) and for five fixed C values (i.e., when SAR is varied). Each of these hysteresis graphs must contain at least five FORCs, with each FORC itself requiring several measurement points. The suggested number of five is, of course, a rule of thumb. If one wishes to obtain weight functions of finer detail, then more hysteresis curves or FORCs per hysteresis could be measured. We do not recommend, however, measuring fewer than four hysteresis curves or fewer than four reversal curves, since this would translate into highly coarse weight functions, which would limit their usefulness.
Given the typical times of equilibration we experienced, we estimate that a full characterization of the soil hysteresis would take approximately 1 month. Soils with higher clay content may require significantly more time, since these soils can be expected to need additional time to equilibrate. Conversely, characterization of soils with a lower clay content may require less time than we have estimated here.
The first thing one should do before measuring the FORCs is to measure the major loop (i.e., the pattern of decay) to determine the range of input values in which Ks varies the most. The points of return of each FORC should be distributed roughly equally along the decay in Ks. For instance, Fig. 4a shows that the return points are more concentrated between concentrations 5 and 20 mmolc L−1 because these return points roughly sample the steepest decline in Ks. Conversely, return points where the slope of the major loop is flat do not yield as much useful information for characterizing the weight function. It is important to note that the measurement points along all FORCs have the same input values. This is necessary because of the way we developed our algorithm for building the weight functions.
In this section, we use the model to demonstrate how consideration of hysteresis is essential when analyzing degradation and rehabilitation in soils. While controlled soil column experiments are expected to behave differently than field conditions, they still provide invaluable information regarding degradation and rehabilitation trends. More importantly, these experiments allow us to generate testable hypotheses on the interactions between soil and environmental conditions over the long run. However, the modeling formalism presented here is general and does not depend on the idiosyncrasies of particular experiments. This is to say that, while we measure FORCs using soil columns, FORCs could also be measured with field experiments. The lab experiments described here are a first attempt to parametrize and constrain the model, and further lab and field experiments are necessary to ensure the model's applicability in real-life conditions. At this stage, no experimental work has addressed the Ks hysteresis in detail. We incorporate results from the well-known McNeal and Coleman (1966) experiments to show how it is possible for soils with the same degradation patterns to exhibit very different levels of reversibility.
The experimental setup employed by McNeal and Coleman (1966) is ideal for this purpose because it resembles the degradation portion of the experimental procedure described in Sect. 3. Specifically, McNeal and Coleman (1966) measured Ks as input water of fixed SAR, but declining salinity, was applied to soil columns. The procedure was repeated for several SAR values, as shown in Fig. 5a, which reproduces the degradation curves presented by McNeal and Coleman (1966) for SAR 15, 25, 50, and 100 mmol L (henceforth, we will omit the units of SAR).
The possibility for different levels of reversibility, given the same degradation curves, is presented in Fig. 5b and c. Here, we use the degradation curve associated with SAR 50 to construct two potential sets of FORCs. The FORCs in Fig. 5b correspond to a low level of reversibility – once the soil has been degraded, an increase to higher Ks values requires a significant increase in salinity. By contrast, the FORCs in Fig. 5c indicate a high level of reversibility, such that Ks values are more easily restored following degradation. It is important to emphasize that we constructed both sets of presented FORCs, and that neither derive from experiments. That is, they are meant to demonstrate the need for detailed experiments on irreversibility, according to the guidelines discussed in Sect. 3, which will allow us to more accurately determine the true extent to which a given soil exhibits hysteresis in degradation and rehabilitation. We also note that even the concave FORCs in Fig. 5c consider a higher level of irreversibility than existing models for degradation and rehabilitation, which allow Ks to increase by returning along the original degradation curve (i.e., the top curve of the major loop).
By extending the approach used in Fig. 5b and c, we produce two potential weight functions based on the degradation curves in McNeal and Coleman (1966). To do this, we first construct two sets of FORCs (low and high reversibility) for each of the respective SAR values (15, 25, 50, and 100). We then use these FORCs to compute weight functions, using the technique and code described in Sect. 3. The resulting weight functions are presented in Fig. 6. The different areas occupied by the limiting triangle in these plots reflect the different values of Cmax – the unique salinity value at which degradation begins – for each respective SAR value. When determining the weight function, salinity values greater than umax are not assigned any weights.
Examination of the weight functions in Fig. 6 shows the crucial way in which our model is able to distinguish between less reversible and more reversible soils. Focusing first on SAR 100, we observe that, for the weight function derived from the less reversible FORCs (Fig. 6a), the heavy weights are concentrated exclusively in the region close to (αmax, 0). For the weight function derived from the more reversible FORCs (Fig. 6b), the heavy weights are similarly concentrated in the region for which β is close to zero. A major contrast can be seen, however, in the distribution of the weights along the α axis, with the heavy weights in the more reversible case distributed along the entire α axis and not only in the region for which α is large, as is the case for the less reversible example. In both examples, the weight distributions along the β axis are identical because this axis is associated with degradation, and both weight functions were derived from the same degradation data. The concentration of weights in the region where β is low reflects the soil's relative resilience to degradation when salinity is high, as seen in Fig. 5. The α axis, by contrast, is associated with rehabilitation. Because the weights in Fig. 6b are distributed along the entire α axis, we can expect any increase in salinity to lead to an increase in Ks, matching the pattern observed in the FORCs from which the weight function was derived. In Fig. 6a, on the other hand, the heavy weights are concentrated in the region for which α is high, such that increases in Ks will only occur when salinity crosses some specific threshold, again reflecting what we expect based on the FORCs. Similar patterns are observed for SAR 50, 25, and 15.
The difference between these two weight functions can be expressed quantitatively by calculating the reversibility index for each. For the weights in Fig. 6a, the reversibility index is Ri=0.33, while for the weights in Fig. 6b the reversibility index is Ri=0.64. In this case, we are comparing two sets of hypothesized FORCs based on the same soil, but we believe the reversibility index's true power is in its ability to compare between different soil types. The index gives us a tool to speak precisely about the degree of reversibility, without having to resort to vague descriptions of soils being “more” or “less” reversible.
It is worth mentioning that the weight “stacks” of Fig. 6 are only half the complete weight function. Their complement is another set of stacks, where the vertical axis is salinity concentration, C, while α and β refer to SAR. For simplicity's sake, these extra stacks are not shown here, but an example can be seen in Fig. C1 of Sect. C.
While the discussion up to now has not considered dynamics, integration of the presented framework with models for salinity and sodicity dynamics would allow us to consider the influence of a soil's previous degradation and rehabilitation cycles on its future state. This is a crucial step in improving our ability to assess the risk of soil degradation because often degradation is triggered by seasonal patterns, which may only lead to significant declines in Ks after a number of years. While models that consider changes in Ks to be reversible are unable to capture the sort of cumulative damage that may result from such cycles, our model is ideally suited for this purpose because it makes future values of Ks dependent on the soil's history. As seen here, we can expect that the soil's weight function will play a crucial role in determining the risk of irreversible soil degradation past critical thresholds. These insights, however, can only be obtained if we widen our attention from traditional degradation experiments and begin to focus on rehabilitation processes through FORC measurements, as demonstrated here.
The model developed here is the first to consider irreversibility in soil degradation as a result of water quality. Specifically focused on changes to Ks as a result of salinity and sodicity, the model demonstrates how hysteresis has the potential to significantly affect our expectations of soil rehabilitation following declines in Ks. The model accounts for soil-specific differences, using a weight function, which allows us to ascertain a given soil's future response to any change in water quality. In this paper, we demonstrated how a soil's weight function can be determined from experimental data, namely the measurement of first-order reversal curves or FORCs. We showcased how to parameterize the model with FORCs obtained from soil column experiments, but in the future, field experiments could increase our understanding of rehabilitation mechanisms in a richer and more realistic environment over longer timescales. Much more than the major loop curves, the FORCs are key to the full characterization of a soil's irreversible behavior. To the best of our knowledge, the scant experimental literature on irreversible decay in Ks (Dane and Klute, 1977) did not exhibit any measurements of FORCs. The experimental results presented here, therefore, are the first to measure FORCs and underscore the need to focus on these measurements, serving as a typical example of the beneficial interaction between theory and experiments.
As a measure of the degree to which a given soil can be rehabilitated following degradation, we introduced a reversibility index. The reversibility index uses the soil's weight functions to quantify the degree to which a given soil is reversible. With this index, one can directly compare the degree to which different soils are susceptible to irreversible degradation. This represents an important advance beyond qualitative statements (e.g., the soil with the greater clay content is more prone to irreversible degradation).
The reversibility index is a convenient quantitative measure; however, it does not determine the risk of irreversible degradation given a particular climate and irrigation regime (Kramer and Mau, 2020). This risk depends on the specific history a soil has been through, and further research is required to unravel the interplay between soil characteristics, environmental drivers, and the actual dynamics of water and solutes in the soil. To this end, a critical next step will be the integration of the model presented here with leading numerical models used for studying the dynamics of soil salinity and sodicity over the long term. At present, all such models consider changes in Ks to be completely reversible (Ri=1), limiting our ability to predict and control the fate of agricultural soils in the face of global climatic change and dwindling water resources.
Finally, while this work focused on the effects of salinity and sodicity on Ks, we believe the conceptual foundation offered by the Preisach framework may have wide-ranging application to other soil processes. As demonstrated here, and in its application to modeling hysteresis in soil moisture patterns (Flynn et al., 2006; McNamara, 2014), a major appeal of the Preisach framework is its generality and lack of major assumptions. This stands in contrast to other models that have been used to describe hysteresis phenomena in soil science, such as those that focus on the response of water retention curves to drying and wetting (Mualem, 1974; Parlange, 1976; Haverkamp et al., 2002). 's 1976 model, for instance, assumes that all hysteresis curves have the same shape (e.g., following van Genuchten's equation). As we argue here, at this point, our understanding on the degree of reversibility in Ks is still based on extremely limited experimental data, and we find no justification for the significant assumption that all reversal curves should have shape similarity with the major hysteresis loop curves. We believe that the Preisach approach is advantageously suited to studying hysteresis in Ks, specifically because of its ability to consider various secondary factors, including soil texture and mineralogy.
Beyond salinity and sodicity, the framework may be adapted to address research questions related to covariates, including composition, clay mineralogy, and soil pH among others. We do note that, while the mathematical formalism admits any number of inputs, that does not mean that the optimal way to use the presented framework would be to increase the number of input variables. We believe that using more than two inputs would obfuscate the understanding of the hysteresis process. Each additional input variable would require a corresponding weight function, such that more than two input variables would require an impractical number of experiments for determining the weights. Nonetheless, we believe that the framework presented here may be of great use to studying the aforementioned soil processes in their own right, especially because many of them are marked by hysteresis. To that end, we encourage readers to make use of our code, all of which is accessible in the repository for this paper.
For the case of a single input variable, we can use the Preisach framework to determine the output value, f(t), for a given input, u(t), by integrating over the limiting triangle such that, in the following:
where is an operator that indicates the state of the hysteron, i.e., whether it is on or off given the value of u(t), and μ(α,β) is the weight function (Mayergoyz, 2003).
The weight function, μ(α,β), can be determined experimentally, as described in Sect. 3, by measuring output along the FORCs. Specifically, the experimental procedure calls for starting at positive saturation (u=umax) and then reducing input so that . The output on the FORCs then corresponds to the values as input is gradually increased from β′ to , as shown in Fig. A1. The output at points (α′, β′) is denoted . The weight at a given point (α′, β′) on the FORC can be found according to the following equation (Mayergoyz, 2003):
For the case when there are two inputs, u(t) and v(t), then Eq. (A1) is expanded so that, in the following:
where μ and ν are the respective weight functions (Mayergoyz, 2003). The dependency of μ and ν on v and u, respectively, serves to couple the inputs. The weight functions μ and ν are given by the following (Mayergoyz, 2003):
where F and G are defined as follows (Mayergoyz, 2003):
Here, denotes the “branching” point of the FORC output for a fixed value of v(t). Likewise, fαβv corresponds to the output as we move along the FORC by increasing the value of α, following the same methodology as in the one input case.
We measured saturated hydraulic conductivity, Ks, in soil columns (length – 5.0 cm; internal diameter – 5.4 cm) within transparent perspex pipes (length – 8.0 cm; external diameter – 6.0 cm), where the hydraulic head H=31 cm was maintained constant using a Mariotte's bottle. The air-dried soil was sieved through a 2 mm sieve and packed to a bulk density of 1.30 g cm−3. For each column, a soil mass of 149 g, divided into five equal parts, was packed into the pipe and bulked with equal pressure to ensure uniform packing across the length of the columns and between all 15 columns. The soil was lined with a geotex fiber, and the pipes were closed at both ends with rubber stoppers. The soil columns were allowed to capillary wet with a solution of concentration 200 mmolc L−1 and SAR equal to 20 (mmolc L−1) for at least 4 h. Upon saturation, leaching started from the top of the soil column. While SAR equal to 20 (mmolc L−1) was kept constant, solutions of decreasing concentrations were made to flow through the columns (200, 80, 50, 20, 15, 10, and 5 mmolc L−1). For each concentration, Ks was measured using Darcy's equation only upon equilibration of the efflux with respect to flow rate, electrical conductivity, and sodium concentration. Electrical conductivity and sodium concentration were measured with a Mettler Toledo InLab 752-6mm and perfectION comb NA, respectively. A total of three replicates were used for each Ks measurement.
The soil used in the experiment is a brown steppe soil (Ravikovitch, 1992) from the Kiryat Gat region of Israel. The measured soil texture is 25 % clay, 25 % silt, and 50 % sand (PARIO soil particle size analyzer; METER Group, Inc.). The dominant clay mineral is montmorillonite (60 %), and the soil's CaCO3 content is 16 %.
It is easiest to start defining the reversibility index Ri for the case of a single input variable as follows:
where I denotes the integral over the limiting triangle as follows:
and, in the following:
is the Euclidean distance of a point (α, β) in the limiting triangle to the diagonal α=β. The reversibility index is , where Ri=1 means total reversibility, and Ri=0 means total irreversibility. In this definition, we assume that the input varies between . In the more general case that the input u varies between , then the denominator in Eq. (C1) should be substituted by umax−umin and the integration limits in Eq. (C2) accordingly.
The intuition behind this definition is as follows. Every point (α, β) in the weight function μ is itself weighted by its distance to the diagonal α=β, given by the function 𝒟(α,β). If all the weights in μ are concentrated along this diagonal, then we have that 𝒟=0, and the reversibility index yields 1, meaning complete reversibility. On the other extreme, if all the weight in μ is as far as possible from the diagonal, namely concentrated at the point (α=umax, β=0), then the reversibility index yields 0, i.e., total irreversibility. The denominator I[μ(α,β)] in Eq. (C1) simply normalizes the weight function to 1, while the factor normalizes the maximal distance to 1. Roughly speaking, the more weight is concentrated near the diagonal α=β, the thinner the major hysteresis loop and the higher the reversibility index, Ri. Conversely, the further the weight is from the diagonal, the wider the major hysteresis loop, and consequently, Ri will be closer to zero. Of course, the factor in Eq. (C1) could be canceled out with the same factor in the distance function 𝒟. We chose to leave 𝒟 as defined to maintain the link between the reversibility index and the distance of the weights from the diagonal α=β.
We can now expand the definition of Ri for two input variables as follows:
where now I denotes the following:
The brackets in Eq. (C4) simply show the arithmetic mean of the same operation done in Eq. (C1), computed for the weight densities μ and ν.
C1 Discrete version of Ri
Often the weight functions μ and ν will be calculated numerically from experimental data, as done in Sect. 4 (see also our full code and widget on GitHub). In this case, μ (ν) will be composed of N (M) discrete stacked layers, each composed by right triangles of side M (N) (see Fig. C1). The reversibility index in the discrete case then becomes the following:
and for the particular case of one input only, it reads as follows:
Note that represents the double summation . Figure C1 depicts an example of the discrete weight functions for N=3 and M=5.
C2 Some simple cases
We calculate here the one input reversibility index of three simple examples in order to build some intuition of how Ri relates to specific weight distributions.
C2.1 Uniform weight distribution
For uniform , the reversibility index reads as follows:
C2.2 Thin diagonal weight distribution
The weight distribution , where corresponds to an infinitesimally thin diagonal line in the (α, β) plane (δ is the Dirac delta). This diagonal line can be offset from the origin by c0umax. Its reversibility index reads as follows:
For c0=0 there is no offset, and the diagonal line sits at α=β; therefore, Ri=1. As the offset is increased, Ri decreases until Ri=0 for c0=1.
C2.3 Dirac delta weight distribution
In the case that all the weight distribution is concentrated at a single point (), the distribution reads , where a0, . The reversibility index then yields the following:
Again, when the weight distribution is on the diagonal line of the limiting triangle (a0=b0), we have Ri=1. As the locus of the distribution moves away from this diagonal, the reversibility index decreases.
All code described in the paper is available in our Zenodo repository at https://doi.org/10.5281/zenodo.4013423 (Kramer et al., 2020). This repository also includes interactive widgets designed to supplement the theoretical overview of the model.
IK and YM developed the mathematical model. IK and YB wrote the code for the model and developed the widgets that accompany this paper. TA and YM designed the experimental setup, and TA conducted the experiment and data analysis. YB and IK did the bibliographic research concerning experimental evidence of irreversible decay in Ks. IK and YM wrote the first version of the paper, and all authors discussed the content, provided feedback, and wrote the final paper.
The authors declare that they have no conflict of interest.
The authors would like to thank Yona Chen, Guy Levy, Naftali Lazarovitch, and Simon Emmanuel for the useful discussion. We would also like to thank the two reviewers, Willem Vervoort and Denis Flynn, for their invaluable feedback and suggestions.
This paper was edited by Gerrit H. de Rooij and reviewed by Willem Vervoort and Denis Flynn.
Ali, A., Biggs, A. J., Šimůnek, J., and Bennett, J. M.: A pH-Based Pedotransfer Function for Scaling Saturated Hydraulic Conductivity Reduction: Improved Estimation of Hydraulic Dynamics in HYDRUS, Vadose Zone J., 18, 1, https://doi.org/10.2136/vzj2019.07.0072, 2019. a
Bennett, J. M., Marchuk, A., Marchuk, S., and Raine, S. R.: Towards predicting the soil-specific threshold electrolyte concentration of soil as a reduction in saturated hydraulic conductivity: The role of clay net negative charge, Geoderma, 337, 122–131, https://doi.org/10.1016/j.geoderma.2018.08.030, 2019. a
Chen, Y. and Banin, A.: Scanning electron microscope (SEM) observations of soil structure changes induced by sodium-calcium exchange in relation to hydraulic conductivity, Soil Science, 120, 428–436, 1975. a, b, c
Cook, F. J., Jayawardane, N. S., Rassam, D. W., Christen, E. W., Hornbuckle, J. W., Stirzaker, R. J., Bristow, K. L., and Biswas, T. K.: The state of measuring, diagnosing, ameliorating and managing solute effects in irrigated systems, CRC for Irrigation Futures Technical Report, CRC for Irrigation Futures, Darling Heights, 2006. a
Dane, J. H. and Klute, A.: Salt Effects on the Hydraulic Properties of a Swelling Soil, Soil Sci. Soc. Am. J., 41, 1043–1049, https://doi.org/10.2136/sssaj1977.03615995004100060005x, 1977. a, b
Dang, A., Bennett, J. M., Marchuk, A., Biggs, A., and Raine, S.: Evaluating dispersive potential to identify the threshold electrolyte concentration in non-dispersive soils, Soil Res., 56, 549–559, https://doi.org/10.1071/SR17304, 2018a. a
Dang, A., Bennett, J. M., Marchuk, A., Biggs, A., and Raine, S.: Quantifying the aggregation-dispersion boundary condition in terms of saturated hydraulic conductivity reduction and the threshold electrolyte concentration, Agr. Water Manage., 203, 172–178, https://doi.org/10.1016/J.AGWAT.2018.03.005, 2018b. a
De Menezes, H. R., De Almeida, B. G., De Almeida, C. D. G. C., Bennett, J. M., Da Silva, E. M., and dos S. Freirer, M. B. G.: Use of threshold electrolyte concentration analysis to determine salinity and sodicity limit of irrigation water, Revista Brasileira de Engenharia Agrícola e Ambiental, 18, s53–s58, 2014. a
Ezlit, Y. D., Bennett, J. M., Raine, S. R., and Smith, R. J.: Modification of the McNeal Clay Swelling Model Improves Prediction of Saturated Hydraulic Conductivity as a Function of Applied Water Quality, Soil Sci. Soc. Am. J., 77, 2149, https://doi.org/10.2136/sssaj2013.03.0097, 2013. a, b, c, d
FAO and ITPS: Status of the World's Soil Resources, Food and Agriculture Organization of the United Nations and Intergovernmental Technical Panel on Soils, Rome, Italy, 2015. a
Flynn, D. and Rasskazov, O.: On the integration of an ODE involving the derivative of a Preisach nonlinearity, J. Phys.: Conf. Ser., 22, 43–55, https://doi.org/10.1088/1742-6596/22/1/003, 2005. a
Flynn, D., McNamara, H., O'Kane, P., and Pokrovskii, A.: Chapter 7 – Application of the Preisach Model to Soil-Moisture Hysteresis, in: The Science of Hysteresis, edited by: Bertotti, G. and Mayergoyz, I. D., Academic Press, Oxford, 689–744, https://doi.org/10.1016/B978-012480874-4/50025-7, 2006. a, b
Goldenberg, L., Magaritz, M., and Mandel, S.: Experimental investigation on irreversible changes of hydraulic conductivity on the seawater-freshwater interface in coastal aquifers, Water Resour. Res., 19, 77–85, 1983. a
Guyer, R. A.: Chapter 6 - Hysteretic Elastic Systems: Geophysical Materials, in: The Science of Hysteresis, edited by: Bertotti, G. and Mayergoyz, I. D., Academic Press, Oxford, 555–688, https://doi.org/10.1016/B978-012480874-4/50024-5, 2006. a
Haverkamp, R., Reggiani, P., Ross, P. J., and Parlange, J.-Y.: Soil water hysteresis prediction model based on theory and geometric scaling, in: Environmental Mechanics: Water, Mass and Energy Transfer in the Biosphere, edited by: Raats, P. A. C., Smiles, D., and Warrick, A. W., American Geophysical Union, Washington, D.C., 213–246, 2002. a
Hillel, D.: Salinity Management for Sustainable Irrigation: Integrating Science, Environment, and Economics, Tech. rep., The World Bank, Washington, D.C., 2000. a
Kramer, I. and Mau, Y.: Soil Degradation Risks Assessed by the SOTE Model for Salinity and Sodicity, Water Resour. Res., 56, 10, https://doi.org/10.1029/2020WR027456, 2020. a, b, c
Kramer, I., Bayer, Y., and Mau, Y.: GitHub Repository – Modeling Irreversible Soil Degradation and Rehabilitation with the Preisach Framework, Zenodo, https://doi.org/10.5281/zenodo.4013423, 2020. a, b
Lagerwerff, J., Nakayama, F., and Frere, M.: Hydraulic Conductivity Related to Porosity and Swelling of Soil, Soil Sci. Soc. Am. J., 33, 3–11, https://doi.org/10.2136/sssaj1969.03615995003300010008x, 1969. a, b
Levy, G. J., Goldstein, D., and Mamedov, A. I.: Saturated Hydraulic Conductivity of Semiarid Soils: Combined Effects of Salinity, Sodicity, and Rate of Wetting, Soil Sci. Soc. Am. J., 69, 653–662, https://doi.org/10.2136/sssaj2004.0232, 2005. a
Mayergoyz, I. D.: Mathematical models of hysteresis, Phys. Rev. Lett., 56, 1518–1521, https://doi.org/10.1103/PhysRevLett.56.1518, 1986. a
Mayergoyz, I. D. (Ed.): Mathematical Models of Hysteresis and Their Applications, Elsevier, New York, New York, https://doi.org/10.1016/B978-0-12-480873-7.X5000-2, 2003. a, b, c, d, e, f, g
McNamara, H.: An estimate of energy dissipation due to soil-moisture hysteresis, Water Resour. Res., 50, 725–735, https://doi.org/10.1002/2012WR012634, 2014. a, b
McNeal, B. L.: Prediction of the effect of mixed-salt solutions on soil hydraulic conductivity, Soil Sci. Soc. Am. J., 32, 190–193, https://doi.org/10.2136/sssaj1968.03615995003200020013x, 1968. a, b
McNeal, B. L. and Coleman, N. T.: Effect of solution composition on soil hydraulic conductivity, Soil Sci. Soc. Am. J., 30, 308–312, https://doi.org/10.2136/sssaj1966.03615995003000030007x, 1966. a, b, c, d, e, f, g, h, i, j, k, l, m, n
Mualem, Y.: A conceptual model of hysteresis, Water Resour. Res., 10, 514–520, https://doi.org/10.1029/WR010i003p00514, 1974. a
O'Kane, J. P. and Flynn, D.: Thresholds, switches and hysteresis in hydrology from the pedon to the catchment scale: a non-linear systems theory, Hydrol. Earth Syst. Sci., 11, 443–459, https://doi.org/10.5194/hess-11-443-2007, 2007. a
Parlange, J.-Y.: Capillary hysteresis and the relationship between drying and wetting curves, Water Resour. Res., 12, 224–228, 1976. a, b, c
Qadir, M., Quillérou, E., Nangia, V., Murtaza, G., Singh, M., Thomas, R. J., Drechsel, P., and Noble, A. D.: Economics of salt-induced land degradation and restoration, Nat. Resour. Forum, 38, 282–295, https://doi.org/10.1111/1477-8947.12054, 2014. a
Quirk, J. and Schofield, R.: The Effect of Electrolyte Concentration on Soil Permeability, J. Soil Sci., 6, 163–178, https://doi.org/10.1111/j.1365-2389.1955.tb00841.x, 1955. a, b, c
Ravikovitch, S.: The Soils of Israel: Formation, Nature and Properties, HaKibbutz Hameuchad Publishing House, Springer, Berlin, Heidelberg, 1992. a, b
Rengasamy, P., Greene, R., Ford, G., and Mehanni, A.: Identification of dispersive behaviour and the management of red-brown earths, Soil Res.h, 22, 413–431, 1984. a
Russo, D.: Numerical analysis of the nonsteady transport of interacting solutes through unsaturated soil: 1. Homogeneous systems, Water Resour. Res., 24, 271–284, https://doi.org/10.1029/WR024i002p00271, 1988. a, b
Shainberg, I. and Singer, M. J.: Soil Response to Saline and Sodic Conditions, in: Agricultural Salinity Assessment and Management, edited by: Wallender, W. W. and Tanji, K. K., American Society of Civil Engineers, Reston, VA, 139–167, https://doi.org/10.1061/9780784411698.ch05, 2011. a
Yaron, B. and Thomas, G. W.: Soil Hydraulic Conductivity as Affected by Sodic Water, Water Resour. Res., 4, 545–552, https://doi.org/10.1029/WR004i003p00545, 1968. a, b, c
Yaron, B., Dror, I., and Berkowitz, B.: Contaminant-induced irreversible changes in properties of the soil–vadose–aquifer zone: an overview, Chemosphere, 71, 1409–1421, 2008. a