Richards Equation Transient.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simple 1D Richards' Equation Model\n",
    "\n",
    "This is a 1D Richards' equation model, written by Andrew Ireson, 3 November 2015.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import all of the basic libraries (you will always need these)\n",
    "%matplotlib inline\n",
    "from matplotlib import pyplot as pl\n",
    "import numpy as np\n",
    "\n",
    "# Import a library that contains soil moisture properties and functions\n",
    "import vanGenuchten as vg\n",
    "\n",
    "# Import ODE solvers\n",
    "from scipy.interpolate import interp1d\n",
    "from scipy.integrate import odeint\n",
    "\n",
    "# Select which soil properties to use\n",
    "p=vg.HygieneSandstone()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Richards equation solver\n",
    "# This is a function that calculated the right hand side of Richards' equation. You\n",
    "# will not need to modify this function, unless you are doing something advanced. \n",
    "# This block of code must be executed so that the function can be later called.\n",
    "\n",
    "def RichardsModel(psi,t,dz,n,p,vg,qTop,qBot,psiTop,psiBot):\n",
    "       \n",
    "    # Basic properties:\n",
    "    C=vg.CFun(psi,p)\n",
    "   \n",
    "    # initialize vectors:\n",
    "    q=np.zeros(n+1)\n",
    "    \n",
    "    # Upper boundary\n",
    "    if qTop == []:\n",
    "        KTop=vg.KFun(np.zeros(1)+psiTop,p)\n",
    "        q[n]=-KTop*((psiTop-psi[n-1])/dz*2+1)\n",
    "    else:\n",
    "        q[n]=qTop\n",
    "    \n",
    "    # Lower boundary\n",
    "    if qBot == []:\n",
    "        if psiBot == []:\n",
    "            # Free drainage\n",
    "            KBot=vg.KFun(np.zeros(1)+psi[0],p)\n",
    "            q[0]=-KBot\n",
    "        else:\n",
    "            # Type 1 boundary\n",
    "            KBot=vg.KFun(np.zeros(1)+psiBot,p)\n",
    "            q[0]=-KBot*((psi[0]-psiBot)/dz*2+1.0)    \n",
    "    else:\n",
    "        # Type 2 boundary\n",
    "        q[0]=qBot\n",
    "    \n",
    "    # Internal nodes\n",
    "    i=np.arange(0,n-1)\n",
    "    Knodes=vg.KFun(psi,p)\n",
    "    Kmid=(Knodes[i+1]+Knodes[i])/2.0\n",
    "    \n",
    "    j=np.arange(1,n)\n",
    "    q[j]=-Kmid*((psi[i+1]-psi[i])/dz+1.0)\n",
    "    \n",
    "    \n",
    "    \n",
    "    # Continuity\n",
    "    i=np.arange(0,n)\n",
    "    dpsidt=(-(q[i+1]-q[i])/dz)/C\n",
    "    \n",
    "    return dpsidt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Richards equation solver\n",
    "# This is a function that calculated the right hand side of Richards' equation. You\n",
    "# will not need to modify this function, unless you are doing something advanced. \n",
    "# This block of code must be executed so that the function can be later called.\n",
    "\n",
    "def RichardsModelTransient(psi,t,dz,n,p,vg,qTfun,qBot,psiTop,psiBot):\n",
    "       \n",
    "    # Basic properties:\n",
    "    C=vg.CFun(psi,p)\n",
    "   \n",
    "    # initialize vectors:\n",
    "    q=np.zeros(n+1)\n",
    "    \n",
    "    # Upper boundary\n",
    "#     if qTop == []:\n",
    "#         KTop=vg.KFun(np.zeros(1)+psiTop,p)\n",
    "#         q[n]=-KTop*((psiTop-psi[n-1])/dz*2+1)\n",
    "#     else:\n",
    "#         q[n]=qTop\n",
    "\n",
    "\n",
    "    if t>100:\n",
    "        q[n]=qTfun(100)\n",
    "    else:\n",
    "        q[n]=qTfun(t)\n",
    "    \n",
    "    # Lower boundary\n",
    "    if qBot == []:\n",
    "        if psiBot == []:\n",
    "            # Free drainage\n",
    "            KBot=vg.KFun(np.zeros(1)+psi[0],p)\n",
    "            q[0]=-KBot\n",
    "        else:\n",
    "            # Type 1 boundary\n",
    "            KBot=vg.KFun(np.zeros(1)+psiBot,p)\n",
    "            q[0]=-KBot*((psi[0]-psiBot)/dz*2+1.0)    \n",
    "    else:\n",
    "        # Type 2 boundary\n",
    "        q[0]=qBot\n",
    "    \n",
    "    # Internal nodes\n",
    "    i=np.arange(0,n-1)\n",
    "    Knodes=vg.KFun(psi,p)\n",
    "    Kmid=(Knodes[i+1]+Knodes[i])/2.0\n",
    "    \n",
    "    j=np.arange(1,n)\n",
    "    q[j]=-Kmid*((psi[i+1]-psi[i])/dz+1.0)\n",
    "    \n",
    "    \n",
    "    \n",
    "    # Continuity\n",
    "    i=np.arange(0,n)\n",
    "    dpsidt=(-(q[i+1]-q[i])/dz)/C\n",
    "    \n",
    "    return dpsidt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "psi = np.linspace(-10,0)\n",
    "theta = vg.thetaFun(psi,p)\n",
    "C=vg.CFun(psi,p)\n",
    "K=vg.KFun(psi,p)\n",
    "\n",
    "pl.rcParams['figure.figsize'] = (5.0, 10.0)\n",
    "pl.subplot(311)\n",
    "pl.plot(psi,theta)\n",
    "pl.ylabel(r'$\\theta$', fontsize=20)\n",
    "pl.subplot(312)\n",
    "pl.plot(psi,C)\n",
    "pl.ylabel(r'$C$',fontsize=20)\n",
    "pl.subplot(313)\n",
    "pl.plot(psi,K)\n",
    "pl.ylabel(r'$K$', fontsize=20)\n",
    "pl.xlabel(r'$\\psi$', fontsize=20)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This block of code sets up and runs the model\n",
    "\n",
    "# Boundary conditions\n",
    "qTop=-0.01\n",
    "qBot=[]\n",
    "psiTop=[]\n",
    "psiBot=[]\n",
    "\n",
    "# Grid in space\n",
    "dz=0.1\n",
    "ProfileDepth=5\n",
    "z=np.arange(dz/2.0,ProfileDepth,dz)\n",
    "n=z.size\n",
    "\n",
    "# Grid in time\n",
    "t = np.linspace(0,10,101)\n",
    "\n",
    "# Initial conditions\n",
    "psi0=-z\n",
    "\n",
    "# Solve\n",
    "psi=odeint(RichardsModel,psi0,t,args=(dz,n,p,vg,qTop,qBot,psiTop,psiBot),mxstep=5000000);\n",
    "\n",
    "print \"Model run successfully\"    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Post process model output to get useful information\n",
    "\n",
    "# Get water content\n",
    "theta=vg.thetaFun(psi,p)\n",
    "\n",
    "# Get total profile storage\n",
    "S=theta.sum(axis=1)*dz\n",
    "\n",
    "# Get change in storage [dVol]\n",
    "dS=np.zeros(S.size)\n",
    "dS[1:]=np.diff(S)/(t[1]-t[0])\n",
    "\n",
    "# Get infiltration flux\n",
    "if qTop == []:\n",
    "    KTop=vg.KFun(np.zeros(1)+psiTop,p)\n",
    "    qI=-KTop*((psiTop-psi[:,n-1])/dz*2+1)\n",
    "else:\n",
    "    qI=np.zeros(t.size)+qTop\n",
    "    \n",
    "# Get discharge flux\n",
    "if qBot == []:\n",
    "    if psiBot == []:\n",
    "        # Free drainage\n",
    "        KBot=vg.KFun(psi[:,0],p)\n",
    "        qD=-KBot\n",
    "    else:\n",
    "        # Type 1 boundary\n",
    "        KBot=vg.KFun(np.zeros(1)+psiBot,p)\n",
    "        qD=-KBot*((psi[:,0]-psiBot)/dz*2+1.0)\n",
    "else:\n",
    "    qD=np.zeros(t.size)+qBot\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot vertical profiles\n",
    "pl.rcParams['figure.figsize'] = (10.0, 10.0)\n",
    "for i in range(0,t.size-1):\n",
    "    pl.subplot(121)\n",
    "    pl.plot(psi[i,:],z)\n",
    "    pl.subplot(122)\n",
    "    pl.plot(theta[i,:],z)\n",
    "\n",
    "pl.subplot(121)\n",
    "pl.ylabel('Elevation [m]',fontsize=20)\n",
    "pl.xlabel(r'$\\psi$ [m]',fontsize=20)\n",
    "pl.subplot(122)\n",
    "pl.xlabel(r'$\\theta$ [-]',fontsize=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot timeseries\n",
    "dt = t[2]-t[1]\n",
    "pl.plot(t,dS,label='Change in storage')\n",
    "pl.hold(True)\n",
    "pl.plot(t,-qI,label='Infiltration')\n",
    "pl.plot(t,-qD,label='Discharge')\n",
    "pl.legend()\n",
    "pl.ylim((0,0.02))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t=np.arange(0,101,1)\n",
    "qT=np.zeros(len(t))-0.01\n",
    "print(qT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This block of code sets up and runs the model\n",
    "\n",
    "# Boundary conditions\n",
    "qTop=-0.01\n",
    "qBot=[]\n",
    "psiTop=[]\n",
    "psiBot=[]\n",
    "\n",
    "# Grid in space\n",
    "dz=0.1\n",
    "ProfileDepth=5\n",
    "z=np.arange(dz/2.0,ProfileDepth,dz)\n",
    "n=z.size\n",
    "\n",
    "# Time grid and transient infiltration flux\n",
    "t=np.arange(0,101,1)\n",
    "qT=np.zeros(len(t))\n",
    "qT[10:20]=-0.01\n",
    "qTfun=interp1d(t,qT)\n",
    "\n",
    "# Initial conditions\n",
    "psi0=-z\n",
    "\n",
    "# Solve\n",
    "psi=odeint(RichardsModelTransient,psi0,t,args=(dz,n,p,vg,qTfun,qBot,psiTop,psiBot),mxstep=5000000);\n",
    "\n",
    "print \"Model run successfully\"   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pl.plot(t,psi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot vertical profiles\n",
    "pl.rcParams['figure.figsize'] = (10.0, 10.0)\n",
    "for i in range(0,t.size-1):\n",
    "    pl.subplot(121)\n",
    "    pl.plot(psi[i,:],z)\n",
    "    pl.subplot(122)\n",
    "    pl.plot(theta[i,:],z)\n",
    "\n",
    "pl.subplot(121)\n",
    "pl.ylabel('Elevation [m]',fontsize=20)\n",
    "pl.xlabel(r'$\\psi$ [m]',fontsize=20)\n",
    "pl.subplot(122)\n",
    "pl.xlabel(r'$\\theta$ [-]',fontsize=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}