On several types of hysteresis phenomena appearing in porous media

1.   Introduction

Before it was recognized that hysteresis effects appear in several processes connected to porous media such as soils in fluid mechanics, geotechnical engineering, agronomy and hydrology, the phenomenon appeared and was studied in other fields as, for example, magnetism and ferroelectricity (see e.g., ^[13]) or mechanics ^[65]. It is now clear, and this article is aimed to underline the importance in various fields, that hysteresis in the response of porous media is not just a minor issue examined in a couple of scattered references, but it is nowadays a widely studied phenomenon whose influence is strongly recognized. For geotechnical engineering it was rather identified as the core of many basic processes. So, hysteresis is at the basis of consolidation, is a key feature of the plastic response of clay soils and of the frictional response of granular media. The literature on hysteresis effects, not only in geotechnical engineering, is extremely comprehensive. To name but a few, see e.g., ^{[12,13,29,34,48,50]} which are also mentioned in ^[4].

The general term hysteresis refers to the dependence of a system not only on its current state but also on its history. A certain effect occurs delayed compared to a change in the cause. This means that a quantity for the same input value can take two or more possible output values depending on the prehistory. The process then is path dependent and if two paths are connected a hysteresis loop is built. For a cyclic process possibly further, so-called inner hysteresis loops, arise. The term hysteresis comes from Greek, means to lag behind and was introduced in 1885 by Ewing in a work on magnetism ^[25].

An actual example examined in the basement in the house of the author^* shows the effect vividly and serves here as an introduction to the topic: the house is old and located directly on the banks of the river Rhine. The floor of the vaulted cellar consists of ground simply covered with paving stones. In the last decade several high and low water levels of the river Rhine could be observed. An increase of the Rhine level yields a rising ground water level. Starting from a certain height the groundwater pushes from below into the basement. A table of clear water then appears on the floor. As shown in Figure 1 in 2013, when in the measuring period the highest Rhine water level appeared, the water level in the basement was around 1.50 meters. The phenomenon which may take place not only inside buildings but also in topographic depressions in the open country is called basic stage flood. Of course, the rising and falling of the water in the basement is driven by the amount of water in the river Rhine and accordingly the ground water level. It is clearly visible that the level in the basement lags behind the level in the river and also that a hysteresis loop arises for the stage during in- and outflow of the water.

^*The author married in 2016 and changed her name from Albers to Detmann.

In Figure 1 on the horizontal axis the Rhine water level in Mainz (decisive water level gauge near the house) and on the vertical axis the water level in the basement is given in cm. For the high levels of the Rhine between 2003 and 2024 the respective curves are plotted. It is obvious that hysteresis loops form. At reach of a certain Rhine water level (which is not each time the same – presumably the initial value of the ground water level plays a role) with further increase water enters the basement and its water level also rises. If the water level of the Rhine begins to fall the basement water level first still rises. Only after some time the level also decreases and, finally, there is no more water on the floor. However, until one can say that it is completely dry several months pass. The vertical axis starts with negative values because at the beginning there is not suddenly a water table visible but the water slowly penetrates through the joints between the paving stones. It is not surprising that the curves partly are not smooth because heavy rainfalls during the rise/fall period influence the water level of the Rhine, therewith of the ground water and finally in the basement.

Even if 2011 to 2013 yearly high water levels were present it is noticeable that in the last years high water levels are occurring more frequently. Generally, climate change is increasingly leading to special weather events. For this reason it is interesting and necessary to understand certain processes and their consequences for the environment. Processes in particularly dry or wet zones are far better understood than processes in which the amount of water and air in the pores of soils change. In dry or completely wet zones the pores between the soil particles are filled with a single pore fluid (water or air) and, thus, they are saturated. But in general, both pore fluids are present and an interaction between them, as for example a redistribution of the air by water, takes place.

Such processes in unsaturated porous media which involve flow and the exchange of components are of interest in this work. The study is not limited to hydraulic hysteresis appearing in drying/wetting processes or the stage/discharge hysteresis mentioned above. The aim of the paper is to show several appearances of hysteresis effects in porous media and to show why they occur. In case of wetting/drying of soils the insights are verified by results of experimental tests carried out specifically for this purpose in the geotechnical laboratory of the University of Duisburg-Essen.

2.   Stage/discharge hysteresis

