Application of the Preisach model in soil-moisture hysteresis
Application of the Preisach model in soil-moisture hysteresis∗ Denis Flynn,
Hugh McNamara,
Philip O’Kane,
Alexei Pokrovskii
Abstract The classical Preisach Model is applied to describe quantitatively the hysteresis effects in the relation between water retention and soil-moisture tension. Special, one parameter, classes of Preisach operators are proposed to construct models of the soil-moisture hysteresis for particular soils. The results of fitting a parameter to imitate experimental data are presented and discussed. Key words: soil-moisture hysteresis, Preisach model, terrestrial hydrology
∗
This study was partially supported by the Boole Centre for Research in Informatics (BCRI).
1
Contents 1
2
Introduction
3
1.1
Hysteresis phenomena . . . . . . . . . . . . . . . . . . . .
3
1.2
Soil-Moisture Hysteresis . . . . . . . . . . . . . . . . . . .
5
1.3
The paper structure . . . . . . . . . . . . . . . . . . . . . .
6
The Preisach Model 2.1
2.2
2.3
3
Basic definitions
7 . . . . . . . . . . . . . . . . . . . . . . .
7
2.1.1
Nonideal relay . . . . . . . . . . . . . . . . . . . .
7
2.1.2
Preisach operators . . . . . . . . . . . . . . . . . .
8
2.1.3
Geometrical description of the evolution rules . . . .
9
Expressions for wetting and drying curves . . . . . . . . . .
11
2.2.1
Main drying curve . . . . . . . . . . . . . . . . . .
11
2.2.2
Primary wetting curve . . . . . . . . . . . . . . . .
13
2.2.3
Secondary drying curve . . . . . . . . . . . . . . .
14
Special classes of Preisach densities . . . . . . . . . . . . .
15
2.3.1
The uniform density . . . . . . . . . . . . . . . . .
15
2.3.2
Wedge densities . . . . . . . . . . . . . . . . . . .
17
2.3.3
Beak densities . . . . . . . . . . . . . . . . . . . .
19
Results of fitting
21
3.1
Soils classification . . . . . . . . . . . . . . . . . . . . . .
21
3.2
Fitting procedure . . . . . . . . . . . . . . . . . . . . . . .
22
3.3
Plots . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
23
3.3.1
Zero-parameter model plots . . . . . . . . . . . . .
23
3.3.2
”Wedge” Model Plots . . . . . . . . . . . . . . . .
25
3.3.3
“Beak” model plots . . . . . . . . . . . . . . . . . .
35
Accumulated table for the fitting accuracy . . . . . . . . . .
46
3.4 4
Discussion
48
5
Appendix: Some Mathematica Programs
49
2
1 Introduction 1.1 Hysteresis phenomena The word hysteresis means etymologically ‘to come behind’. It was introduced into the scientific vocabulary about 120 years ago by the Scottish physicist Alfred Ewing as follows “when there are two quantities M and N , such that cyclic variations of N cause cyclic variation of M , then if the changes of M lag behind those of N , we may say that there is hysteresis in the relation of M and N ” [3]. In modern science it is the collective name for a class of strongly nonlinear natural phenomena. (The wording ‘strongly nonlinear’ means that linearisation cannot encapsulate the observed phenomena.) Hysteresis occurs in many engineering, physical and economic systems. For example, in the ICT area hysteresis: “stores” information in dynamic memory; allows neural-networks to learn; affects the orbits of communication satellites via ferro-magnetic friction. The principal role of hysteresis in the post-Keinessian economics is splendidly highlighted in [2]. Omitting a few technicalities, hysteresis nonlinearities are formally defined as deterministic, rate independent operators [7]. Rate independence means that hysteresis operators are invariant with respect to the action of the group G of affine transformations of the time-scale. That is, the blue and the red pairs of functions belong simultaneously to the totality of possible input/output correspondence. Input x(t)
x(at+b)
Time Output y(t)
y(at+b)
Time F IGURE 1: Rate-independence
3
The mathematical concept of the hysteresis operator was suggested by M. Krasnosel’skii and his co-workers in the seventies [8]. The structure of these concepts is often the following: 1. Choosing elementary hysteresis nonlinearities, so-called hysterons. 2. Treating complex hysteresis nonlinearities as block-diagrams of hysterons. 3. Establishing identification principles. Nowadays the above approach to hysteresis contains a variety of ‘branches’, depending on the choice of hysterons in item (1) and/or the basic type of the blockdiagrams in item (2). The most important is the Preisach model which embraces the phenomena mentioned as the examples in the overview section, as well as many others (see, for example, [13]). In the Preisach model the role of hysterons is played by nonideal relay nonlinearities. Block-diagrams are the standard parallel connection of a number of hysterons. Essentially, the Preisach model was first suggested over 60 years ago [14]. The modern explanation of the pivotal role of the Preisach nonlinearity is based on an operational approach and on special identification theorems. The idea of this theorem is the following. In many physical situations, experiments suggest that the hysteresis relation exhibits the return point memory: every minor loop returns back to its starting point. Moreover, all closed loops with endpoints (x1 , θ1 ), (x2 , θ2 ) and (x1 , θ10 ), (x2 , θ20 ) are congruent.
F IGURE 2: Return point memory
By a classical result of Mayergoyz [10], the output y(t) of a hysteresis with return point memory and congruent loops can be represented as Preisach nonlinearity with a suitable measure µ. 4
1.2 Soil-Moisture Hysteresis In this paper we are interested in the so-called soil-moisture hysteresis. This hysteresis is important in terrestrial hydrology. The essence of this phenomenon is that “Less mechanical work is required to insert water into unsaturated soil than to remove it - good for agriculture!” ([5]). It has been shown experimentally that there is a hysteresis effect in the relation between water retention and soil-moisture tension as far back as 1941 [15]. The hysteresis effect is evident in the Soil-Water Characteristic Curve (SWCC) , in that it is different depending on whether it is for wetting or drying. The origin of the hysteretic effect may be attributed to several factors [5]: 1. Geometric nonuniformity of individual pores, resulting in the so called “ink bottle effect”. 2. Different spatial connectivity of pores during drying and wetting. 3. Variations in liquid-solid contact angle. 4. Air entrapment. Important and successful mathematical descriptions of soil-moisture hysteresis have been suggested. We will just refer to the fundamental papers [1, 6, 11, 12] and the further references therein. To describe the SWCC we require an appropriate equation. The equation we choose for our study is the van-Genuchten equation, which takes the following form: · θ = F (h) = θr + (θs − θr ) 1 +
µ
h hg
¶n ¸−m (1)
3
where θ is the volumetric water content [ LL3 ]; θr is the residual water content (in most cases this is zero); θs is the volumetric water content at natural saturation which is chosen as the water content scale parameter. h [L] is the soil water pressure head (also known as the matric potential), this variable is taken to be negative and expressed in cm of water; and hg is the van Genuchten pressure head scale parameter. The two dimensionless water retention shape parameters, m and n, are related by: m = 1 − 1/n. 5
All of these parameters are taken from [4], where experimental and field data was successfully fitted with Van Genuchten type functions. In this paper we will build on the work in [4] by using the Van Genuchten type Main Drying curves as the basis for a Preisach model, which we will examine quantitatively and compare with [4] in terms of agreement with data. Note that the Van Genuchen form for the curves discussed is not the only possible equation to describe them, we would refer the reader to [4] and the bibliography therein for more details. Main drying
curve with primary
and secondary
scanning
curves
1 MDC
cm3 cm3
0.8
SDC 0.6
Water Content
0.4 PWC 0.2
0 -1750
-1500
-1250
-1000 Water Pressure
-750 Head h
-500
-250
0
cm
F IGURE 3: Schematic diagram of the hysteresis model with the main drying curve (MDC), a primary wetting curve (PWC) and a secondary drying curve (SDC) branching off the PWC.
1.3 The paper structure We will discuss in some detail the Preisach model in Section 2. In particular Subsection 2.1 is devoted to the related analytical and geometrical construction (for a more rigorous discussion see e.g. [8, 10]). In Subsection 2.2 we present explicit expressions for analogs to the wetting and drying curves in the Preisach model context. The central place in the paper is occupied by Subsection 2.3. Here we introduce some special, one parameter, classes of Preisach operators, which will be used afterwards as models of the soil-moisture hysteresis for particular soils. The results of fitting are presented in Section 3. In Section 4 we discuss some advantages and disadvantages of the Preisach model as a description of the soilmoisture hysteresis. In Appendix we relegated some Mathematica programs. 6
2 The Preisach Model 2.1 Basic definitions 2.1.1 Nonideal relay Consider now the nonideal relay. It is characterised by its threshold values α < β and is the simplest hysteretic transducer. Its output y(t) can take one of two values 0 or 1: at any moment the relay is either ‘switched off’ or ‘switched on’. The dynamics of the relay are described by the picture in Figure 4. y
1
α
β
0
x
F IGURE 4: The nonideal relay
The variable output y(t) y(t) = Rα,β [t0 , η0 ]x(t),
t ≥ t0 ,
(2)
depends on the variable input x(t) (t ≥ t0 ) and on the initial state η0 . Here the input is an arbitrary continuous scalar function; η0 is either 0 or 1. The scalar function y(t) has at most a finite number of jumps on any finite interval t0 ≤ t ≤ t1 . The output behaves rather ‘lazily’: it prefers to be unchanged, as long as the phase pair (x(t), y(t)) belongs to the union of the bold lines in the picture above. The values of function (2) at a moment t are defined by the following explicit formula:
7
η0 , 1, y(t) = Rα,β [t0 , η0 ]x(t) =
0,
if α < x(τ ) < β for all τ ∈ [t0 , t]; if there exists t1 ∈ [t0 , t] such that x(t1 ) ≥ β, x(τ ) > α for all τ ∈ [t1 , t]; if there exists t1 ∈ [t0 , t] such that x(t1 ) ≤ α, x(τ ) < β for all τ ∈ [t1 , t].
The equalities y(t) = 1 for x(t) ≥ β and y(t) = 0 for x(t) ≤ α always hold for t ≥ t0 . 2.1.2 Preisach operators The main assumption made in the Preisach model is that the system can be thought of as a parallel summation of a continuum of such weighted non ideal relays Rαβ , where the weighting of each relay is µ(α, β). Such a summation can be uniquely represented as a collection of nonideal relays as points on the two-dimensional half-plane Π = {(α, β) : β > α} (see [8]), which is also known as the Preisach plane which is typically shown as in Figure 5. Here the colored area S = S(t) is the set of the threshold values (α, β) for which the corresponding relays Rαβ are switched on at a given moment t. The output of the Preisach model is then represented by the following formula: Z y(t) =
Z p(α, β)Rα,β [t0 , η0 (α, β)]x(t) dαdβ =
α
Hugh McNamara,
Philip O’Kane,
Alexei Pokrovskii
Abstract The classical Preisach Model is applied to describe quantitatively the hysteresis effects in the relation between water retention and soil-moisture tension. Special, one parameter, classes of Preisach operators are proposed to construct models of the soil-moisture hysteresis for particular soils. The results of fitting a parameter to imitate experimental data are presented and discussed. Key words: soil-moisture hysteresis, Preisach model, terrestrial hydrology
∗
This study was partially supported by the Boole Centre for Research in Informatics (BCRI).
1
Contents 1
2
Introduction
3
1.1
Hysteresis phenomena . . . . . . . . . . . . . . . . . . . .
3
1.2
Soil-Moisture Hysteresis . . . . . . . . . . . . . . . . . . .
5
1.3
The paper structure . . . . . . . . . . . . . . . . . . . . . .
6
The Preisach Model 2.1
2.2
2.3
3
Basic definitions
7 . . . . . . . . . . . . . . . . . . . . . . .
7
2.1.1
Nonideal relay . . . . . . . . . . . . . . . . . . . .
7
2.1.2
Preisach operators . . . . . . . . . . . . . . . . . .
8
2.1.3
Geometrical description of the evolution rules . . . .
9
Expressions for wetting and drying curves . . . . . . . . . .
11
2.2.1
Main drying curve . . . . . . . . . . . . . . . . . .
11
2.2.2
Primary wetting curve . . . . . . . . . . . . . . . .
13
2.2.3
Secondary drying curve . . . . . . . . . . . . . . .
14
Special classes of Preisach densities . . . . . . . . . . . . .
15
2.3.1
The uniform density . . . . . . . . . . . . . . . . .
15
2.3.2
Wedge densities . . . . . . . . . . . . . . . . . . .
17
2.3.3
Beak densities . . . . . . . . . . . . . . . . . . . .
19
Results of fitting
21
3.1
Soils classification . . . . . . . . . . . . . . . . . . . . . .
21
3.2
Fitting procedure . . . . . . . . . . . . . . . . . . . . . . .
22
3.3
Plots . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
23
3.3.1
Zero-parameter model plots . . . . . . . . . . . . .
23
3.3.2
”Wedge” Model Plots . . . . . . . . . . . . . . . .
25
3.3.3
“Beak” model plots . . . . . . . . . . . . . . . . . .
35
Accumulated table for the fitting accuracy . . . . . . . . . .
46
3.4 4
Discussion
48
5
Appendix: Some Mathematica Programs
49
2
1 Introduction 1.1 Hysteresis phenomena The word hysteresis means etymologically ‘to come behind’. It was introduced into the scientific vocabulary about 120 years ago by the Scottish physicist Alfred Ewing as follows “when there are two quantities M and N , such that cyclic variations of N cause cyclic variation of M , then if the changes of M lag behind those of N , we may say that there is hysteresis in the relation of M and N ” [3]. In modern science it is the collective name for a class of strongly nonlinear natural phenomena. (The wording ‘strongly nonlinear’ means that linearisation cannot encapsulate the observed phenomena.) Hysteresis occurs in many engineering, physical and economic systems. For example, in the ICT area hysteresis: “stores” information in dynamic memory; allows neural-networks to learn; affects the orbits of communication satellites via ferro-magnetic friction. The principal role of hysteresis in the post-Keinessian economics is splendidly highlighted in [2]. Omitting a few technicalities, hysteresis nonlinearities are formally defined as deterministic, rate independent operators [7]. Rate independence means that hysteresis operators are invariant with respect to the action of the group G of affine transformations of the time-scale. That is, the blue and the red pairs of functions belong simultaneously to the totality of possible input/output correspondence. Input x(t)
x(at+b)
Time Output y(t)
y(at+b)
Time F IGURE 1: Rate-independence
3
The mathematical concept of the hysteresis operator was suggested by M. Krasnosel’skii and his co-workers in the seventies [8]. The structure of these concepts is often the following: 1. Choosing elementary hysteresis nonlinearities, so-called hysterons. 2. Treating complex hysteresis nonlinearities as block-diagrams of hysterons. 3. Establishing identification principles. Nowadays the above approach to hysteresis contains a variety of ‘branches’, depending on the choice of hysterons in item (1) and/or the basic type of the blockdiagrams in item (2). The most important is the Preisach model which embraces the phenomena mentioned as the examples in the overview section, as well as many others (see, for example, [13]). In the Preisach model the role of hysterons is played by nonideal relay nonlinearities. Block-diagrams are the standard parallel connection of a number of hysterons. Essentially, the Preisach model was first suggested over 60 years ago [14]. The modern explanation of the pivotal role of the Preisach nonlinearity is based on an operational approach and on special identification theorems. The idea of this theorem is the following. In many physical situations, experiments suggest that the hysteresis relation exhibits the return point memory: every minor loop returns back to its starting point. Moreover, all closed loops with endpoints (x1 , θ1 ), (x2 , θ2 ) and (x1 , θ10 ), (x2 , θ20 ) are congruent.
F IGURE 2: Return point memory
By a classical result of Mayergoyz [10], the output y(t) of a hysteresis with return point memory and congruent loops can be represented as Preisach nonlinearity with a suitable measure µ. 4
1.2 Soil-Moisture Hysteresis In this paper we are interested in the so-called soil-moisture hysteresis. This hysteresis is important in terrestrial hydrology. The essence of this phenomenon is that “Less mechanical work is required to insert water into unsaturated soil than to remove it - good for agriculture!” ([5]). It has been shown experimentally that there is a hysteresis effect in the relation between water retention and soil-moisture tension as far back as 1941 [15]. The hysteresis effect is evident in the Soil-Water Characteristic Curve (SWCC) , in that it is different depending on whether it is for wetting or drying. The origin of the hysteretic effect may be attributed to several factors [5]: 1. Geometric nonuniformity of individual pores, resulting in the so called “ink bottle effect”. 2. Different spatial connectivity of pores during drying and wetting. 3. Variations in liquid-solid contact angle. 4. Air entrapment. Important and successful mathematical descriptions of soil-moisture hysteresis have been suggested. We will just refer to the fundamental papers [1, 6, 11, 12] and the further references therein. To describe the SWCC we require an appropriate equation. The equation we choose for our study is the van-Genuchten equation, which takes the following form: · θ = F (h) = θr + (θs − θr ) 1 +
µ
h hg
¶n ¸−m (1)
3
where θ is the volumetric water content [ LL3 ]; θr is the residual water content (in most cases this is zero); θs is the volumetric water content at natural saturation which is chosen as the water content scale parameter. h [L] is the soil water pressure head (also known as the matric potential), this variable is taken to be negative and expressed in cm of water; and hg is the van Genuchten pressure head scale parameter. The two dimensionless water retention shape parameters, m and n, are related by: m = 1 − 1/n. 5
All of these parameters are taken from [4], where experimental and field data was successfully fitted with Van Genuchten type functions. In this paper we will build on the work in [4] by using the Van Genuchten type Main Drying curves as the basis for a Preisach model, which we will examine quantitatively and compare with [4] in terms of agreement with data. Note that the Van Genuchen form for the curves discussed is not the only possible equation to describe them, we would refer the reader to [4] and the bibliography therein for more details. Main drying
curve with primary
and secondary
scanning
curves
1 MDC
cm3 cm3
0.8
SDC 0.6
Water Content
0.4 PWC 0.2
0 -1750
-1500
-1250
-1000 Water Pressure
-750 Head h
-500
-250
0
cm
F IGURE 3: Schematic diagram of the hysteresis model with the main drying curve (MDC), a primary wetting curve (PWC) and a secondary drying curve (SDC) branching off the PWC.
1.3 The paper structure We will discuss in some detail the Preisach model in Section 2. In particular Subsection 2.1 is devoted to the related analytical and geometrical construction (for a more rigorous discussion see e.g. [8, 10]). In Subsection 2.2 we present explicit expressions for analogs to the wetting and drying curves in the Preisach model context. The central place in the paper is occupied by Subsection 2.3. Here we introduce some special, one parameter, classes of Preisach operators, which will be used afterwards as models of the soil-moisture hysteresis for particular soils. The results of fitting are presented in Section 3. In Section 4 we discuss some advantages and disadvantages of the Preisach model as a description of the soilmoisture hysteresis. In Appendix we relegated some Mathematica programs. 6
2 The Preisach Model 2.1 Basic definitions 2.1.1 Nonideal relay Consider now the nonideal relay. It is characterised by its threshold values α < β and is the simplest hysteretic transducer. Its output y(t) can take one of two values 0 or 1: at any moment the relay is either ‘switched off’ or ‘switched on’. The dynamics of the relay are described by the picture in Figure 4. y
1
α
β
0
x
F IGURE 4: The nonideal relay
The variable output y(t) y(t) = Rα,β [t0 , η0 ]x(t),
t ≥ t0 ,
(2)
depends on the variable input x(t) (t ≥ t0 ) and on the initial state η0 . Here the input is an arbitrary continuous scalar function; η0 is either 0 or 1. The scalar function y(t) has at most a finite number of jumps on any finite interval t0 ≤ t ≤ t1 . The output behaves rather ‘lazily’: it prefers to be unchanged, as long as the phase pair (x(t), y(t)) belongs to the union of the bold lines in the picture above. The values of function (2) at a moment t are defined by the following explicit formula:
7
η0 , 1, y(t) = Rα,β [t0 , η0 ]x(t) =
0,
if α < x(τ ) < β for all τ ∈ [t0 , t]; if there exists t1 ∈ [t0 , t] such that x(t1 ) ≥ β, x(τ ) > α for all τ ∈ [t1 , t]; if there exists t1 ∈ [t0 , t] such that x(t1 ) ≤ α, x(τ ) < β for all τ ∈ [t1 , t].
The equalities y(t) = 1 for x(t) ≥ β and y(t) = 0 for x(t) ≤ α always hold for t ≥ t0 . 2.1.2 Preisach operators The main assumption made in the Preisach model is that the system can be thought of as a parallel summation of a continuum of such weighted non ideal relays Rαβ , where the weighting of each relay is µ(α, β). Such a summation can be uniquely represented as a collection of nonideal relays as points on the two-dimensional half-plane Π = {(α, β) : β > α} (see [8]), which is also known as the Preisach plane which is typically shown as in Figure 5. Here the colored area S = S(t) is the set of the threshold values (α, β) for which the corresponding relays Rαβ are switched on at a given moment t. The output of the Preisach model is then represented by the following formula: Z y(t) =
Z p(α, β)Rα,β [t0 , η0 (α, β)]x(t) dαdβ =
α
