Instruction of hystfit

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()

Main drying curve

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'.

Hysteresis parameters

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.

Predicting change in h from change in θ

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(γ₀).

By default, the initial contact angle γ₀ 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, γ₀ is updated to the last state, enabling the reproduction of the hysteresis behavior through repeated usage of the h method. γ₀ is reset to 0° when set_vg or set_fx methods are called.

To set γ₀ 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.

Optimizing hysteresis parameters from changes in (h, θ)

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 γ₀ is set at the initial state before optimization and updated to its last state after this operation.