Let us come back to the example shown at the beginning. The phenomenon shown there is known in River Engineering as stage/discharge hysteresis. The determination of the discharge of a river by direct measurements is laborious. It is therefore normal practice to establish a relationship, called rating curve, between the stage at a particular gauging station and the discharge. From observation of the stage, the discharge can then be estimated ^[39]. If the hydraulic channel control is not changing and the stage-discharge relationship is not influenced by unsteadiness, it is relatively easy to derive a single-valued rating curve which reflects the one-to-one correspondence between stage and discharge for steady and uniform flow ^[54].

Unsteady flow, in contrast, is connected with flood wave propagation. In ^[49] is pointed out what happens if a flood wave passes a reach. If the wave front arrives, in the upstream section the velocity of approach is increased. When the flood peak passes into the downstream section, the rear of the wave increases the backwater condition and so reduces the velocity at a given discharge at the measuring cross-section. The result is that on rising discharges the stage at the measuring cross-section generally is lower at a given discharge than when discharges are falling. Thus, when stages are plotted against corresponding discharges measured during the passage of the wave a looped curve is obtained. This can be seen in Figure 2 (right) where $Q$ denotes the discharge, $h$ the stage and $U$ the flow velocity. The situation is comparable to this in the basement. If the Rhine level is very high, also the current is strong and the flow velocity is high. This directly gets obvious by the large amount of floating debris. When the peak is over, the water level decreases and simultaneously the flow velocity reduces. The stage (vertical axis) - discharge (horizontal axis) relationship shown in Figure 2 (right) shows a similar behavior as the Rhine water level - basement water level connection shown in Figure 1. Additionally, on the left hand side the well-known fact is shown that the maximum water discharge arrives before the maximum stage. It should also be mentioned in this context that another process pair goes along with the flow in rivers, namely sedimentation and erosion. If the water flows faster not only a bigger amount of floating debris can be observed but also that particles of the river bed will be transported. If the water after that flow is slower again, then particles will be deposited. Of course, this does not only depend on the flow velocity but also on the geometry of the river.

Free-surface flows in rivers and channels are described by the Saint-Venant Equations which are also called Shallow Water Equations. Under the assumptions that the flow is one-dimensional, unsteady, non-uniform and almost rectilinear, that the bed slope is mild, the water is incompressible and the river width constant, the mass and momentum conservation laws can be formulated in the following way

and

These equations include the two unknown functions $h(x, t)$, the stage which is the vertical distance between free surface $S_{w}$ and bed slope $S_{0} = -(\partial z/\partial x)$ and $U(x, t),$ the velocity averaged over the cross-section $A\left(x, t\right) = Bh,$ with width $B$ assumed to be constant. The discharge $Q = UA,$ $g$ is the earth acceleration and $S_{e}$ the energy slope, also called friction slope. The latter depends on the wall shear stress which is in turn dependent on the flow velocity $U.$ These quantities can be related, e.g., by the Darcy–Weisbach formula (after Henry Darcy [1803–1858] and Julius Weisbach [1806–1871])

which originally was established for uniform, steady flow and contains a friction factor $f$ and the hydraulic radius $R_{h}$.

Another well-known (and earlier introduced) relation is the Chézy formula (after Antoine Chézy [1718–1798])

where $C$ is the Chézy coefficient also considering friction. Robert Manning [1816–1897] showed a better fit to observations if the hydraulic radius has the exponent 2/3 and proposed using the roughness value $n$

If Chézy's version is used, for unsteady flows follows from (2.2)

a non-unique relation with a loop, shown in Figure 2 (right) in black. The width of the loop is an indication of the influence of inertia and pressure terms ^[34,35]. For steady and uniform flows (2.6) simplifies to

which is presented by the unique, gray curve in Figure 2 (right). The relations $Q = Q\left(h\right)$ are referred to as rating or gauging curves at a certain section for unsteady or steady flow (see e.g., ^[39]). According to ^[35] several observations can be summarized. The most important which also qualitatively describe the relation between the water level in the basement and the water level of the Rhine (see Figure 1) in the above mentioned example are:

● For steady flow the discharge $Q$ has a single value. For unsteady flow two different values for the same flow depth, $h$, depending on the question whether the water level is increasing or decreasing occur. The emerging loop has counter-clockwise direction.

