Soil sodicity originating from marginal groundwater

Soil salinity and sodicity are among the oldest soil and groundwater pollution problems and are widespread across the globe. Where salinity affects crop water uptake and yield, sodicity may additionally cause poorly reversible soil structure degradation and a severely reduced hydraulic conductivity. We use the model HYDRUS‐1D to simulate sodicity development in soils with shallow, Na‐rich groundwater under a normal weather regime with distinct dry seasons. Attention is given to the impact of a sudden fresh water input on the formation of a sodic layer. The complex interplay between soil chemistry, soil physics, soil mechanics (as far as swell–shrink behavior is concerned), and fluctuating atmospheric conditions results in a remarkably regular relation between depth, location, and severity of a sodic layer that forms within the soil as a function of rainfall intensity. A threshold behavior is observed: sodic layer formation is absent at rainfall intensities below this threshold, whereas sodic layer thickness and hydraulic conductivity reduction increase rapidly with intensities exceeding this threshold. This is the case even for different soil types and groundwater depths. Field observations agree with our simulations: the properties of the layer with sodicity‐induced structure degradation are more strongly developed, as this layer is situated at a shallower depth. The implementation of hydraulic conductivity reduction as a function of exchangeable Na percentage and ionic strength in HYDRUS‐1D can be improved towards a smooth reduction function, changing soil physical parameters due to swelling and dispersion of clay and reconsideration of the reversibility of sodicity development.


