Sample code of unsatfit

Easiest way to start learning how to use unsatfit is to run sample codes as instructed in this page. You can optimize parameters of WRF (water retention function) and HCF (hydraulic conductivity function) of various available models to measured data set. You can test with sample data provided in this page. For optimizing only WRF parameters, please use SWRC Fit.

List of sample codes

• Models
  • VGM (van Genuchten - Mualem) model with full instruction
  • Bimodal model with general HCF
    • KBC (KO₁BC₂-CH) model
    • DVC (dual-VG-CH) model
    • DBC (dual-BC-CH) model
  • Peters model
• Multiple curves
• Contour plot

VGM (van Genuchten - Mualem) model

\[\begin{cases} \theta(h) &=& (\theta_s - \theta_r)S(h) + \theta_r\\ K(h) &=& K_s {S(h)}^p \bigl[1-(1-S(h)^{1/m})^m\bigr]^2\\ \end{cases}\]

where

\[\begin{eqnarray} S(h) &=& \biggl[\dfrac{1}{1+(ah)^n}\biggr]^m\\ m &=& 1-1/n \end{eqnarray}\]

• Install Python 3 and unsatfit. Install pandas by python -m pip install pandas.
• Prepare datafile as csv (comma-separated values) format for (h, θ) data in the filename “swrc.csv” with a header of “h, theta”. See sample data of clay (Unsoda 2362).
• Prepare (h, K) data in the filename “hcc.csv” with a header of “h, K”. See sample (h, K) data of clay (Unsoda 2362).
• Download a sample code for VGM model.
• Read the sample code and edit initial parameters, bound of parameters and figure setting as written in the comment of the sample code if necessary. For more figure options, see source code of unsatfit.
• Run the sample code at the same directory with data files (swrc.csv and hcc.csv). For running the code on Mac or unix-like system, edit the first line (shebang) and mark the file executable by chmod +x VGM.py. For running on Windows, please refer to Python on Windows FAQ.
• It first optimizes WRF parameters (θ_s, θ_r, a, n) of VG model, and then optimizes HCF parameters (K_s, p) of VG Mualem model, or modified VG model (h_b=2) Vogel et al. (2000) when n > 1.1. h_b value can be changed as HB value in the sample code.
• Fitted parameters are shown at the standard output, where qs and qr means θ_s and θ_r respectively, and R2 q means R² for θ of water retention curve and R2 logK means R² for log(K) of hydraulic conductivity curve.
• Note that the program is unit independent, meaning that the unit of the parameters depends on the unit of the input data. Unit of pressure head is assumed as cm for a (cm^-1) and h_b.
• Figure file is produced as VG.png. For use in papers, pdf file can be produced as instructed in the sample code.

Result with sample data of clay (Unsoda 2362) is shown below.

Result with sample data of Gilat loam is shown below. In this case, bimodal model is appropriate, as shown below.

Bimodal model with general HCF

Bimodal model with general HCF (Seki et al., 2022) can represent water retention and hydraulic conductivity of various types of soil in a wide range of pressure head, as verified in Seki et al., 2023. You can conduct the same fitting as written in the paper by using the following sample codes.

KBC (KO₁BC₂-CH) model (θ_r=0, r=1)

\[\begin{cases} \theta(h) &=& \theta_s S(h)\\ K(h) &=& K_s {S(h)}^p \gamma^{-1} \Biggl[ b_1 Q \biggl[\dfrac{\ln(h/H)}{\sigma_1} + q\sigma_1\biggr]+ b_2 (h/H)^{-\lambda_2 - q} \Biggr] \end{cases}\]

where

\[\begin{eqnarray} S(h) &=& \begin{cases}w S_1(h) + (1-w)\left(h/H\right)^{-\lambda_2} & (h>H)\\ w S_1(h) + 1-w & (h \le H)\end{cases}\\ S_1(h) &=& Q \biggl[\dfrac{\ln(h/H)}{\sigma_1}\biggr]\\ Q(x) &=& \frac{1-\mathrm{erf}(x/\sqrt{2})}{2}\\ b_1 &=& w \exp\biggl(\frac{q^2 \sigma_1^2}{2}\biggr)\\ b_2 &=& (1-w)\biggl(\frac{q}{\lambda_2} + 1\biggr)^{-1}\\ \gamma &=& b_1+b_2 \end{eqnarray}\]

• Use sample code for KBC model. See instruction in the VGM model above.
• It first optimizes WRF parameters (θ_s, w, H, σ₁, λ₂) of KBC model (θ_r=0), and then optimizes general HCF parameters (K_s, p, q) of KBC model (r=1), or modified KBC model (h_b=2, r=1) when σ₁ > 2.
• Result with sample data of Gilat loam is shown below.