● When the discharge in a section reaches its maximum $Q_{\mathrm{max}}$, thus, $\partial Q/\partial t = 0,$ the flow depth, $h$, is still increasing (see Figure 2).

● Since $Q = UBh$, for constant $B$, $\partial Q/\partial t = B\left(U\partial h/\partial t+h\partial U/\partial t\right)$. When $\partial Q/\partial t = 0,$ either the two derivatives cancel themselves or they take opposite signs. If $\partial h/\partial t > 0,$ then $\partial U/\partial t < 0,$ thus the velocity decreases. This implies that $h_{\mathrm{max}}$ arrives after $Q_{\mathrm{max}}$ preceded by $U_{\mathrm{max}}$ (see Figure 2 (left)).

Of course, there are additional effects like heavy rainfalls during the process which let the curves in Figure 1 look not as smooth as in Figure 2.

While the above mentioned relations between velocity and energy slope at the beginning were set up for steady flows, it should be mentioned that already in 1916 Benjamin Jones (born in 1893) dealt also with unsteady flows. He found a hysteretic relation correcting the river discharge for a changing stage ^[41]. The method is still used in actual studies ^[53,54].

3.   Soil-moisture hysteresis or dependence of capillary pressure on saturation change–drying/wetting

The probably most well known hysteresis effect in soils is the dependence of the capillary pressure on changes of the saturation for drying/wetting processes. The saturation may change between full water saturation $S_{w}$ (when the pores are completely filled with water) and full air saturation $S_{a}$ (when the pores are completely filled with air). In nature these limit values will not occur since always residual or irreducible saturations of the pore fluids remain. They arise when the pore fluid is entrapped in the complex pore system and cannot exit completely. For water the order of the residual value is of 15%, for air it is around 5%. In nature such changes of the degree of saturation appear e.g., in drying/wetting or drainage/imbibition processes.

It is well known that often not even cycles of the degree of saturation are necessary to cause failures but that e.g., on slopes the occurrence of very high degrees of saturation may lead already to landslides. Studying artificial rainfalls this topic was investigated, e.g., in ^[37,60] and, particularly, in ^[24] it was shown that hysteresis effects are important also for this problem. In these approaches hysteresis is typically considered by modeling only the secondary scanning curve, while the entire hysteresis cycle is modeled in only a few cases. It is pointed out that experimental techniques which are currently used neglect the hydraulic hysteresis and that they represent the greatest limitation to landslide forecasting.

It was shown already in earlier works (see e.g., ^{[4,5,6,7,8,9,23]}) that for wetting and drying of soils (absorbing and desorbing water) the physical properties are different: If an initially dry soil sample is wetted, a different capillary pressure/saturation curve arises as if a wet sample is drained. If after a first wetting/drying cycle, reflected by main drying/wetting curves, further cycles follow then inner, so-called scanning curves, arise (see left panel of Figure 3).

It is well known that the pressure of fluids depends on the type and and the physical properties of the individual fluid. While water is only slightly compressible and often is regarded to be incompressible, gases like air are characterized by very high compressibility. This leads to a difference in pressure of the pore fluids. This discontinuity in pressure across the interface of the two immiscible fluids is called capillary pressure $p_{c}$ and is defined by the microscopic relation $p_{c} = p_{a}-p_{w},$ where $p_{a}$ and $p_{w}$ are the true pressures in air and water. In hydrology often the pressure head is used instead of the capillary pressure with

containing the true mass density of the fluid $\rho^{FR}$ and the earth acceleration $g$. The water saturation $S_{w}\in [0, 1]$ often is replaced by the water contents $\theta_{w}\in [0, n]$ where $n$ is the porosity defining the pore ratio.

The phenomenon of this hydraulic hysteresis is mostly attributed to a number of causes like the ink-bottle effect which results from the shape of the pore space with interchanging narrow and wide passages, the raindrop effect which is due to the fact that the contact angle at the advancing trace of an interface on a solid differs from that at the receding one or the snap-off effect which occurs if the void space diameter is very different from the void opening diameter (see e.g., ^[18] or ^[59]). In the newer literature also other causes are discussed. In ^[19], a model of pore expansion and contraction caused by hydraulic loading is established to explain hysteresis.

