{ "cells": [ { "cell_type": "markdown", "id": "92fe4842", "metadata": {}, "source": [ "# From k value to Gardner via HYPAGS + Curve Fitting Routine\n", "\n", "*Martin Vonk (2025)*\n", "\n", "## Overview\n", "\n", "This notebook demonstrates an integrated workflow that combines parameter estimation techniques:\n", "\n", "1. **HYPAGS estimation** - Get initial van Genuchten parameters from a single saturated conductivity measurement\n", "2. **Data generation** - Create water retention and conductivity data from the estimated model\n", "3. **Curve fitting** - Fit a different soil model (Gardner) to this data using least-squares optimization\n", "\n", "This workflow is useful for:\n", "- **Model conversion** - Converting between different soil model formulations (van Genuchten → Gardner, etc.)\n", "- **Hypothesis testing** - Comparing how well different models fit the same data\n", "- **Parameter refinement** - Using HYPAGS as an initial guess, then refining with additional data\n", "- **Workflow demonstration** - Showing how `pedon` components integrate together" ] }, { "cell_type": "markdown", "id": "9c52e27f", "metadata": {}, "source": [ "## Workflow Steps\n", "\n", "This notebook follows these steps:\n", "1. Import libraries and define helper plotting functions\n", "2. Estimate van Genuchten parameters using HYPAGS from a k value\n", "3. Generate synthetic water retention and conductivity data from the van Genuchten model\n", "4. Set up and visualize initial Gardner model parameters\n", "5. Fit the Gardner model to the synthetic data\n", "6. Compare the fitted Gardner model with the original van Genuchten curve\n", "\n", "## Setup" ] }, { "cell_type": "markdown", "id": "e8c2b6e8", "metadata": {}, "source": [ "### Imports" ] }, { "cell_type": "code", "execution_count": 1, "id": "6c5e347e", "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "import pedon as pe\n", "\n", "# pe.show_versions()" ] }, { "cell_type": "markdown", "id": "cff77c90", "metadata": {}, "source": [ "### Helper plot" ] }, { "cell_type": "code", "execution_count": 2, "id": "14ca8132", "metadata": {}, "outputs": [], "source": [ "def plot_compare(\n", " soilsample: pe.SoilSample, soilmodel: pe.SoilModel\n", ") -> np.typing.NDArray[plt.Axes]:\n", " f, ax = plt.subplots(1, 2, sharey=True, figsize=(7.0, 6.0))\n", " ax[0].scatter(soilsample.theta, soilsample.h, c=\"k\", s=10, label=\"Soil Sample\")\n", " _ = pe.soilmodel.plot_swrc(\n", " soilmodel, ax=ax[0], label=f\"Soil Model {soilmodel.__class__.__name__}\"\n", " )\n", " ax[0].set_yscale(\"log\")\n", " ax[0].set_xlim(0, 0.5)\n", " ax[0].set_yticks(soilsample.h)\n", " ax[0].set_xticks(np.linspace(0, 0.5, 6))\n", " ax[0].set_ylim(min(soilsample.h), max(soilsample.h))\n", " ax[1].scatter(soilsample.k, soilsample.h, c=\"k\", s=10)\n", " _ = pe.soilmodel.plot_hcf(soilmodel, ax=ax[1])\n", "\n", " ax[1].set_yscale(\"log\")\n", " ax[1].set_xscale(\"log\")\n", "\n", " k_left = 10 ** (np.floor(np.log10(min(soilsample.k))) - 1)\n", " k_right = 10 ** (np.ceil(np.log10(max(soilsample.k))) + 1)\n", " ax[1].set_xlim(k_left, k_right)\n", " ax[0].set_ylabel(r\"|$\\psi$| [cm]\")\n", " ax[0].set_xlabel(r\"$\\theta$ [-]\")\n", " ax[1].set_xlabel(r\"$K_s$ [cm/d]\")\n", " ncol = 3\n", " ax[0].legend(\n", " loc=(-0.02, 1),\n", " fontsize=6,\n", " frameon=False,\n", " ncol=ncol,\n", " columnspacing=0.8,\n", " handlelength=2.5,\n", " )\n", "\n", " f.align_xlabels()\n", " return ax" ] }, { "cell_type": "markdown", "id": "364af67c", "metadata": {}, "source": [ "## Step 1: Estimate van Genuchten Parameters via HYPAGS\n", "\n", "We start with a single measurement: saturated hydraulic conductivity of 100 cm/d (typical for sandy soil). HYPAGS will estimate the complete van Genuchten parameter set from this single value.\n", "\n", "Note: HYPAGS expects k in m/s, so we convert from cm/d." ] }, { "cell_type": "code", "execution_count": 3, "id": "24f8d4a4", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Genuchten(k_s=100.0, theta_r=0.108, theta_s=np.float64(0.260064701261806), alpha=np.float64(3.2905002387762474), n=np.float64(1.7228201552429758), l=0.5)" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "k_cmd = 100.0 # cm/d\n", "k_ms = k_cmd / 86400 / 100 # m/s since HYPAGS expects that\n", "genuchten = pe.SoilSample(k=np.array([k_ms])).hypags()\n", "genuchten.k_s = k_cmd # convert back to cm/d\n", "genuchten" ] }, { "cell_type": "markdown", "id": "f8d6a59a", "metadata": {}, "source": [ "## Step 2: Convert to Gardner Model via Curve Fitting\n", "\n", "Now we'll fit a Gardner model to the van Genuchten data. This demonstrates model conversion: taking one model's predictions and fitting a different model to match them.\n", "\n", "### Why This Workflow?\n", "\n", "Sometimes you need to convert between models because:\n", "- Your simulation software only accepts a specific model\n", "- You want to compare model performance on the same data\n", "- You're testing whether a simpler model (like Gardner) works as well as a more complex one (like van Genuchten)\n", "\n", "### Sample Pressure Heads\n", "\n", "We create synthetic data by evaluating the van Genuchten model at various pressure head values. Note the careful range selection to avoid numerical issues where conductivity becomes extremely small." ] }, { "cell_type": "markdown", "id": "a794e0c7", "metadata": {}, "source": [ "### Sample Genuchten curve for data points" ] }, { "cell_type": "code", "execution_count": 4, "id": "79faf394", "metadata": {}, "outputs": [], "source": [ "h = np.logspace(\n", " -3, 2, 100\n", ") # carefull with sampling h since k goes to zero very fast for large h which gives trouble with logs\n", "k = genuchten.k(h)\n", "theta = genuchten.theta(h)\n", "soilsample = pe.SoilSample(h=h, k=k, theta=theta)" ] }, { "cell_type": "markdown", "id": "e4ebeda5", "metadata": {}, "source": [ "### Setting Parameter Bounds\n", "\n", "Before fitting, we define reasonable bounds for each Gardner parameter. The bounds have three values:\n", "- **p_ini**: Initial guess for the optimizer\n", "- **p_min**: Minimum allowed value\n", "- **p_max**: Maximum allowed value\n", "\n", "Good bounds are crucial for successful fitting - they guide the optimization and prevent unrealistic values." ] }, { "cell_type": "code", "execution_count": 5, "id": "7006da85", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "WARNING:root:No default parameter bounds for SoilModel type Gardner\n" ] }, { "data": { "text/html": [ "
| \n", " | p_ini | \n", "p_min | \n", "p_max | \n", "
|---|---|---|---|
| k_s | \n", "100.000000 | \n", "0.0000 | \n", "110.0 | \n", "
| theta_s | \n", "0.260065 | \n", "0.2000 | \n", "0.3 | \n", "
| m | \n", "0.300000 | \n", "0.0001 | \n", "0.5 | \n", "
| c | \n", "0.150000 | \n", "0.0001 | \n", "0.5 | \n", "