Hysteresis in soil hydraulic conductivity as driven by salinity and sodicity: a modeling framework

Declines in soil saturated hydraulic conductivity (Ks) as a result of saline and sodic irrigation water are 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 5 Ks lies at the center of our 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 10 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 novel 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 15 measure “reversal curves,” 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. 1 https://doi.org/10.5194/hess-2020-455 Preprint. Discussion started: 9 November 2020 c © Author(s) 2020. CC BY 4.0 License.

1 Introduction 20 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 hectares of farmland are irretrievably lost due to salinity and sodicity induced degradation, while crop output is decreased on an additional 20 to 46 million hectares (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). Increased strain on freshwater resources, growing human populations, and 25 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, K s , 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 K s through a combination of different 30 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 K s (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 K s are characterized by smooth curves as water quality is changed, rather than by 35 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. 40 While numerous experimental studies have led to models describing declines in K s 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 K s can be restored simply by changing the chemical properties of the soil water (McNeal and in K s , 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., 50 1984; Cook et al., 2006;Bennett et al., 2019). Likewise, the experimental evidence that has considered degradation and rehabilitation together suggests that changes in K s 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 our 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 55 considered, it has been shown that seasonally changes in K s may be compounded, such that initially small declines in K s 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 to 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 60 understanding how hysteresis affects the processes of degradation and rehabilitation.
In this paper, we present a model for the effect of salinity and sodicity on K s that considers hysteresis. That is, in contrast to existing models, ours accounts for how a soil's history of degradation and rehabilitation will affect K s 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, our model has the potential to significantly enhance 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 et al., 2006;McNamara, 2014) and the effect of stress and strains on soil and rock structures (Guyer, 2006). A comprehensive 70 overview of the Preisach model can be found in Mayergoyz (1986Mayergoyz ( , 2003. We begin this paper by introducing the Preisach framework and how it can be applied to model hysteresis in soil K s . 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 parametrized 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 75 future assessments of the risk of soil degradation. Our modeling of hysteresis in K s is based on Preisach's framework. We will briefly present this framework, and then apply it to our specific needs for modeling K s . 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 80 of the figures in the following sections.

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 85 values, α and β, such that: 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 90 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 Figure 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) ≥ α.