The incorporation of the soil-moisture hysteresis into models is a tricky matter. Mostly in engineering works it is completely ignored that the behavior for drying and wetting is different and only one single curve is used. Several researchers dealt with the theoretical formulation of the capillary pressure/saturation relation of which the most recognized are van Genuchten ^[66] and Brooks/Corey ^[15]. However, these investigations are limited to one branch of the hysteresis curve. But it became clear very quickly that the hysteresis in the mechanical response of soils in not a minor issue but that it is rather the core of a whole discipline, namely geotechnical engineering. One of the first works in which hysteretic two-phase flows were studied experimentally and numerically is ^[46]. Therein, the authors describe how experimental water contents are measured with a gamma radiation attenuation system and experimental water pressures are measured with porous ceramic probes connected to pressure transducers. Simulations are conducted employing a proposed hysteretic model and a corresponding nonhysteretic model. Henceforth, more and more researchers were concerned with the hysteresis effects in soils. In ^[7] the nonhysteretic procedure has been expanded using separately two curves which, however, are connected to each other. Only one of the curves is measured, the other one is predicted using the measured data. The original prediction method has been proposed by Haverkamp et al. ^[38]. Several other researchers dealt with the modeling of the hysteresis effect. In ^[57], a hysteretic hydraulic constitutive model has been proposed in which main wetting and main drying curves are modelled using a modified version of the van Genuchten model. The new version improves the behavior at low degrees of saturation and was originally introduced in ^[26]. Also closed-form expressions for the scanning curves are included. Another water retention model for deformable soils ^[31] takes into account hydraulic hysteresis caused by both wetting-drying cycles at constant void ratio and compression-swelling cycles at constant suction. Main wetting and drying surfaces are defined and, additionally, the scanning behavior inside the domain of attainable soil states is described. An approach which accounts for both the hydraulic hysteresis and the specific volume dependence of the retention relationship in a three-dimensional formulation is proposed in ^[64]. The authors introduce a relationship which models the soil compressibility with suction as a function of the degree of saturation. The mathematically sophisticated approach, of course, is the use of a hysteresis operator. The main problem is the adjustment of the physical data. This issue has been addressed first by a group of researchers around Flynn and led to the works ^[27,28]. Without precise adjustment of physical parameters to a special case, the incorporation of hysteresis operators into soil-moisture models and the mathematical solution was topic of joint works of Krejčí and Albers ^[8,9]. One important point which has to overcome is that models containing hysteresis operators are formulated in Lagrangian description while classically soil moisture and porous media approaches are in Eulerian form. There are first approaches to describe multi-phase models for porous media in Lagrangian way (see e.g., Wilmanski ^[69] or Coussy ^[21]) but the task remains complicated until today.

In the following the above mentioned modeling approaches with two separate curves and with a hysteresis operator are summarized.

3.1.   Engineering-like approach to capture soil-moisture hysteresis

The engineering way to describe soil-moisture hysteresis is to describe the drying and the wetting curve separately. Mostly, it is presumed that both the main drying curve (MDC) and the main wetting curve (MWC) or alternatively the primary wetting curve are measured and scanning curves (i.e., inner curves for further drying/ wetting cycles) are interpolated from the data of both curves. Parlange ^[52], on the contrary, presented an approach to predict the second boundary and scanning curves from only one boundary curve. Since the measurements of these curves are laborious and time consuming this offers advantages compared with other methods. Since the Parlange approach had still some flaws (see e.g., the predicted wetting curve did not have an inflection point) Haverkamp et al. ^[38] introduced an extension which implies that all drying and wetting curves (regardless of the scanning order) have the shape of the van Genuchten equation ^[66] which reads in normalized form

Therein $\theta$ is the volumetric water content, $\theta_{r}$ the residual water content, $\theta_{S}$ the water content at natural saturation, $h$ the soil water pressure head or matric potential, $h_{g}$ a van Genuchten pressure head scale parameter and $n$ and $m$ are van Genuchten shape parameters. In the Haverkamp et al. procedure one of the main curves is fitted to measured data while the other is predicted. This leads to the following equations for the MWC and the MDC

where the subscripts $mw$ and $md$ refer to main wetting and main drying, respectively. The relations between the specific wetting parameters $m_{mw}$, $n_{mw}$ and $h_{gmw}$ and the drying parameters $m_{md}$, $n_{md}$ and $h_{gmd}$ are specified in ^[38] as follows

