{
 "cells": [
  {
   "cell_type": "raw",
   "metadata": {
    "vscode": {
     "languageId": "raw"
    }
   },
   "source": [
    "---\n",
    "title: \"Pedotransfer functions\"\n",
    "author: \"Andres Patrignani\"\n",
    "source: \"pynotes agriscience\"\n",
    "date: \"1/19/2024\"\n",
    "keywords: [\"soil properties\", \"pedotransfer function\", \"Saxton\", \"Rawls\"]\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Pedotransfer functions (PTFs) are predictive models used in soil science to estimate soil properties that are usually hard to measure directly. They are particularly valuable for predicting soil hydraulic properties such as field capacity, permanent wilting point, and available water capacity from more easily measured soil characteristics like soil texture, bulk density, and organic matter content. The PTFs developed by Rawls et al. (1982) and Saxton & Rawls (2006) use soil texture and organic matter as inputs to predict key soil hydraulic parameters and have become widely used in hydrological modeling, irrigation management, and soil physics due to their simplicity and minimal data requirements.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import modules\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ptf(clay, sand, om=0.02, ec=0, rho=np.nan):\n",
    "    \"\"\" \n",
    "    Pedotransfer functions to determine soil physical properties from \n",
    "    soil texture and organic matter.\n",
    "\n",
    "    Parameters:\n",
    "    clay (interger, float): clay content\n",
    "    sand (integer, float): sand content\n",
    "    om (integer, float): Fraction of organic matter.\n",
    "\n",
    "    Returns:\n",
    "    Outputs: S is a dictionary containing scalar values.\n",
    "             V is a Pandas dataframe containing vector values across the range from 1500 to 0 kPa.\n",
    "\n",
    "    \"\"\" \n",
    "\n",
    "\n",
    "    # Section numbering does NOT follow the Saxton and Rawls paper\n",
    "    # 1. Permanent wilting point\n",
    "    # Calculation of soil water retention at -1500 kPa of tension (Table 1. Eq. 1)\n",
    "    theta_1500t = -0.024*sand + 0.487*clay + 0.006*om + 0.005*(sand*om) - 0.013*(clay*om) + 0.068*(sand*clay) + 0.031 \n",
    "    theta_1500 = theta_1500t + (0.14*theta_1500t - 0.02)\n",
    "\n",
    "    # 2. Field capacity\n",
    "    # Calculation of soil water retention at -33 kPa of tension (Table 1. Eq. 2)\n",
    "    theta_33t= -0.251*sand + 0.195*clay + 0.011*om + 0.006*(sand*om) - 0.027*(clay*om) + 0.452*(sand*clay) + 0.299\n",
    "    theta_33 = theta_33t + 1.283*(theta_33t)**2 - 0.374*(theta_33t) - 0.015\n",
    "\n",
    "    # 3. Saturation-Field capacity.\n",
    "    # Volumetric soil water content between 0 and -33 kPa of tension (Table 1. Eq. 3). \n",
    "    thetaS_33t= 0.278*sand + 0.034*clay + 0.022*om - 0.018*(sand*om) - 0.027*(clay*om) - 0.584*(sand*clay) + 0.078\n",
    "    thetaS_33= thetaS_33t + (0.636*thetaS_33t - 0.107)\n",
    "\n",
    "    # 4. Tension at air entry suction.\n",
    "    # Estimation of the tension at air entry (Table 1. Eq. 4).\n",
    "    psi_et = -21.67*sand - 27.93*clay - 81.97*thetaS_33 + 71.12*sand*thetaS_33 + 8.29*clay*thetaS_33 + 14.05*sand*clay + 27.16\n",
    "    psi_e = max(psi_et + 0.02*psi_et**2 - 0.113*psi_et - 0.7,0.1) # A minimum air-entry suction of 1 kPa was assumed to avoid negative values in sandy soils. \n",
    "\n",
    "    # 5. Saturation.\n",
    "    # Calculation of soil water retention at -0 kPa of tension (Table 1. Eq. 5)\n",
    "    thetaS = theta_33 + thetaS_33 - 0.097*sand + 0.043 # Saturation at normal density\n",
    "\n",
    "    # 6. Density effects (allow user to input of bulk density in case is known)\n",
    "    if ~np.isnan(rho):\n",
    "        rho_DF = rho\n",
    "    else:\n",
    "        rho_N = (1-thetaS)*2.65 # Particle density is assumed to be 2.65 g cm^-3. (Table 1. Eq. 6)\n",
    "        DF = 1 # Compaction degree.\n",
    "        rho_DF = rho_N*DF # (Table 1. Eq. 7).\n",
    "\n",
    "    thetaS_DF = 1 - (rho_DF/2.65) # (Table 1. Eq. 8).\n",
    "    theta_33_DF = theta_33 - 0.2*(thetaS-thetaS_DF) # (Table 1. Eq. 9).\n",
    "    thetaS_33_DF = thetaS_DF - theta_33_DF; # (Table 1. Eq. 10).\n",
    "    # thetaS_33_DF is used with equations that account for gravel, not present\n",
    "    # in this version. For more information see Table 1 Eqs. 19-22 in Saxton and Rawls, 2006.\n",
    "\n",
    "    # 7. Moisture tension\n",
    "    B = (np.log(1500) - np.log(33))/(np.log(theta_33) - np.log(theta_1500)) # (Table 1. Eq. 15).\n",
    "    A = np.exp(np.log(33) + B*np.log(theta_33)); # (Table 1. Eq. 14).\n",
    "\n",
    "    # Generate points of volumetric water content.\n",
    "    theta_piece1 = np.linspace(theta_1500, theta_33_DF, 1000) # theta from theta_1500 to theta_33\n",
    "    theta_piece2 = np.linspace(theta_33_DF, thetaS_DF, 1000) # theta from theta_33 to theta_0\n",
    "\n",
    "    # Generate soil matric potential values\n",
    "    psi_piece1 = A*theta_piece1**(-B) # Matric potential from 1500 to 33 kPa (Table 1. Eq. 11). \n",
    "    psi_piece2 = 33 - ((theta_piece2 - theta_33)*(33 - psi_e)/(thetaS_DF - theta_33)) # Matric potential from 33 kPa to air-entry suction (Table 1. Eq. 12).\n",
    "\n",
    "    # Merge water content and matric potential pieces to generate the soil water retention curve\n",
    "    theta = np.concatenate((theta_piece1, theta_piece2, [thetaS_DF]))\n",
    "    psi_m = np.concatenate((psi_piece1, psi_piece2, [0.1]))\n",
    "\n",
    "    # 8. Calculate volumetric water content at 10 kPa (perhaps a better definition of field capacity)\n",
    "    psi_diff10 = (psi_m - 10)**2;\n",
    "    theta_10_DF = theta[psi_diff10.argmin(axis=0)]\n",
    "\n",
    "    # 9. Moisture conductivity\n",
    "    par_lambda = 1/B; # (Table 1. Eq. 18).\n",
    "    Ks = 1930*(thetaS_DF - theta_33)**(3 - par_lambda) # (Table 1. Eq. 16) in [mm/h].\n",
    "    K_theta = Ks*(theta/thetaS_DF)**(3 + (2/par_lambda)) # (Table 1. Eq. 17).\n",
    "\n",
    "    # 10. Salinity effects\n",
    "    psi_o = 36*ec # (Table 1. Eq. 23).\n",
    "    psi_o_theta = thetaS/theta*psi_o # (Table 1. Eq. 24).\n",
    "\n",
    "    # 11. Set outputs\n",
    "    props = {'theta_sat': thetaS_DF, 'theta_10':theta_10_DF, 'theta_33':theta_33_DF, \n",
    "             'theta_1500':theta_1500, 'K_sat':Ks, 'psi_o':psi_o, 'rho':rho_DF,\n",
    "             'theta':theta, 'psi_m':psi_m, 'psi_o':psi_o_theta, 'K_theta':K_theta}\n",
    "\n",
    "    return props \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'theta_sat': 0.4339622641509434, 'theta_10': 0.43279196976870016, 'theta_33': 0.3504534007323006, 'theta_1500': 0.20392632400000002, 'K_sat': 1.5160753706941434, 'psi_o': array([5.85295059, 5.84874389, 5.84454323, ..., 2.75093194, 2.75040204,\n",
      "       2.75040204]), 'rho': 1.5, 'theta': array([0.20392632, 0.204073  , 0.20421967, ..., 0.43387867, 0.43396226,\n",
      "       0.43396226]), 'psi_m': array([1.50000000e+03, 1.49248108e+03, 1.48500521e+03, ...,\n",
      "       9.69247111e+00, 9.66865433e+00, 1.00000000e-01]), 'K_theta': array([4.09372980e-06, 4.14400991e-06, 4.19487078e-06, ...,\n",
      "       1.51112465e+00, 1.51607537e+00, 1.51607537e+00])}\n"
     ]
    }
   ],
   "source": [
    "# Example inputs\n",
    "clay = 0.34 # 34% clay\n",
    "sand = 0.22 # 22% sand\n",
    "om = 0.015  # 1.5% OM\n",
    "ec = 0.075 # dS/m\n",
    "rho_b = 1.5 # g/cm^3\n",
    "\n",
    "# Call function\n",
    "props = ptf(clay, sand, om, ec, rho_b)\n",
    "\n",
    "print(props)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make two plots next to each other  \n",
    "- The left one shows theta on x axis and psi_m on y axis. The right one shows K_theta on y axis and theta on x axis. \n",
    "- Put the curves from all three soil types in the same plots with different color and legend."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Your code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "Rawls, W.J., Brakensiek, D.L. and Saxton, K.E., 1982. Estimation of soil water properties. Transactions of the ASAE, 25(5), pp.1316-1320.\n",
    "           \n",
    "Saxton, K.E., Rawls, W., Romberger, J.S. and Papendick, R.I., 1986. Estimating generalized soil-water characteristics from texture 1. Soil science society of America Journal, 50(4), pp.1031-1036.\n",
    "\n",
    "Saxton, K.E. and Rawls, W.J., 2006. Soil water characteristic estimates by texture and organic matter for hydrologic solutions. Soil science society of America Journal, 70(5), pp.1569-1578."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}