95
Decreases in input are modeled by turning "off" hysterons for which u(t) ≤ β (see Widget 1 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 Figure 2, where the panels represent the following steps: The hysteron forms the foundation of the Preisach model, making output (f ) dependent on threshold input values (α and β) and the previous output value. (b) The study of hysteresis on the system level is enabled by considering an infinitude of hysterons spread over a limiting triangle, in which each coordinate pair corresponds to specific threshold values (αi, βi).
(a) We begin at t = 0 with input u = u 0 . 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.
(b) Input u is decreased to u 1 . We swipe left on the limiting triangle, turning off all hysterons for which β > u 1 , which causes a decrease in f .
(c) Input is increased to u 2 . We swipe up on the limiting triangle, turning on all hysterons for which α < u 2 , 105 which leads to an increase in f .
(d) A further decrease in input to u 3 will again cause a swipe left, turning off all hysterons for which β > u 3 , and leading to a decline in f .

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 K s as the output, 110 while the inputs are given by the chemical composition of the soil water. In our model, the two inputs are electrolyte Hysteron Off Hysteron On concentration (C, mmol c L −1 ) and Sodium Adsorption Ratio (SAR, mmol 1/2 c L −1/2 ). In doing so, the framework presented in subsection 2.1 must be expanded to accommodate two input variables, but this is easily accomplished (see Appendix A). In addition to our two input variables (C and SAR), many covariates can have an important role in the K s response, such as clay content, clay mineralogy, soil pH, organic matter, etc. The Preisach framework is 115 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 (and highly elucidating) example of the role of the weighting function can be found in Widget 2 of our GitHub repository.

120
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, Figure 3 presents three weight functions together with their corresponding hysteresis curves. Every weight function (lefthand column) is associated with a unique set of hysteresis curves showing how K s responds to changes in C (right-hand column). The hysteresis graphs depict major loop curves (solid lines), as well as first-order reversal 125 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 K s is fully degraded. That is, when the gap is thin, as in Figure 3a, degradation in K s is highly reversible. When the gap between the two curves of the major hysteresis 130 loop is wider, as in Figure 3c, degradation in K s 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 K s is limited. The closer the FORCs are to the top curve of the major loop, the higher the degree of reversibility. The FORCs in Figure 3 show high, medium, and low degrees of reversibility, for panels (a)-(c) respectively. Because the weight function uniquely defines a hysteresis graph, characterizing the degree to which a given soil is susceptible to irreversible 135 changes in K s can be done by finding that soil's weight function µ(α, β).

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, R i . The full development of the reversibility index is contained in Appendix C, but the fundamental idea behind it is that the closer the weight distribution is to the diagonal of the limiting triangle 140 (i.e., the line α = β), the more reversible the soil, as in Figure 3a. Conversely, as the weight distribution becomes closer to the vertex of the triangle on the top left (α max , 0), the more irreversible the soil, as in Figure 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.
3 Determination of the weight function 145 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 procedure required 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 150 show limited FORC measurements and use these to construct a partial weight function.

Experimental setup
Experiments to measure FORCs should begin with the same setup commonly used to study degradation of K s 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 K s 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 K s are recorded as the quality of the input water is changed. In the context of the Preisach framework, these experiments begin from "positive saturation," the state when all hysterons are 160 "on." This corresponds to the point where both K s and the input variable are maximal, as seen in the annotation in Figure 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 K s , we "reverse" the trend in water quality (now either towards higher C or lower SAR) and 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 Figure 3, two FORCs (the dashed lines) are shown for each hysteresis graph. Because our model considers changes in K s as a result of both C and SAR, calculation of a particular soil's complete weight function requires the measurement of two sets of 170 FORCs, one in which changes in K s 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., K s ) 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 in the GitHub repository. We note that in the case where SAR is varied while C is kept constant, positive saturation occurs for low input values, 175 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: K s declines as the inputs are decreased.

Initial results
To showcase how FORCs can be used to calculate a soil's weight function, we conducted a limited experiment to 180 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, our experiments do not constitute the whole set of procedures delineated above, 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 K s , as influenced by salinity and sodicity.

185
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 = 20 (mmol c L −1 ) −1/2 constant, and measured K s as the value of C was repeatedly decreased and then increased. Measurements of K s 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 was three hours, after which determination of K s 190 was conducted, and the solution was changed for the next step of the procedure. For more details of the soil properties and experimental setup, see Appendix B. In what follows, we will refer to each FORC by its lowest concentration value.
The degradation in K s experienced by the soil before salinity was decreased below C = 20 mmol c L −1 was found 195 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 K s not returning to its original value, even for very concentrated solutions of 400 mmol c L −1 (not shown in the graph).
We used the data from the FORCs in Figure 4a to determine a limited weight function, shown in Figure 4b.
This is a limited weight function because it contains information on one SAR value only. Because the decay in The expansion of the triangle leaves some weight beyond the experimental bounds (see dotted line in Figure 4b), thus allowing for incomplete rehabilitation.
In order to make sure that the weight triangle contains indeed the information of the hysteresis graph, we modeled the soil's response according to this triangle. The modeled results, shown in Figure 4c, are quite similar 210 to the experimental data of Figure 4a. It is worth noting that the modeled FORCs return to higher K s values, i.e., there is a slight overestimation of rehabilitation, and several FORCs converge together as C increases. This is due to our 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 Figure 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 215 upscaling involved in the process. In our estimation, the algorithm devised to reconstitute the weight function from experimental data (see the GitHub repository) is able to impart both the qualitative and quantitative behavior of the soil under salinity and sodicity changes.

Insights for future experiments
This limited experiment yields a number of important insights, and we hope that these will prove useful to future 220 efforts to characterize trends of partial reversibility in soil K s .  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 Figure 4a) and for five fixed C values (i.e., when SAR is varied). Each of these hysteresis 225 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 impair their usefulness.