Exemplarily, in Figure 3 (right) measured data of the MDC, the fitted MDC and the predicted MWC for Del Monte sand (soil parameters are available from ^[47]) are reproduced. The irreducible water contents is neglected here. The approach is used e.g., in ^[7] to analyse the wave propagation for drying and wetting processes in unsaturated soils using the continuum model introduced in ^[3]. The result is that the shear wave and the fastest longitudinal waves are only slightly vary for drying/wetting while the second and third longitudinal waves which are driven by the pore fluids are stronger affected.

3.2.   Mathematical approach to soil-moisture hysteresis using a Preisach operator

More elegant is the incorporation of an operator to describe the hysteresis in closed-form. Such a mathematical treatment of soil-moisture hysteresis based on the Preisach model was proposed by Flynn. In ^[27], Flynn states that for the mathematical treatment of hysteresis first an elementary rate independent hysteresis nonlinearity, a so called hysteron has to be chosen. Second, the complex rate independent hysteresis nonlinearities have to be treated as block-diagrams of hysterons and, finally, identification principles have to be established.

One frequently used type of hysteron is the non-ideal relay – it is the basic idea of the Preisach model ^[56]. It is characterized by its threshold values $\alpha < \beta$ and internal memory state $\eta\left(t\right)$. Its output can take one of the two values 0 or 1 which means that at any moment the relay is either 'switched off' or 'switched on'. The output $y(t) = R_{\alpha, \beta}\left[t_{0}, \eta_{0}\right]x\left(t\right)$, depends on the input $x\left(t\right)$ and on the initial state $\eta_{0}$ which is either 0 or 1. The main assumption made in the Preisach model is that the system can be thought of as a parallel summation of a continuum of weighted non-ideal relays $R_{\alpha, \beta}$, where the weighting of each relay is $\mu\left(\alpha, \beta\right)$. Such a summation can be uniquely represented as a collection of non-ideal relays as points on the two-dimensional half-plane $\Pi = \left\{ \left(\alpha, \beta\right):\beta > \alpha\right\}$, which is also known as the Preisach plane. The output of the Preisach model is then represented by the following formula:

where $p\left(\alpha, \beta\right)$ is an integrable positive function in $\Pi$, called the Preisach density.

Flynn ^[27] studies a family of ordinary differential equations with Preisach hysteresis and uses them to describe the flow of liquids through porous media. First the mass balance equation is considered

where $\rho_{w}$ is the density of water (obviously, for the case of uniform density of water $\rho_{w}$ can be canceled out), $\dot{\theta} = \frac{d\theta}{dt}$, $\theta_{S}$ is the moisture content at natural saturation and $q$ is the net inflow rate of water volume per unit time into a representative elementary volume (REV). It is assumed that the moisture content is related to the matric potential $\psi$ by the Preisach operator, i.e., $\theta\left(\psi\right) = P\left(\psi\right)$. Thus, hysteresis enters Eq (3.5) in the following way

In order to close this equation, a law is needed, relating the rate of the liquid flow through the porous medium to the potential difference. Darcy's law is used which states that the flow past two points in the medium is proportional to the difference in their potentials ($\psi_{1}$ and $\psi_{2}$), and is inversely proportional to the distance, $L$, between them

For simplicity, $L = 1$ m and the conductivity $k_{c}$ is used as proportionality constant, so that

If now one of the potentials is considered as a reference potential, the hysteretic differential equation takes the following form

where $\mathcal{P}$ is the Preisach operator.

In the book chapter of Flynn et al. ^[28] and in his PhD thesis ^[27] three different types of Preisach densities are investigated and compared. One of them, the so-called 'Wedge model' containing a one-parameter density, together with the classical van Genuchten equation (3.1) has proven to be particularly suitable to determine the main wetting and the scanning curves from the main drying curve. For the 'Wedge model' and each of the other models Flynn searched for the best possible parameter values to fit experimental data. He used data given in ^[38]. As can be seen on the example of Del Monte sand in Figure 4, the theoretical predictions fit well the experimental curves.

4.   Capillary rise/infiltration–drying/wetting

