Fit soil water retention function with hysteresis

To begin, install Python 3 along with the `hystfit`

library. For installing on Debian system, APT repository is available.

Once installed, create a fitting object `f`

from the `hystfit.Fit`

class, which is a subclass of the `unsatfit.Fit`

class. For more detailed information, refer to the unsatfit documentation.

```
import math
import numpy as np
import hystfit
f = hystfit.Fit()
```

The main drying curve can be modeled using either the van Genuchten (VG) or the Fredlund and Xing (FX) model. Once a model is chosen, parameters can be set in two ways:

(1) **Provide given parameters**:

- For the VG model, parameters (θ
_{s}, θ_{r}, α, n) can be set as follows:`f.set_vg(0.33, 0.05, 1 / 180, 1.65)`

- For the FX model, parameters (θ
_{s}, θ_{r}, a, m, n) can be set as follows:`f.set_fx(0.35, 0.02, 45, 1.25, 7.23)`

(2) **Optimize parameters from measured data (h, θ)**:

As `hystfit.Fit`

is a subclass of `unsatfit.Fit`

, all methods from `unsatfit`

are available in `f`

. You can obtain fitting parameters for the main drying curve using `unsatfit`

methods such as `get_init_vg`

and `get_wrf_vg`

. Note that for the VG function, parameter m should be converted to n=1/(1-m).

Set the water retention data (h, θ) for the main drying curve in `f.swrc`

as described in the `unsatfit`

documentation. After that, to obtain and set VG parameters:

```
qs, qr, a, m, q = f.get_wrf_vg()
n = 1 / (1 - m)
f.set_vg(qs, qr, a, n)
```

To obtain and set FX parameters:

```
f.set_fx(*f.get_wrf_fx())
```

If θ_{r} is equal to 0, you can use the `init_hyst`

method. For instance, for the VG model:

```
f.model_name = 'VG'
f.init_hyst()
print(f.message)
```

To use the FX model, simply set `f.model_name = 'FX'`

.

The hysteresis behavior in soil is determined by the advancing contact angle (γ_{A}) and the parameter b. These parameters are collectively set in a single tuple `p`

, defined as p=(cos γ_{A}, b). For instance, to set γ_{A} = 75° and b = 0.24:

```
p = (math.cos(math.radians(75)), 0.24)
```

Please note that in Zhou (2013), the contact angle is denoted by θ; however, we use γ instead to avoid confusion with volumetric water content. A detailed description of Zhou’s model is omitted here; please refer to Zhou’s paper for more information.

The receding contact angle γ_{r} is fixed at 0°, and its cosine value is initially set to `f.cos_gr = 1`

. This value remains unchanged unless manually altered by the user.

The change in h can be predicted from changes in θ using the `h`

method as shown below:

```
h = f.h(p, theta)
```

Where:

`p`

represents (cos γ_{A}, b).`theta`

is a numpy array representing the change in θ.`f.cos_g0`

is set to cos(γ_{0}).

By default, the initial contact angle γ_{0} is set at 0°, meaning `f.cos_g0 = 1`

. Since it is the same as the default γ_{R}, it indicates that the initial (h,θ) data point is on the main drying curve. After executing the `h`

method, γ_{0} is updated to the last state, enabling the reproduction of the hysteresis behavior through repeated usage of the `h`

method. γ_{0} is reset to 0° when `set_vg`

or `set_fx`

methods are called.

To set γ_{0} to a specific state of (h, θ):

```
f.cos_g0 = f.contact(h_ini, theta_ini)
```

Here, `h_ini`

and `theta_ini`

represent the initial values of h and θ, respectively.

To optimize hysteresis parameters based on changes in (h, θ)—for instance, from the main wetting curve—use the `opt`

method. Ensure that each (h, θ) dataset adheres to the data structure conventions of `unsatfit`

, and that the order of the data reflects the sequence of time events. The optimization can be performed with the following code:

```
f.opt(h, theta)
p = f.hyst
print(f'{f.message} R2 = {f.r2:.3}')
```

This code obtains the hysteresis parameter p=(cos γ_{A}, b) and outputs the result. In the optimization, the residual of ln(h) is used as the cost function. Note that the contact angle γ_{0} is set at the initial state before optimization and updated to its last state after this operation.