230
Given the typical times of equilibration we experienced, we estimate that a full characterization of the soil hysteresis would take approximately one 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 235 decay), to determine the range of input values where K s varies the most. The points of return of each FORC should be distributed roughly equally along the decay in K s . For instance, Figure 4a shows that the return points are more concentrated between concentrations 5 and 20 mmol c L −1 , because there return points roughly sample equally the steepest decline in K s . Conversely, return points where the slope of the major loop is mild, do not yield as much useful information for reconstituting the weight function. It is important to note that the measurement points 240 along all FORCs have the same input values. This is necessary because of the way we developed our algorithm for building the weight functions.

Model applications
In this section, we use the model to demonstrate how consideration of hysteresis is essential when analyzing degradation and rehabilitation in soils. While no experimental work has addressed this issue in detail, we 245 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 section 3. Specifically, McNeal and Coleman (1966) measured K s as input water of fixed SAR, but declining salinity, was applied to soil columns. The procedure was repeated for 250 several SAR values, as shown in Figure 5a, which reproduces the degradation curves presented by McNeal and Coleman (1966) for SAR 15, 25, 50, and 100 mmol 1/2 c L −1/2 (henceforth we will omit the units of SAR).
The possibility for different levels of reversibility, given the same degradation curves, is presented in Figure 5b-c.
Here, we use the degradation curve associated with SAR 50 to construct two potential sets of FORCs. The FORCs in Figure 5b correspond to a low level of reversibility -once the soil has been degraded, an increase to higher 255 K s values requires a significant increase in salinity. By contrast, the FORCs in Figure 5c indicate a high level of reversibility, such that K s 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 section 3, which will allow us to more accurately determine the true extent to which a given soil exhibits hysteresis 260 in degradation and rehabilitation. We also note that even the concave FORCs in Figure 5c consider a higher level of irreversibility than existing models for degradation and rehabilitation, which allow K s to increase by returning along the original degradation curve (i.e., the top curve of the major loop).
By extending the approach used in Figure 5b-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 265 high reversibility) for each of the respective SAR values (15, 25, 50, 100). We then use these FORCs to compute weight functions, using the technique and code described in section 3. The resulting weight functions are presented in Figure 6. The different areas occupied by the limiting triangle in these plots reflect the different values of C max -the unique salinity value at which degradation begins -for each respective SAR value. When determining the weight function, salinity values greater than u max are not assigned any weights.

270
Examination of the weight functions in Figure 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 (Figure 6a), the heavy weights are concentrated exclusively in the region close to (α max , 0). For the weight function derived from the more reversible FORCs (Figure 6b), the heavy weights are similarly concentrated in the region for which β is close to zero. A major contrast can be seen, however, 275 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 280 seen in Figure 5. The α-axis, by contrast, is associated with rehabilitation. Because the weights in Figure 6b are distributed along the entire α-axis, we can expect any increase in salinity to lead to an increase in K s , matching the pattern observed in the FORCs from which the weight function was derived. In Figure 6a, on the other hand, the heavy weights are concentrated in the region for which α is high, such that increases in K s will only occur when salinity crosses some specific threshold, again reflecting what we expect based on the FORCs. Similar patterns are 285 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 Figure 6a the reversibility index is R i = 0.33, while for the weights in Figure 6b the reversibility index is R i = 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 290 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 "more" or "less" reversible soils.
It is worth mentioning that the weight "stacks" of Figure 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 Appendix C, 295 Figure C1.
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 300 declines in K s after a number of years. While models that consider changes in K s 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 K s 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

