{ "cells": [ { "cell_type": "raw", "metadata": { "vscode": { "languageId": "raw" } }, "source": [ "---\n", "title: \"Soil water retention curves\"\n", "author: \"Andres Patrignani\"\n", "source: Pynotes in Agriscience\n", "date: \"1/19/2024\"\n", "keywords: [\"soil water retention\", \"matric potential\", \"soil moisture\"]\n", "---" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "A soil water retention curve represents the relationship between the amount of water in the soil typically expressed on a volume basis (i.e., volumetric soil water content) and the energy-state of the soil water expressed in units of energy per unit volume ($J/m^3$), energy per unit mass ($J/kg$), or units of pressure/tension ($kPa$). This curve is fundamental for understanding soil water storage and flow in porous media in applications like irrigation management, drainage, ground water recharge, and soil water availability to plants. \n", "\n", "The shape of the curve is largely dictated by soil physical properties such as texture, structure, and organic matter that control the surface area and pore size distribution of the soil. Soil water retention curves are usually determined empirically by collecting soil samples and using a series of laboratory instruments to quantify the soil water content at different tension levels. In this exercise we will fit common models used in the scientific literature to a dataset determined in our soil physics lab.\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "# Import modules\n", "import pandas as pd\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.optimize import curve_fit\n", "from scipy.special import erfc\n" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
matrictheta
00.0010.448
10.0010.448
20.0560.448
30.1730.447
40.2860.447
\n", "
" ], "text/plain": [ " matric theta\n", "0 0.001 0.448\n", "1 0.001 0.448\n", "2 0.056 0.448\n", "3 0.173 0.447\n", "4 0.286 0.447" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# Read dataset\n", "df = pd.read_csv('soil_water_retention_curve.csv')\n", "df['theta'] = df['theta']/100\n", "df.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Define soil water retetnion models" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Brooks and Corey model (1964)\n", "\n", "$$ \\frac{\\theta - \\theta_r}{\\theta_s - \\theta_r} = \\Bigg( \\frac{\\psi_e}{\\psi} \\Bigg)^\\lambda $$\n", "\n", "where:\n", "\n", "$\\theta$ is the volumetric water content $(cm^3/cm^3)$
\n", "$\\theta_r$ is the residual water content $(cm^3/cm^3)$
\n", "$\\theta_s$ is the saturation water content $(cm^3/cm^3)$
\n", "$\\psi$ is the matric potential $(kPa)$
\n", "$\\psi_e$ is the air-entry suction $(kPa)$
\n", "$\\lambda$ is a parameter related to the pore-size distribution
\n" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def brooks_corey_model(x,alpha,lmd,psi_e,theta_r,theta_s):\n", " \"\"\"\n", " Function that computes volumetric water content from soil matric potential\n", " using the Brooks and Corey (1964) model.\n", " \"\"\"\n", " theta = np.minimum(theta_r + (theta_s-theta_r)*(psi_e/x)**(lmd), theta_s)\n", " return theta" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### van Genuchten model (1980)\n", "\n", "$$ \\frac{\\theta - \\theta_r}{\\theta_s - \\theta_r} = [1 + (-\\alpha \\psi)^n]^{-m} $$\n", " \n", "where:\n", "\n", "$\\alpha$ is a fitting parameter inversely related to $\\psi_e$ $(kPa^{-1})$
\n", "$n$ is a fitting parameter
\n", "$m$ is a fitting parameter that is often assumed to be $M=1-1/n$, but in this example it was left as a free parameter (see article by Groenevelt and Grant (2004) for mo details.
\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def van_genuchten_model(x,alpha,n,m,theta_r,theta_s):\n", " \"\"\"\n", " Function that computes volumetric water content from soil matric potential\n", " using the van Genuchten (1980) model.\n", " \"\"\"\n", " theta = theta_r + (theta_s-theta_r)*(1+(alpha*x)**n)**-(1-1/n)\n", " return theta" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Kosugi model (1994)\n", "\n", "$$ \\theta = \\theta_r + \\frac{1}{2} (\\theta_s - \\theta_r) \\ erfc \\Bigg( \\frac{ln(\\psi/ \\psi_m)}{\\sigma \\sqrt(2)} \\Bigg)$$\n", "\n", "where:\n", "\n", "$erfc$ is the complementary error function
\n", "$\\sigma$ is the standard deviation of $ln(\\psi_{med})$
\n", "$\\psi_m$ is the median matric potential\n", "\n" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def kosugi_model(x,hm,sigma,theta_r,theta_s):\n", " \"\"\"\n", " Function that computes volumetric water content from soil matric potential\n", " using the Kosugi (1994) model. Function was implemented accoridng \n", " to Eq. 3 in Pollacco et al., 2017.\n", " \"\"\"\n", " theta = theta_r + 1/2*(theta_s-theta_r)*erfc(np.log(x/hm)/(sigma*np.sqrt(2))) \n", " return theta" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Groenvelt-Grant model (2004)\n", "$$ \\theta = k_1 \\Bigg[ exp\\Bigg( \\frac{-k_0}{\\big(10^{5.89}\\big)^n} \\Bigg) - exp\\Bigg( \\frac{-k_0}{\\psi ^n} \\Bigg) \\Bigg]$$\n", "\n", "where:\n", "\n", "$k_0$, $k_1$, and $n$ are fitting parameters. The value $10^{5.89}$ is the matric potential in $kPa$ at oven-dry conditions (105 degrees Celsius)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def groenevelt_grant_model(x,k0,k1,n,theta_s):\n", " \"\"\"\n", " Function that computes volumetric water content from soil matric potential \n", " using the Groenevelt-Grant (2004) model.\n", " \"\"\"\n", " theta = k1 * ( np.exp(-k0/ (10**5.89)**n) - np.exp(-k0/(x**n)) ) # Eq. 5 in Groenevelt and Grant, 2004\n", " return theta\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Define error models" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# Error models\n", "mae_fn = lambda x,y: np.round(np.mean(np.abs(x-y)),3)\n", "rmse_fn = lambda x,y: np.round(np.sqrt(np.mean((x-y)**2)),3)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Fit soil water retention models to dataset" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# Define variables\n", "x_obs = df[\"matric\"]\n", "y_obs = df[\"theta\"]\n", "x_curve = np.logspace(-1.5,5,1000)\n", "\n", "# Empty list to collect output of each model\n", "output = []\n", "\n", "# van Genuchten model\n", "p0 = [0.02,1.5,1,0.1,0.5]\n", "bounds = ([0.001,1,0,0,0.3], [1,10,25,0.3,0.6])\n", "par_opt, par_cov = curve_fit(van_genuchten_model, x_obs, y_obs, p0=p0, bounds=bounds)\n", "y_curve = van_genuchten_model(x_curve, *par_opt)\n", "output.append({'name':'van Genuchten',\n", " 'y_curve':y_curve,\n", " 'mae':mae_fn(van_genuchten_model(x_obs, *par_opt), y_obs),\n", " 'rmse':rmse_fn(van_genuchten_model(x_obs, *par_opt), y_obs),\n", " 'par_values':par_opt,\n", " 'par_names':van_genuchten_model.__code__.co_varnames,\n", " 'color':'black'})\n", " \n", "# Brooks and Corey\n", "p0=[0.02, 1, 10, 0.1, 0.5]\n", "bounds=([0.001, 0.1, 1, 0, 0.3], [1, 10, 100, 0.3, 0.6])\n", "par_opt, par_cov = curve_fit(brooks_corey_model, x_obs, y_obs, p0=p0, bounds=bounds)\n", "y_curve = brooks_corey_model(x_curve, *par_opt)\n", "output.append({'name':'Brooks and Corey',\n", " 'y_curve':y_curve,\n", " 'mae':mae_fn(brooks_corey_model(x_obs, *par_opt), y_obs),\n", " 'rmse':rmse_fn(brooks_corey_model(x_obs, *par_opt), y_obs),\n", " 'par_values':par_opt,\n", " 'par_names':brooks_corey_model.__code__.co_varnames,\n", " 'color':'tomato'})\n", "\n", "\n", "# Kosugi\n", "p0=[50, 1, 0.1, 0.5]\n", "bounds=([1, 1, 0, 0.3], [500, 10, 0.3, 0.6])\n", "par_opt, par_cov = curve_fit(kosugi_model, x_obs, y_obs, p0=p0, bounds=bounds)\n", "y_curve = kosugi_model(x_curve, *par_opt)\n", "output.append({'name':'Kosugi',\n", " 'y_curve':y_curve,\n", " 'mae':mae_fn(kosugi_model(x_obs, *par_opt), y_obs),\n", " 'rmse':rmse_fn(kosugi_model(x_obs, *par_opt), y_obs),\n", " 'par_values':par_opt,\n", " 'par_names':kosugi_model.__code__.co_varnames,\n", " 'color':'darkgreen'})\n", "\n", "\n", "# Groenevelt-Grant\n", "p0=[5, 1, 2, 0.5]\n", "bounds=([1, 0.1, 0.1, 0.3], [2000, 10, 5, 0.6])\n", "par_opt, par_cov = curve_fit(groenevelt_grant_model, x_obs, y_obs, p0=p0, bounds=bounds)\n", "y_curve = groenevelt_grant_model(x_curve, *par_opt)\n", "output.append({'name':'Groenevelt-Grant',\n", " 'y_curve':y_curve,\n", " 'mae':mae_fn(groenevelt_grant_model(x_obs, *par_opt), y_obs),\n", " 'rmse':rmse_fn(groenevelt_grant_model(x_obs, *par_opt), y_obs),\n", " 'par_values':par_opt,\n", " 'par_names':groenevelt_grant_model.__code__.co_varnames,\n", " 'color':'lightblue'})\n" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "van Genuchten\n", "alpha = 0.03982152189968467\n", "n = 1.4290698696422959\n", "m = 3.916110721036148e-06\n", "theta_r = 1.9220461763010268e-12\n", "theta_s = 0.4376342588946559\n", "\n", "Brooks and Corey\n", "alpha = 0.02\n", "lmd = 0.2559256472395729\n", "psi_e = 10.068395061444875\n", "theta_r = 2.186268889571107e-21\n", "theta_s = 0.43043106156880795\n", "\n", "Kosugi\n", "hm = 122.18916557122672\n", "sigma = 1.9621133806177067\n", "theta_r = 5.009929658632887e-22\n", "theta_s = 0.44660152443876516\n", "\n", "Groenevelt-Grant\n", "k0 = 8.104570461072237\n", "k1 = 0.44363145314935387\n", "n = 0.5029795445930895\n", "theta_s = 0.5\n", "\n" ] } ], "source": [ "# Print parameters\n", "for k,model in enumerate(output):\n", " print(model['name'])\n", " for par_name, par_value in list(zip(model['par_names'][1:], model['par_values'])):\n", " print(par_name, '=', par_value)\n", " print('')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Make a figure showing volumetric water content on the y axis and soil matric potential on the x axis. \n", "- It includes the observations (data from the csv file) as dots and the four fitted curves with different colors\n", "- Axis lables, units, and legends \n", "- The x axis is on log scale" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "# Your code here" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## References\n", "\n", "Brooks, R. H. (1965). Hydraulic properties of porous media. Colorado State University.\n", "\n", "Groenevelt, P.H. and Grant, C.D., 2004. A new model for the soil‐water retention curve that solves the problem of residual water contents. European Journal of Soil Science, 55(3), pp.479-485.\n", "\n", "Kosugi, K. I. (1994). Three‐parameter lognormal distribution model for soil water retention. Water Resources Research, 30(4), 891-901.\n", "\n", "Pollacco, J. A. P., Webb, T., McNeill, S., Hu, W., Carrick, S., Hewitt, A., & Lilburne, L. (2017). Saturated hydraulic conductivity model computed from bimodal water retention curves for a range of New Zealand soils. Hydrology and Earth System Sciences, 21(6), 2725-2737.\n", "\n", "van Genuchten, M.T., 1980. A closed form equation for predicting hydraulic conductivity of unsaturated soils: Journal of the Soil Science Society of America." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.13" } }, "nbformat": 4, "nbformat_minor": 4 }