INTRODUCTION
Water deficiency for crops and natural vegetation is widespread, not just in arid and semiarid regions, but also in humid regions with a distinct dry season. The latter has 2 of 14 VAN DE CRAATS ET AL. Vadose Zone Journal example. If leaching is not pursued or is insufficient, salinization of the root zone will occur that at some time makes primary (e.g., crop) production impossible (Qadir, Noble, Schubert, Thomas, & Arslan, 2006).
Whereas soil salinity may be prevented, or remediated, rather easily by appropriate leaching of salts, a more stealthy process of soil degradation is the development of soil sodicity. Sodicity generally occurs hand in hand with salinity problems, especially in areas with Na-rich ground or surface water. Tanji and Wallender (2011) and Wicke et al. (2011) estimate that around 1 billion ha of land is salt affected, of which 40-60% is threatened by sodicity problems as well. A major part of salt-affected lands is found in semiarid regions, where the options for mitigation, such as soil amendments or a change in irrigation water quality or quantity, can be limited due to environmental or socioeconomic factors (Qadir et al., 2006).
The main physicochemical aspects of sodicity have been documented already by Bolt (1979) and Bresler, MacNeal, and Carter (1982); Richards (1954). Sodicity, its name derived from "sodium," is characterized by a relative accumulation of monovalent Na + at the expense of, for example, divalent cations of Ca 2+ and Mg 2+ at adsorption sites of soil particles. These divalent ions normally comprise the vast majority of the adsorbed ions (Richards, 1954;van der Zee, Shah, & Vervoort, 2014). Sodicity not only implies an accumulation of Na in the soil, which in itself can be toxic to plants or lead to nutrient imbalances (Qadir & Schubert, 2002), but it also influences the electrochemical behavior of the soil. Adsorbed cations screen off the negative charge of solid soil material (e.g., clay colloids). At a large ionic strength of the soil solution (i.e., a high salinity), screening is very effective and independent of the valence of cations. At a small ionic strength, however, this screening is less effective, and the charge plays a role; accordingly, multivalent cations (Ca 2+ , Mg 2+ ) screen off better than monovalent ones (Na + , K + ) and are preferentially adsorbed in the electrical double layer of charged soil minerals.
This process plays an important role in regions with distinct dry and wet seasons (or irrigation seasons) in semiarid regions, but also in more humid regions, where capillary rise of, or irrigation with, Na-rich water occurs. For these situations the exchangeable sodium percentage (ESP) of the adsorption sites of solid soil materials will gradually increase, even if year-to-year salinity remains constant (Gonçalves et al., 2006;van der Zee, Shah, van Uffelen, Raats, & dal Ferro, 2010). During dry periods, a large ionic strength within the soil water develops due to evapotranspiration such that the sorption complex does not favor cations of higher valence. Thus, Na (present in high concentrations in soil water) replaces Ca (present in lower concentrations) on the exchange sites. In view of the high pH in sodic soils, Ca might precipitate as gypsum, removing it from the soil solution, further enhancing relative Na concentrations in the Core Ideas • Sodic layer formation by fresh water infiltration is modeled with HYDRUS-1D. • Infiltration intensity affects sodic layer thickness, depth, and severity. • Upon exceeding a threshold rain intensity, soil structure degradation rapidly intensifies. • Groundwater depth and rainfall events are important for soil structure deterioration. • Feedbacks in sodicity-induced structure degradation require better model implementation.
soil solution (Sparks, 2003). During a subsequent wet or irrigation season with sufficient leaching water of a good quality, the ionic strength of the soil solution decreases, and these conditions favor divalent ion sorption. However, as both Ca and Na cations are leached, replacement of Na by Ca is limited, leading to a net increase in ESP over time.
The screening efficiency as a function of soil solution concentration also influences the physical and mechanical behavior of soils (Qadir & Schubert, 2002;Quirk & Schofield, 1955). Depending on the composition of the solid-phase mineralogy, soils may swell and shrink. The tendency to swell increases if the ionic strength of the soil solution decreases, for instance due to infiltration of fresh precipitation water (Cornelissen, Leijnse, Joekar-Niasar, & van der Zee, 2019;Minhas & Sharma, 1986). This tendency is much stronger when a larger portion of adsorbed ions consists of monovalent, instead of divalent ions. Upon swelling, larger pores are compressed, whereas a much larger abundance of smaller pores grow slightly in size. In addition, dispersion of clay colloids due to a weak cohesion between the soil particles can result in the formation of a more dense layer, clogging pores and leading to surface crusting (on top of the soil) or hard setting (within the soil) (Qadir & Schubert, 2002).
The outcome of these processes is an overall distinct reduction in hydraulic conductivity (Quirk & Schofield, 1955;van der Zee et al., 2014). A layer with dispersed clay can obtain a very low hydraulic conductivity and will be mechanically very hard, impeding water flow as well as biological activity. Such layers can only be disrupted by mechanical means (Shainberg & Letey, 1984), whereas sodic layers, which still provide some degree of hydraulic conductivity, can also be remediated by chemical amendments, phytoremediation, or tillage (Minhas, Qadir, & Yadav, 2019).
Because of the poor reversibility of sodicity-induced structure degradation (SISD), prevention is the motto. It requires that upcoming sodicity development is recognized at a sufficiently early stage (although it is not very visible), or that the logic that it will develop can be communicated convincingly and in time. A slow buildup of ESP can be halted and reversed by taking the appropriate measures if identified in time, before any serious decrease in hydraulic conductivity or dispersion of clay has occurred (Qadir & Schubert, 2002).
One way to recognize the risk of SISD is the use of models to predict changes in the water and salt balance, including the effects of sodicity. Several model approaches were developed (Mau & Porporato, 2015;Šimůnek & Suarez, 1997;van der Zee et al., 2014), which differ in complexity in, for instance, how the processes described above are incorporated, as well as in spatial and temporal time scales involved. In particular, the numerical model HYDRUS-1D (Šimůnek, Šejna, Saito, Sakai, & van Genuchten, 2008) has frequently been used as a benchmark to simulate the buildup and remediation of salinity (e.g., Sun et al., 2019), and in combination with the major ion chemistry module, sodicity (Gonçalves et al., 2006;Ramos et al., 2011; in soil profiles. As follows from the above, a sudden fresh water input can lower the ionic strength of the soil solution and cause swelling of clay, dispersion of soil particles, and changes in hydraulic conductivity (Minhas & Sharma, 1986). The scope of this paper is therefore to simulate the development of salinity and sodicity under field conditions by capillary rise of marginal (too brackish) Na-rich groundwater and in response to sudden fresh water infiltration. The results help us to anticipate pending soil structure degradation and threshold behavior of this process. Our experiences with modeling the processes with HYDRUS-1D are used to suggest improvements regarding the software and process understanding.

Conceptual model
Our conceptual model is based on the presence of a sufficient quantity of swelling and shrinking clayey particles in the soil to make SISD likely. This implies that sandy soils as such are of less interest, as soil structure degradation is not probable. We consider the situation that the salts originate from the groundwater; we do not consider the use of marginal irrigation water in this paper. Rather, infiltrating water concerns good quality rainfall. Such conditions occur, for instance, in Australia (George, McFarlane, & Nulsen, 1997) and Hungary (Van Asten, 1996).
The soil is assumed to be vertically homogeneous. Different from the parsimonious approach by van der Zee et al.
(2014), we take a distributed model approach, where the onedimensional soil column is a buildup of a large number of soil cells in series, ranging between the soil surface, where rainfall occurs, and the groundwater level, where capillary rise starts. Our motivation to step away from simpler approaches, which consider a perfectly mixed plow layer yet represent the subsoil between this layer and the groundwater level by simple functions (Rodriguez-Iturbe & Porporato, 2004;Vervoort & van der Zee, 2008, is that field observations demonstrate the typical sequence of different soil layers or "horizons" for sodic soils. These are the A horizon (or E horizon, from which plants extract their transpiration water and from which clay particles can be dispersed and transported downward), the textural B horizon (with enriched, compressed, and densified clay and enhanced Na concentrations; also known as sodic or natric horizon; IUSS Working Group WRB, 2015), and the C horizon (subsoil, or parent material). Examples of such sequences are found in, for instance, Hungary (Toth & Jozefaciuk, 2002;Van Asten, 1996). Of course, our analysis assumes initially a homogeneous soil but aims to assess how the three-horizon composition of swelling and shrinking soils may be developing.
Typical for the sodic soil is the textural B horizon. Important properties of this B horizon are the depth in the soil profile where this layer starts, its intensity (or how high ESP is), and its thickness. For sodic soils that develop not by using marginal irrigation water, but by capillary rise of marginal groundwater, it is plausible that the depth of groundwater is a dominant genetic feature: groundwater that is too deep may not provide an adequate amount of capillary rise to facilitate transport of sufficient salts for the development of a sodic layer (Schofield, Thomas, & Kirkby, 2001;Sun et al., 2019).
In our analysis, we have to admit that coupled soil physical (i.e., soil hydrology) with soil mechanical models appear to be lacking. This also appears to be the case for software that accommodates the gradual change in soil physical parameters as a function of changing state parameters (volumetric water fraction or pressure, salinity, and water composition). The only numerical model that is, to the authors' knowledge, capable of including the impacts of soil chemistry on saturated hydraulic conductivity is HYDRUS-1D (Reading, Baumgartl, Bristow, & Lockington, 2012).

Theory
In our simulations, we use the software HYDRUS-1D (Šimůnek et al., 2008) including the major ion chemistry module UNSATCHEM (Šimůnek & Suarez, 1994) (see Table 1 for descriptions of variables). Water flow in the vadose zone is described with the equation Empirical parameter for van Genuchten functions α 1 -Reduction function for T p due to water stress α 2 -Reduction function for T p due to osmotic stress is a sink/source term for root water uptake [L 3 L −3 T −1 ], and K is the hydraulic conductivity function [L T −1 ] for which the van Genuchten expression (van Genuchten, 1980) is used: where θ s is the saturated and θ r is the residual volumetric water fraction [L 3 L −3 ], respectively, K s is the saturated [L T −1 ] and K r the relative [-] hydraulic conductivity, respectively, 66S e is the relative saturation [-], α [L −1 ], m and n [-] are empirical parameters, and r 1 and r 2 are reduction functions, with values ranging between 0 and 1, that represent the effects of the ESP, total solute concentration (C t ) [mmol c L −1 ], and pH on hydraulic conductivity. Although we are not interested in (pH dependent) chemical precipitation and dissolution reactions, the HYDRUS-1D model allows the use of either both or none of the reduction functions. The pH reduction function r 2 (pH) is given by 1 (pH < 6.83), 3.46 − 0.36 × pH (6.83 ≤ pH ≤ 9.3) or 0.1 (pH > 9.3). The reduction in hydraulic conductivity due to sodicity is the combined effect of ESP and C t for which we use an empirical expression, following McNeal (1968): where c and p are empirical parameters and x is a swelling factor that depends on the weight fraction of the clay mineral montmorillonite (f Mont , set to 0.1), d * is the adjusted interlayer spacing [L] and ESP* is the adjusted ESP, given by where (CNa + ) a is the concentration of Na + at the cation exchange complex, and CEC the cation exchange capacity (both in mmol c kg −1 ).
The adjusted interlayer spacing is d* = 356.4C t −0.5 + 1.2 if C t < 300 mmol c L −1 , and otherwise zero. The value of p is 1 (ESP < 25), 2 (25 < ESP < 50), or 3 (ESP > 50), and for these ESP ranges; the parameter c takes the values of 35, 932, and 25,000, respectively. The reduction in hydraulic conductivity as function of ionic strength and ESP has been approximated by many expressions (see, for instance, the summary given in Ezlit, Bennett, Raine, and Smith, 2013). We adopted this parameterization of McNeal (1968), instead of the smoothed version of van der Zee et al. (2014) or the modified version by Ezlit et al. (2013), because it is encoded this way in HYDRUS-1D. It was already recognized that the experimental foundation of this parameterization is rather shallow and dependent on soil specific properties as mineralogy, but also on the saturation of the soil (Ezlit et al., 2013;Šimůnek & Suarez, 1997;Šimůnek et al., 2008). This parameterization has the additional disadvantage that the reduction function becomes (i) non-smooth and (ii) has a large reduction in hydraulic conductivity even for negligible ESP if the salinity (total concentration) is low, as is apparent from Figure 1 and also identified by Ezlit et al. (2013). Both disadvantages are unrealistic and affected our simulations. Transport of solutes is described by the convectiondispersion equation, given by The tortuosity encoded in HYDRUS-1D is τ = θ 7/3 /θ s 2 , following Millington and Quirk (1961). We assume the local equilibrium assumption (LEA) to be valid and describe the exchange of major cations with the empirical Gapon-type of equation (Bolt, 1979): where v i and v j represent the valences of cation i and j, respectively, and K ex i,j is the Gapon exchange coefficient. Actually, the concentrations in Equation 6 are the activities (i.e., concentrations corrected for the ionic strength of the soil solution, as explained by Bolt, 1979). In the present case, we only consider exchange between Na and Ca, as this captures the main impacts of the salt concentration and its composition on exchange and transport relevant for sodicity, although differences exist between different types of cations (Alperovitch, Shainberg, & Keren, 1981).
In our simulations, we examine root water uptake with a macroscopic approach that considers reduced uptake in case of water or salinity stress (Feddes, Kowalik, & Zaradny, 1978;Maas & Hoffman, 1977). Both stressors are assumed to reduce uptake in a multiplicative sense (van Genuchten, 1987): where the actual sink term S that comprises the actual root water uptake for transpiration is reduced by a correction for drought (α 1 ), osmotic stress (α 2 ), and a normalization of potential transpiration T p [L T −1 ] based on the root density distribution ρ d (Raats, 2007). The correction for drought involves four h i values for the piecewise linear reduction function where h 1 represents the point where anaerobic conditions occur, h 4 is the permanent wilting point (both leading to a transpiration rate of zero, as α 1 = 0), and the range from h 2 to h 3 represents the conditions with optimal root water uptake (i.e. water is transpired at the potential transpiration rate, α 1 = 1). Between h 1 and h 2 , and h 3 and h 4 , the actual transpiration rate is reduced salinities. r 1 is multiplied by the saturated hydraulic conductivity Equation 2 to obtain the reduction in hydraulic conductivity due to sodicity effects. C t is total solute concentration linearly depending on the actual water pressure. The h 3 itself is dependent on the potential transpiration rate, varying linearly between h 3 H and h 3 L (with superscripts H and L denoting high and low) for potential transpiration rates between T p H and T p L , respectively. For potential transpiration rates beyond these limits, h 3 takes the corresponding high or low h 3 value.
In a very similar approach, osmotic stress is accounted for with a critical value of the osmotic pressure head (h πc ) and slope of the reduction factor s π , such that α 2 (h π , z, t) = 1 + (h π − h πc )s π for |h π | > |h πc |, and is equal to one elsewhere. The EC e (electrical conductivity [EC] in saturation extract) values, as frequently prescribed for the critical value and slope, are transformed to osmotic head using h π = −3.8106 × 2EC e + 0.5072 [m], if EC e is in dS m −1 .
Besides ESP Equation 4 as an indicator of soil sodicity and the hazard of structure degradation, we also use the so-called electrochemical stability index (ESI), given by where the saturation percentage (SP) is by weight and given as [θ s (100/ρ)]/(θ s × 100), and dry bulk density is in g cm −3 .

T A B L E 2 Hydraulic function parameters for the loam and clay
loam soils corresponding to Figure 1 and used in the simulations

Simulation data
We simulate two soil types using HYDRUS-1D, a loam and clay loam, which differ in soil hydraulic properties such as saturated hydraulic conductivity and retention function parameters. These parameters are provided in Table 2 for both soil types. Their corresponding pF curves are displayed in Figure 2. All simulations use a vertical cell size (dz) of 0.5 cm. Simulations of the loam soil type account for reduction in hydraulic conductivity due to chemistry (r 1 and r 2 in Equation 2a) during the simulation, whereas we neglect this reduction for the clay soil and only obtain the reduction which would have occurred, by applying Equations 3 and 4. The reason for this is that pH was affecting hydraulic conductivity in the clay loam soil simulation (through r 2 ), whereas this was not the case for the loam soil, impeding water flow and limiting the effect of sodicity on hydraulic conductivity. Initial soil solution and exchangeable cation compositions, as well as solution compositions at the domain boundaries, are given in Table 3. These boundary solution compositions are constant in time. Groundwater is given a composition similar to that of seawater (i.e., a NaCl-type groundwater, with a [lower than seawater] total salt concentration of 1.8 g L −1 ). Initially, no precipitated compounds are present. Pressure heads are initialized by a hydrostatic equilibrium, with groundwater (h = 0) at the profile bottom. This groundwater level is fixed for the entire simulation period. Both soil types have a dry bulk density of 1.54 g cm −3 and a CEC of 40 mmol c kg −1 . The molecular diffusion coefficient is set at 1.28736 cm 2 d −1 , and dispersivity is set at 1 cm. The Gapon exchange coefficient for exchange of Ca and Na is fixed at 2.2794.
The root zone extends to 90 cm below surface for both soil types, with a root distribution function decreasing linearly from one at the surface to zero at the root zone bottom. The crop considered is wheat (Triticum aestivum L.), which is grown between April and November. For the remaining period, the soil is considered bare. Daily precipitation and potential transpiration rates are obtained from the Paynes Find station, Western Australia (BOM, 2018). This station is located in a semiarid region and experiences distinct dry and wet seasons, with an average, but highly variable, precipitation amount of 210 mm yr −1 . The quality and duration of the rainfall and temperature measurements at this location in combination with the specific climate and occurrence of sodic soils in Western Australia makes this data series suitable for our analysis. However, data from other locations matching these criteria could have been used as well. Minimum and maximum temperatures are converted to reference crop evapotranspiration (ET 0 ) following Hargreaves method (Hargreaves & Samani, 1985). We use crop coefficients k c and k cb to calculate potential transpiration (k cb ET 0 ) and evaporation [(k ck cb )ET 0 ] following Allen, Pereira, Raes, and Smith (1998) during the growing season. Coefficients used are k c = 0.7 and k cb = 0.15 for 1-30 April, k c = 1.15 and k cb = 1.1 for 1 May to 31 October, and k c = 0.25 and k cb = 0.15 for 1-30 November. For the remaining period, no transpiration occurs, and potential evaporation is given by 0.5ET 0 .
The main part of the simulation period of 20 yr (using data from 1982 to 2001) is used to build up a salinity profile in the unsaturated zone for a situation with a nonirrigated, groundwater-fed wheat field in a semiarid climate. A rainfall event in the 17th year (26-29 May 1999, 128 mm of precipitation in total) is used to observe effects of a significant fresh water input on the development of a sodic layer in a saltaffected soil. By changing this particular downpour event's intensity, we observe the effect of rainfall intensity (or more general, fresh water input intensity) upon sodic layer development. We do this for a loam soil type with a groundwater depth of 150 (shallow) and 200 cm (deep), as well as for the clay loam with a depth of 200 cm. Furthermore, we examine the effect of groundwater depth on the development of the sodic layer, using reference rainfall intensity and a sequence of depths ranging between 150 and 200 cm.

Field observations
An example of SISD can be found in the Great Hungarian Plain, or Hortobágy area. In this region, shallow water tables are present that provide Na to the root zone by capillary rise (Schofield et al., 2001). Although not classified as semiarid, this subhumid region is characterized by hot, dry summers, where potential evapotranspiration (700 mm yr −1 ) exceeds the total precipitation (550 mm yr −1 ), and intense summer downpours can occur (Schofield et al., 2001;Tóth, Novák, & Rakonczai, 2015). In addition, surface runoff frequently occurs through ephemeral streams due to snowmelt in combination with frozen soil, as well as due to the presence of poorly permeable sodic layers (Novák & Tóth, 2016), effectively increasing the net precipitation deficit. This combination of factors resulted in topsoils enriched in Na where sodicity, as well as structure degradation, occurred  (Toth & Jozefaciuk, 2002). Accordingly, degraded soils (socalled Solonetz) developed with a very poor soil structure. In essence, a hard, compacted, textural B horizon with a low hydraulic conductivity developed (Tóth et al., 2015), which effectively forms a barrier between the fine-textured, thin A (or E) horizon, from which the vegetation can use the stored water during the distinct dry season, and deeper layers that harbor the groundwater table (Van Asten, 1996). Despite huge costs and efforts, remediation has remained unsuccessful during the past 150 yr.
As the significant amounts of Na in groundwater, the soil types, and the range of groundwater levels in the Hortobágy region are similar to those in simulated cases, an available dataset of the Hortobágy area from the Agricultural University of Debrecen in Hungary (Van Asten, 1996) is suitable to compare with our simulations. However, our aim is not to simulate a local situation as the Hortobágy soils, but to assess SISD for seasonal weather in general.
The dataset used consists of 28 soil profiles with sodic conditions, taken in a toposequence. Each of the profiles consist of several genetic horizons (i.e., horizons that are significantly different from horizons above and below, regarding texture, organic matter content, color, etc.). For each of these horizons, CEC, base saturation, and adsorbed Na are available. Soil profiles were described to a depth of 1.6 m; since the lower boundary of sodic layers was often below this limit, we only define the top of the sodic layer. The soils were classified according to the Russian-Hungarian classification. The dataset comprised Meadow Chernozems (higher in the toposequence and nondegraded), Meadow Solonetz (sodic soils at the middle position in the toposequence), eroded or non-eroded Solonetzic Meadow soils, and a Meadow Soil that is Solonetzic in the subsoil (all in the lowest positions of the toposequence). Groundwater quality was not constant; Na concentrations in the lowest part of the toposequence were lower than those higher up in the sequence. Figure 3 shows the development of the soil solution concentration (C t ) and ESP over time, averaged over the root zone, for the reference precipitation intensity. Both C t and ESP increase faster and obtain a higher value for shallower groundwater depths, as observed from the two loam simulations. Comparing the two soils with a groundwater depth of 200 cm, C t shows a slightly stronger response on dry and wet periods for the loam soil, but buildup of ESP within the root zone is less. All simulations seem to have reached a quasi-steady state in the root zone with respect to C t around the time of the downpour event. Despite the quasi-steady state in C t , ESP still increases, especially for the situations with deeper groundwater levels. Nonetheless, the downpour event causes clear reductions in both C t and ESP (Figure 3). The relative reduction in C t over the root zone after the downpour is more pronounced than for ESP (e.g., roughly 45 vs. 30% for the clay loam simulation). This fact poses a risk of structure degradation, as adsorbed monovalent ions in a low-salinity environment are responsible for SISD. Figure 4 shows C t , ESP, reduction function r 1 , and hydraulic conductivity in the unsaturated zone before (thin lines) and immediately after (thick lines) the downpour event. A peak in C t is evident (Figure 4a) within the root zone, governed by capillary rise of groundwater (supplying salts) and root water uptake (condensing salts). This peak is located highest in the profile for the shallow loam soil, followed by the clay loam. The deep loam soil has its peak located at the deepest position, in correspondence with the lowest average concentrations over the root zone (Figure 3). Also note the differences in width and intensity of the peaks.

Simulation results
Peaks in ESP are visible within the root zone as well (Figure 4b). As high salt concentrations exist within the root zone, valence of ions is less relevant for the process of adsorption, leading to a relative increase of monovalent Na + (thus a higher ESP) as compared with a situation with a lower ionic strength. During the downpour event, a downward movement of infiltrating water pushes peaks in C t downward. Adsorption of divalent Ca 2+ is now preferred in the upper regions of the root zone, resulting in a decrease in ESP where fresh water infiltrates. However, as Ca concentrations in infiltrating water are low, not all Na is replaced instantaneously, resulting in a relatively high ESP in combination with low solution concentrations. This combination, in turn, results in swelling (and possibly dispersion) of clay, and thus the risk of soil structure degradation. This is reflected in reduction function r 1 in Figure 4c. The layer in the root zone in which reduction occurs is the developed sodic layer. Figure 4c shows a 25-cm-thick sodic layer with significant reduction in hydraulic conductivity for the shallow loam soil, as well as a less strongly developed sodic layer in the clay loam soil. No reduction is present in the root zone for the deep loam soil under reference rainfall intensity. This pattern is in correspondence with the root zone average ESP (Figure 3).
The effect of the reduction function upon hydraulic conductivity is evident from Figure 4d, for the shallow loam soil. Most reduction in hydraulic conductivity is, however, caused by nonsaturated conditions (K r in Equation 2a), rather than reduction due to sodicity. However, reduction due to nonsaturated conditions is reversed by rewetting of the soil, whereas reduction due to structure degradation is poorly reversible.
Also notable is the reduction in relative hydraulic conductivity below the root zone. This zone of reduction is mainly dependent on the lower boundary, for which a high Na/Ca ratio is prescribed (Table 3), leading to a high ESP in combination with moderate salt concentrations. As our focus is on sodic layer formation due to a sudden fresh water input, we do not consider this zone any further. Non-smoothness of the reduction functions (Figure 1) used in HYDRUS-1D is evident by the sudden shift in reduction of hydraulic conductivity for the shallow loam (110 cm below surface), deep loam (150 cm below surface), and clay loam (140 cm below surface) soils and cannot currently be avoided with HYDRUS-1D. Figure 5 shows how changing the rainfall intensity of the downpour event affects the main sodic layer properties: the minimum reduction factor r 1 , the thickness of the sodic layer, defined as the thickness of the layer in the root zone where a reduction in hydraulic conductivity occurs, the mean depth of the sodic layer, defined as the depth of the top plus bottom of the sodic layer divided by two, and the minimum ESI in the sodic layer. Each simulated point in Figure 5 with the same soil type and groundwater depth has exactly the same conditions prior to the event. As noted, the reference event intensity (128 mm over 3 d) induces the formation of a sodic layer within the root zone for the shallow loam and clay loam, but this occurs already for less intense events. Reduction occurs from approximately 0.55 and 0.85 times the reference intensity, for the shallow loam and clay loam soil, respectively (Figure 5a), and is remarkably abrupt. The threshold after which reduction occurs for the deep loam is at rainfall intensities exceeding 1.4 times the reference intensity. With intensities exceeding this threshold, reduction in hydraulic conductivity rapidly becomes more severe. Apparent also is that a deeper groundwater table reduces the risk of sodic layer formation due to a decrease of capillary rise of groundwater and corresponding transport of salts. Soil type plays a role as well; given the current parameterization of soils (Table 2), the risk of SISD is larger for the clay loam than the loam soil. The shape of the curves in Figure 5a is comparable with the shape of the second curve in Figure 1, which displays reduction for a solute concentration of 1 mmol c L −1 , as a function of ESP. In fact, the prescribed salt concentration of rainwater in the model (Table 3) corresponds to the situation of this curve in Figure 1. Apparently, as rainfall intensity increases, the infiltration front penetrates deeper into the soil. Here, it encounters a higher ESP (Figure 4b), resulting in a stronger reduction in hydraulic conductivity. Hence, the relation between rainfall intensity and ESP displays a similar curvature as the relation between reduction function r 1 and ESP.
The thickness of the sodic layer ( Figure 5b) increases sharply with an increasing rainfall intensity when it exceeds the threshold intensity, which, in combination with a further reducing hydraulic conductivity, severely hampers water movement through the soil. For the three combinations of soil type and groundwater depth considered, a similar pattern in sodic layer development appears, the main difference being the threshold rainfall intensity at which SISD commences. Close observation of the three simulations reveals very minor differences at which the thickness increases. These originate from differences in the actual infiltration speed of water, but also the distribution of C t and ESP prior to the event, influencing relative hydraulic conductivity. Nonetheless, regardless of the differences in properties of the soils and the deviating conditions prior to the precipitation event, they behave remarkably similarly. Mean depth of the sodic layer with respect to field surface (Figure 5c) increases linearly with event intensity, with a certain initial depth at the onset of formation of a sodic layer, which depends on the C t and ESP profiles at the start of the rainfall event ( Figure 4). Also, slight differences in the angle at which the mean depth increases are observed, similar to the thickness.
Comparison of ESI and r 1 (Figures 5d and 5a, respectively) reveals that reduction of hydraulic conductivity through r 1 starts at an ESI of approximately 0.075, rather than 0.05 as suggested by McKenzie (1998). The shapes of the curves below the threshold ESI are nearly identical (confirmed by correlating the two) for ESI and the r 1 function.  (Van Asten, 1996). Distinction is made between observations with shallow (square), moderate (plus), and deep (cross) groundwater levels. The same relation is also shown for simulated results for the loam (150-to 170-cm groundwater depth, triangle) and clay loam (190-to 200-cm groundwater depth, circle) soils using the reference rainfall intensity

Field observations
Field observations are shown in Figure 6. It shows a clear relation in the ESP found in the sodic layer and the depth of the start of this layer. A distinction is made between profiles higher up in the toposequence (with deep groundwater levels), in the middle positions, and in the lowest positions (with shallow groundwater levels). In the middle positions, the sodic layer is found at a shallower depth and has a higher ESP, as compared with the deep groundwater levels. Comparison with the lowest positions is not possible, since analyses of the groundwater composition has shown that Na concentrations in these positions were lower than for locations with moderate and deep groundwater levels. The same relation between groundwater depth and ESP is found in our simulations; the simulation with the sodic layer closest to the surface (shallow loam soil) also has the highest ESP in the sodic layer (Figure 3). Apart from field observations, Figure 6 also gives five simulated relations between the depth at which the sodic layer starts and ESP. These results were obtained for groundwater depths of 190 and 200 cm, and 150, 160, and 170 cm for the clay loam and loam soil, respectively. The largest simulated groundwater depths yield the lowest ESP and the deepest top of the sodic layer in Figure 6. That the field observations for a comparable, yet slightly different, situation to that in which simulations were performed correspond to the simulations corroborates the validity of the simulations performed.

DISCUSSION
Simulation of two soil types using HYDRUS-1D with major ion chemistry module UNSATCHEM has provided insight into the relation between rainfall event intensity (or irrigation water gift) and the development of a sodic layer, in soils with shallow, Na-rich groundwater. This research confirms the risk of sodification in such situations, hence a sudden fresh water input can result in SISD, and therefore in a decrease in hydraulic conductivity, reducing the potential of the soil for (crop) growth. Striking, however, is the regularity with which the mean depth, thickness, and the peak in hydraulic conductivity reduction can be described as a function of rainfall intensity: these relations appear to be very similar for the three combinations of soil type and groundwater depth described in this study, with the major difference being the rainfall intensity at which SISD starts. Sodicity formation as a function of rainfall intensity thus shows a threshold behavior, regardless of groundwater depth or soil type. This difference in threshold intensity between these different simulations is not related to weather conditions prior to the downpour event, nor to groundwater quality, as these were the same for all simulations. Instead, mainly soil hydraulic properties and groundwater depth determine the potential for capillary rise, and thus the buildup of salinity and sodicity in the unsaturated zone. This corroborates the findings by Schofield et al. (2001) in Hungary that groundwater depth and soil texture seem to determine the occurrence and depth of sodic layers. For other soil types, crops, and historical weather conditions, these relations may be similar as well.
The quality of infiltrating water plays an important role in the relation between rainfall intensity and hydraulic conductivity reduction, as demonstrated by the similarities in shapes between Figures 1 and 5a. Additional simulations (results not shown) in which the quality of the infiltrating water was altered confirm this assessment. Lower Ca concentrations (and thus a lower ionic strength and limited replacement of adsorbed Na) resulted in (a) a decrease in depth at which the sodic layer formed, (b) a lower rainfall intensity at which reduction starts, and (c) a more severe reduction in hydraulic conductivity. This also is in correspondence with laboratory results (Morshedi & Sameni, 2000).
In irrigation practices, a frequently used standard for water quality is the threshold electrolyte concentration, which represents the minimum requirements for irrigation water to ascertain a hydraulic conductivity reduction of, at most, 10-25% (Dang, Bennett, Marchuk, Biggs, & Raine, 2018;Ezlit et al., 2013;Quirk & Schofield, 1955). At this reduction percentage, dispersion of clay is unlikely to occur. Based on the current findings, however, this practice can be treacherous. The steepness of the relations between rainfall intensity and both hydraulic conductivity reduction and sodic layer thickness as found in this study (Figure 5a and 5b) shows that SISD quickly becomes more severe at rainfall intensities (or irrigation gifts) only slightly exceeding the threshold intensity. This is alarming, as upon manifestation of the effects of sodicity (such as a reduced infiltration capacity), the impact of SISD might already be severe and poorly reversible.
To determine the risk of SISD, it is therefore advisable to combine information on, for example, groundwater depth and quality and soil properties with the likelihood of the occurrence of intense precipitation events or flooding by rivers and information on irrigation water quality and quantity. Since intense precipitation events or flooding might be unavoidable, measures to reduce soil ESP by soil amendments such as gypsum, or irrigation strategies that limit capillary rise of Na-rich groundwater, may be needed.
Simulations have also confirmed that the current state of the numerical model limits possibilities to extend simulations to other soil types or different conditions, as required for further risk analysis. Numerical instability was a frequent problem of the simulations and could be handled only with very specific conditions with regard to time and spatial discretization, iteration parameters, boundary conditions, and even output times. Changing one of these factors could result in nonconvergence. This was also noted by Reading et al. (2012), who added some alkalinity to the solution to increase numerical stability; we did not do this in the current study. Regardless of the choices made, UNSATCHEM encountered nonconvergence errors in specific nodes at specific times during all simulations. These were likely the result of sharp fronts in water (quality) propagating downward, affecting pressure heads, precipitation of salts, and hydraulic conductivity.
Model improvements are necessary regarding the reduction function that is currently encoded in HYDRUS-1D. The effect of the discontinuous reduction function r 1 clearly shows in the modeled profile of the reduction function (Figure 4c), and this affects the modeled flow of both moisture and salt. This was also noted by Ezlit et al. (2013) and van der Zee et al. (2014). Reading et al. (2012) suggest making reduction function parameters user accessible or implementing alternative functions. As ESI is much simpler (both conceptually and computationally) than the r 1 function used in HYDRUS-1D, yet produces almost identical results, the ESI might be considered as an alternative. The more elaborate, modified r 1 function presented by Ezlit et al. (2013) might also be implemented for situations with sufficient data availability to permit more complexity.
A conceptual complexity that deserves attention in models is that degradation of soil structure due to sodicity effects is unlikely to be reversed "simply" by a decrease in ESP or increase in soil solution concentration alone. This has been demonstrated with laboratory conditions by Minhas and Sharma (1986), for example. A significant decrease in hydraulic conductivity will, in reality, result in isolation of the A horizon from the groundwater and the destruction of roots within the affected layer. Changes in soil structure due to clay particle displacement and swell-shrink processes take place, which not only affect the saturated hydraulic conductivity, but also result in changes in the pore size distribution. This inevitably affects the retention function parameters and properties such as CEC. These processes are not currently specifically modeled by any model, although the r 1 function is an empirical approximation. Restoring these properties requires the formation of larger pores, which can be achieved by biological activity or tillage. Therefore, an increase in r 1 may not occur spontaneously if sodicity and ionic strength return to more favorable values. The process described is basically the process of horizon formation, which is not well captured by any of the current models.

CONCLUSIONS
This research has confirmed, through model simulations with HYDRUS-1D, that SISD can occur for soils with shallow brackish groundwater due to a sudden fresh water input (for instance, precipitation or irrigation). Simulations show that over time, the percentage of monovalent Na + ions at adsorption sites (ESP) in the root zone increases as a result of capillary rise (supplying salts) and evapotranspiration (condensing salts). Sudden fresh water inputs lead to a fast downward movement of solutes in the root zone, while a high ESP remains. This results in the formation of a sodic layer with a severely decreased hydraulic conductivity. As the formation of this sodic, textural B horizon is poorly reversible, SISD should be avoided. The severity of structure degradation depends on, among other things, the intensity of fresh water input, its ionic strength, groundwater depth, and the initial conditions with regard to ESP within the soil. However, for fresh water input intensities only slightly higher than the threshold intensity at which SISD starts to occur, hydraulic conductivity is already severely reduced. Also, the thickness of the sodic layer increases fast with increasing input intensities as the infiltration water front moves deeper into the soil, hampering water flow even more. Regardless of the soil type or groundwater depth, simulations have shown similar relations in the formation of the sodic layer as function of rainfall intensity, the main difference being the threshold rainfall intensity at which the onset of reduction in hydraulic conductivity takes place.