Conclusions
The model developed here is the first to consider irreversibility in soil degradation as a result of water quality.
Specifically focused on changes to K s 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 K s . The model 310 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. 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 K s (Dane and Klute, 1977) did not 315 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 320 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, such as "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 325 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 K s to be completely reversible (R i = 1), limiting our ability to predict and control the fate of 330 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 K s , 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 335 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). Parlange's (1976) model, for instance, assumes that all hysteresis curves have the same shape (e.g., follow van Genuchten's equation). As we argue here, at this point, our understanding on the degree of reversibility in K s is still based on extremely limited experimental data, and we find no justification for the 340 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 K s , 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 345 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 that 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 350 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.
Code availability. All code described in the manuscript is available in our GitHub repository: https://github.com/yairmau/ hysteresis-python (Kramer et al., 2020). This repository also includes interactive widgets designed to supplement the theoretical 355 overview of the model.

Appendix A: The Preisach framework
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,

360
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 section 3, by measuring output along the FORCs. Specifically, the experimental procedure calls for starting at positive saturation (u = u max ) and then reducing input so that u = β . The output on the FORCs then corresponds to the values as input is gradually 365 increased from β to u = α , as shown in Figure A1. The output at points (α , β ) is denoted f α β . The "weight" at a given point (α , β ) on the FORC can be found according to the equation (Mayergoyz, 2003) µ(α , β ) = 1 2 For the case when there are two inputs, u(t) and v(t), then Equation A1 is expanded so that

370
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 (Mayergoyz, 2003) µ where F and G are defined as follows (Mayergoyz, 2003) 375 Here, f + βv 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.

Appendix B: Experimental setup
We measured saturated hydraulic conductivity, K s , 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 5 equal parts, was 385 packed into the pipe and bulked with equal pressure to ensure uniform packing across the columns length and between all the 15 columns. The soil was lined with a geotex fibre, and the pipes were closed at both ends with The soil used in the experiment is a brown steppe soil (Ravikovitch, 1992) from the Kiryat Gat region of Israel.

Appendix C: Reversibility index
It is easiest to start defining the reversibility index R i for the case of a single input variable: where I denotes the integral over the limiting triangle and is the Euclidean distance of a point (α, β) in the limiting triangle to the diagonal α = β. The reversibility index 405 0 ≤ R i ≤ 1, where R i = 1 means total reversibility, while R i = 0 means total irreversibility. In this definition we assume that the input varies between 0 ≤ u ≤ u max . In the more general case that the input u varies between u min ≤ u ≤ u max , then the denominator in Eq. (C1) should be substituted by u max − u min , 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 410 by its distance to the diagonal α = β, given by the function D(α, β). If all the weights in µ are concentrated along this diagonal, then we have that D = 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 (α = u max , β = 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 √ 2 umax 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, R i . Conversely, the further the weight is from the diagonal, the wider the major hysteresis loop, and consequently R i will be closer to zero. Of course, the factor √ 2 in Eq. (C1) could be cancelled out with the same factor in the distance function D. We chose to leave D as defined to maintain the link between the reversibility index and the distance of the weights from the diagonal α = β. 420 We can now expand the definition of R i for two input variables:

C2 Some simple cases
We calculate here the 1-input reversibility index of three simple examples, in order to build some intuition of how R i relates to specific weight distributions.

445
Its reversibility index reads: In the case that all the weight distribution is concentrated at a single point (α = a 0 u max , β = b 0 u max ), the distribution reads: µ(α, β) = δ(α − a 0 u max )δ(β − b 0 u max ), where a 0 , b 0 ∈ [0, 1]. The reversibility index then Again, when the weight distribution is on the diagonal line of the limiting triangle (a 0 = b 0 ), we have R i = 1. As 455 the locus of the distribution moves away from this diagonal, the reversibility index decreases.
Author contributions. 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 manuscript, and all authors discussed the content, provided feedback and wrote the final paper.