Recently, the one-dimensional infiltration and capillary rise in initially dry uniform sand samples was investigated ^[22]. The focus of that work was the understanding of the two processes arising if water enters a soil sample: Capillary rise if the sample is watered from the bottom and infiltration if the water penetrates from the top. For a constant amount of water during a single passage both processes are, of course, not hysteretic. But if the focus lies on the flow of a certain amount of water through a sample, then a certain water contents occurs two times in a wetting/drying or drying/wetting process. To examine this more closely, a flow test of an initially dry sand sample has been carried out^*. As shown in Figure 5 (left) an initially dry sample of fine sand put into a cylinder of 1 m length with sensors measuring the water contents at 0.1, 0.3, 0.5, 0.7 and 0.9 meters was infiltrated from above. The 5 sensors are ECH₂O EC-5 Volumetric Water Content sensors from METER Group with 5 cm sensor length and they are connected to an Em5b datalogger. In the simple test only one flow direction was investigated and suctions have not been monitored. The water entered the sample from a reservoir at the top in which a constant water level was held. After the entire sample was wetted, the supply was turned off and the the sample began to dry. It gets obvious from the middle panel of Figure 5 that the wetting period with around 30 minutes was rather short while the drying progress depended very much on the depth. The higher the sensor, the faster the region became drier. However, it gets obvious, that for each depths a certain residual water contents remained. For the top it is of order between 0.05 and 0.15, at the bottom the sample hardly dried at all (and this also after a period of three weeks). The situation after three weeks time is shown in the right panel of Figure 5. It should be mentioned that for the qualitative result it is of interest whether the sample is wetted from top or from the bottom. In ^[10], it is pointed out that besides differences of laboratory and field characterizations also differences resulting from different wetting procedures occur. Differences are interpreted as a result of different hysteretic paths being observed.

^*Thanks to Jörg Nolzen and Bjarne Kroll from the Geotechnical Engineering Laboratory of the University of Duisburg-Essen for preparing and performing several tests of this kind.

More interesting for the present article is, that in this process a directly obvious hysteresis effect appears. If one looks at the water contents in the different depths at certain instants of time during the drying process (see left hand side of Figure 6) the curves for an instant of time at the beginning of the drying process and a later instant of time form a loop. The same was observed for the reverse process. In ^[33] an initially saturated sample was studied which first was drained and then rewet. In Figure 6 (right) for two instants of time during the rewetting process the water content over the depth is presented and exhibits also a loop of the curves. The authors solve the problem also numerically. To this aim the Richards equation for one-dimensional vertical flow in pressure head form, as given in ^[33], is solved

In this equation $C$ denotes the water capacity which is a hysteretic function of the soil water pressure head $h$ and the position $z$. The hydraulic conductivity $K$ is a function of water content $\theta$ and position $z$. Based on ^[32] the hysteretic $\theta\left(h\right)$ relation was represented by an empirical equation of the form

where $\theta_{0},$ $b,$ $h_{0}$ and $\theta_{r}$ curve fitting parameters and the hydraulic conductivity-water content relation is given by

with empirical parameters $a$ and $\eta.$ In ^[22] a more recent predictor-corrector method to solve Richards equation had been mentioned and applied. On the webpage ^[14] the matlab code developed to this solution method is downloadable.

5.   Adsorption/desorption

Also in other fields of civil engineering than geotechnics the moisture movement through porous media is examined. In ^[55] adsorption and desorption curves are determined by use of the instantaneous profile method. In experiments initially dry material (calcium silicate) is exposed nearly to saturation and then dried to about 8%. Similarly to the above mentioned experiments in the sand column it is observed how quickly the moisture front moves. Therefore, a chamber is equipped with five vertically placed sensors which measure the water content. While in the time dependence of the relative humidity shown on the left hand side and in the middle of Figure 7 the hysteresis of these processes is not directly palpable, on the right hand side where the moisture conductivity for both processes is given in dependence on the relative humidity the hysteresis is clearly visible. It is remarkable that the moisture conductivity as function of the water content [m³/m³] is found to be non-hysteretic. But this observation is not new, it only confirms findings reported in ^[33].