DVC (dual-VG-CH) model (θ_r=0, q=1)

\[\begin{cases} \theta(h) &=& \theta_s S(h)\\ K(h) &=& K_s {S(h)}^p \bigl[w\Gamma_1(h) + (1-w)\Gamma_2(h)\bigr]^r\\ \end{cases}\]

where

\[\begin{eqnarray} S(h) &=& w\bigl[1+(ah)^{n_1}\bigr]^{-m_1} + (1-w)\bigl[1+(ah)^{n_2}\bigr]^{-m_2}\\ m_i&=&1-1/{n_i}\\ \Gamma_i(h) &=& 1-\biggl[1-\big[1+(ah)^{n_i}\bigr]^{-1}\biggr]^{m_i} \end{eqnarray}\]

• Use sample code for DVC model. See instruction in the VGM model above.
• It optimizes WRF parameters (θ_s, w, a, n₁, n₂) , and then optimizes general HCF parameters (K_s, p, r) of DVC model or modified DVC model (h_b=2) when n₁ or n₂ is smaller than 1.1.
• Result with sample data of Gilat loam is shown below.

DBC (dual-BC-CH) model (θ_r=0, r=1)

\[\begin{eqnarray} \theta(h) &=& \theta_s S(h)\\ K(h) &=& \begin{cases}K_s {S(h)}^p \gamma^{-1} \bigl[ wB_1 \Gamma_1(h) + (1-w)B_2 \Gamma_2(h) \bigr] & (h>H)\\ K_s & (h \le H)\end{cases}\\ \end{eqnarray}\]

where

\[\begin{eqnarray} S(h) &=& \begin{cases}w \left(h / H\right)^{-\lambda_1} + (1-w)\left(h / H\right)^{-\lambda_2} & (h>H)\\ 1 & (h \le H)\end{cases}\\ B_i &=& \biggl(\frac{q}{\lambda_i} + 1\biggr)^{-1}\\ \gamma &=& wB_1+(1-w)B_2\\ \Gamma_i(h) &=& \left(h / H\right)^{-\lambda_i-q} \end{eqnarray}\]

• Use sample code for DBC model. See instruction in the VGM model above.
• It optimizes WRF parameters (θ_s, w, λ₁, λ₂), and then optimizes general HCF parameters (K_s, p, q) of DBC model.
• Result with sample data of Gilat loam is shown below.

Peters model (θ_r=0)

\[\begin{cases} \theta(h) &=& \theta_s S(h)\\ K(h) &=& K_s \bigl[ (1-\omega)K_1(h)+\omega K_2(h) \bigr]\\ \end{cases}\]

where

\[\begin{eqnarray} S(h) &=& w S_1(h) + (1-w)S_2(h)\\ S_1(h) &=& Q \biggl[\dfrac{\ln(h/H)}{\sigma}\biggr]\\ Q(x) &=& \frac{1-\mathrm{erf}(x/\sqrt{2})}{2}\\ S_2(h) &=& \begin{cases}\dfrac{L(h_0)-L(h)}{L(h_0)-L(H)} & (h>H)\\ 1 & (h \le H)\end{cases}\\ L(h) &=& \ln(1+h/H)\\ K_1(h) &=& {S_1(h)}^p \Biggl[ Q \biggl[\dfrac{\ln(h/H)}{\sigma} + \sigma \biggr] \Biggr]^2\\ K_2(h) &=& \begin{cases}(h/H)^a & (h>H)\\ 1 & (h \le H)\end{cases}\\ \end{eqnarray}\]

• Use sample code for Peters model. See instruction in the VGM model above.
• h₀=6.3×10⁶ is constant and can be edited as H0 value in the sample code.
• It optimizes WRF parameters (θ_s, w, H, σ), and then optimizes general HCF parameters (K_s, p, a, ω) of Peters model or modified Peters model (h_b=2) when σ > 2.
• Result with sample data of Gilat loam is shown below.

Multiple curves

• Use sample code for multiple curves. See instruction in the VGM model above.
• Draw DBC, DVC, KBC and Peters model in the same figure.
• Result with sample data of Gilat loam is shown below.

Contour plot

• Use sample code for contour plot. See instruction in the VGM model above.
• It draws a contour plot of (p, q) for RMSE of estimated log₁₀(K) for KBC (r=1) optimization.
• Result with sample data of Gilat loam is shown below.