The notions adsorption and desorption describe the sticking and detachment of particles on/of surfaces (see e.g., ^[36]). In geophysics they are often used to describe the pollution transport with water in soils. In ^[1], a continuum model for adsorption-diffusion processes in soils has been introduced. However, since in uniform, cohesionless soil types generally low concentrations and single-layer adsorption are presumed, in this model a Langmuir-type isotherm without hysteresis is used. Hysteresis only appears if capillary condensation in fine pores and multi-layer adsorption is taken into account. Similarly to the hysteresis of the capillary pressure/saturation curve in engineering approaches concerning soils hysteresis is rarely considered.

The adsorption/desorption behavior is described by adsorption isotherms. They relate the amount of adsorbed material to pressure or concentration. For a constant temperature they reflect equilibrium conditions. The shape of the isotherm is significantly influenced by the particle size. This, of course, is related to the porosity, i.e., the pore ratio. The predominant factor for the adsorbability is the available area. Since in porous media each particle offers space for the process, the enormous inner surface offers best conditions.

Probably the first who set up a formula for an isotherm describing gas adsorption on homogeneous surfaces was Freundlich in 1909 ^[30]. He used the empirical formula

to relate the amount of adsorbed gas $q$ to the partial pressure of the adsorbate $p$ and contains constants $a$ and $n > 1.$ However, the function is not adequate at high pressures because it is well known that the amount adsorbed possesses an asymptotic value for high pressures which is not reflected by (5.1). This was accounted for by Langmuir in his extensive work on adsorption processes around 1915 (exemplarily see ^[44]). It resulted in the very famous adsorption isotherm which under certain assumptions is equivalent to Freundlich's isotherm (see ^[2]). Langmuir assumed that the particles only in a very thin layer (a mono-layer) stick to the surface of the adsorbent. Under the assumption that the particles to be adsorbed do not interact with each other and build together a new layer at specified sites on the surface the following isotherm holds

Therein $q$, again, is the amount of adsorbed material, $p$ is partial pressure of the adsorbate, $q_{0}$ the maximum amount of adsorbate and $b$ is a constant exponentially related to the heat of adsorption. As gets obvious from Figure 8 (left) which gives a qualitative impression of the shapes of adsorption isotherms, in this formulation the amount adsorbed is limited and tends to an asymptotic value (the full occupation of available sites on the mono-layer). Although Freundlich's isotherm (5.1) is not restricted to mono-layer adsorption due to the low pressures for which it is valid, it only holds for a low degree of coverage.

In 1938 the multi-layer theory of Brunauer, Emmett and Teller ^[17] which is based on Langmuir's considerations was published. The graphic representation of

is a S-shaped isotherm. Therein, $v$ denotes the total adsorbed volume, $v_{m}$ the adsorbed volume if the first layer is adsorbed, $p_{0}$ is a reference pressure and $c > 1$ a constant. In ^[16], several isotherms of this type are analyzed and compared to the capillary condensation theory. If the first layer is completely filled and multi-layer adsorption starts the surface becomes inhomogeneous and the passages are not equally wide. Depending on the width and shape, i.e., on the geometry of the pore space the pressure changes and yields adsorption and desorption curves which do not coincide. A detailed explanation for the occurrence of hysteresis and an assignation of the very different types of hysteresis loops to various materials (depicted on the right hand side of Figure 8) can be found in ^[63].

6.   Freezing/thawing

Freezing/thawing processes in soils occur on the one hand under natural conditions in regions where the water in the pores undergoes a phase change between liquid and ice. Additionally, they occur when using ground freezing techniques which are temporarily employed in construction works to stabilize the subsoil. In ^[67], an iterative procedure for the calculation of mechanical properties during thawing/freezing in saturated two-component media is shown. The focus was on the damage which can occur due to cryo-swelling, i.e., volume changes caused by the transition from water to ice. Already in 1966 Koopmans and Miller in ^[43] experimentally showed that analogously to drying/wetting processes also freezing/thawing curves exist and that they in a very similar way exhibit hysteresis. Exemplarily, the curves for a soil considered to be free of colloidal material from ^[43] are shown in Figure 9 (left).

In the meantime also the theoretical description has been tackled. In ^[70], the formulation and calibration of freezing and thawing curves on the basis of the well-known soil moisture curves has been carried out. One of the investigated approaches is the van Genuchten approach which was used above to describe the hysteresis of the drying/wetting curves. Analogously to the above described approach to obtain the wetting curve from the drying curve, here the thawing curve is obtained from the freezing curve. The shapes of the curves follow from the van Genuchten approach. Instead of the capillary pressure, here the temperature is decisive. Thus, Eq (3.1) becomes

where $t$ is the temperature and $\alpha_{f},$$n_{f}$ and $m_{f}$ are van Genuchten parameters for freezing/thawing defined as follows

The pairs $(t_{1}, \overline{\theta}_{1})$ and $(t_{2}, \overline{\theta}_{2})$ are two points on the measured thawing curve. The parameter $m_{f}$ in ^[70] is assumed to be the same for the thawing and the freezing curve. Using the data of ^[43] (Figure 9, left panel) ^[70] fitted the freezing and thawing curves presented in Figure 9 (right). It should be mentioned that depending on the soil type, the initial water content and other factors the hysteresis curves may have another shape ^[62].

7.   Poro-elastoplasticity

Above several hysteresis effects in soils have been shown which are primarily connected to the behavior and physical characteristics of the pore fluids. However, also the solid properties can cause a hysteresis loop in the constitutive law. Although in classical soil mechanics often elastic behavior of the solid is presumed, it is well known that the region of elastic behavior is very small and soils in principle behave mostly irreversibly, e.g., according to a plastic constitutive law. Several material models have been developed particularly to describe the complex behavior of soils (see e.g., hypoplasticity ^[42] or the hardening soil model ^[58]). These models intend to consider features like densification and dilatation, i.e., decrease or expansion of the pore volume, stress dependent stiffness, stress history and cyclic behavior, plastic yielding or hardening. For general information on poromechanics, see e.g., ^[20].

In a triaxial test of a soil sample it can be shown that the behavior in general for the first loading, unloading and reloading is different and a hysteresis loop in the strain-stress relation occurs (Figure 10, left). If further cycles of unloading and reloading are observed, one possible behavior is that a shakedown of the elastic-plastic structure occurs, i.e., the hysteresis loops for further cycles become increasingly narrow and finally match and represent a steady state (accumulation of strain). For extensive studies about the behavior of soils under monotonic and cyclic loading, see e.g., ^{[45,51,61,68]}.

The hysteretic behavior of soils gets even more obvious if the shear stress over the shear strain is qualitatively plotted for a strain-controlled undrained cyclic simple shear test (Figure 10, right). However, it has to be considered that it is an undrained test which means that in the pores fluid exists which, of course, influences the result. Thus, it is an effect which can not only traced back to the properties of the solid. In porous media the components influence each other and some effects cannot clearly attributed to one single phase.

In ^[8,9] the authors considered two types of hysteresis in unsaturated porous media: the soil moisture hysteresis of the capillary pressure and, moreover, the elastoplastic response of the solid. Two independent hysteresis operators describing hysteresis phenomena in solid and pore fluids have been introduced and incorporated into a model for unsaturated soils and mathematically analyzed.

8.   Conclusions and outlook

In this article several types of hysteresis appearing in soils and porous materials have been shown. Mostly they are connected to the behavior of the pore fluid(s) but also the solid component exhibits hysteretic reactions. Hysteretic features occurring in many natural processes and engineering applications are driven by the properties and interaction of the components of a porous medium. In soils their characteristics strongly depends on the soil and pore fluid types, on the geometry of the material and the process considered.

The modeling of hysteresis effects in engineering approaches resulted mostly from empirical observations. Thus, it is no wonder that the theoretical description mostly focused on one of the two branches of a hysteresis loop in which the processes often can be splitted. Often for opposing processes like loading/unloading, drying/wetting, adsorption/desorption, freezing/thawing only one direction is studied. This means, for simplicity, the hysteretic behavior is neglected. If both branches are studied, then they are investigated separately. The mathematical description by use of operators is, of course, much more elegant. However, sometimes the more skillful representation comes at the expense of accuracy. Often, in the more theoretical treatises the processes and their results are not as precisely reflected as is the case in engineering approaches.

In times of challenging climate changes the description of such processes becomes more and more necessary. It is to be expected that processes as described above will increasingly take place. In order to be prepared for the environmental effects it will not be enough to reproduce results of experiments or natural occurrences but also to predict what will happen in the future. To this aim the investigation and theoretical description of above described processes is an important concern.

Use of Generative-AI tools declaration

The author declares she has not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

The author declares no conflict of interest.