{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulation of ground temperature profile evolution\n",
    "\n",
    "*This notebook creates and runs a one-dimensional simulation model of temperature profile evolution under sinusoidal forcing. It assumes heat transport occurs only by conduction (and not, for example, by water infiltration or groundwater flow).*\n",
    "\n",
    "(Author: Greg Tucker, University of Colorado, Boulder. Latest update November 2021.)\n",
    "\n",
    "## How do I run the model?\n",
    "\n",
    "First, run the complete notebook (Cell => Run All).\n",
    "\n",
    "Next, page down and read through the sections \"An example model run with default parameters\" and \"An example of changing parameters to simulate seasonal temperature cycles.\" These sections provide examples of how to run the model and plot its output.\n",
    "\n",
    "Then, continue on to the section \"Your turn\". Use the cells below that section to type your commands.\n",
    "\n",
    "## What to see and do?\n",
    "\n",
    "- Watch the movies in the example runs. The movies show a temperature profile through a hypothetical column of soil. With diurnal temperature variation, how deep do you have to go before the temperature variations are only half of what they are at the surface?\n",
    "\n",
    "- See the example of seasonal forcing. With seasonal forcing, how deep do you have to go before variations are only about half of what they are at the surface?\n",
    "\n",
    "- What happens when the thermal diffusivity is twice its default value? What if it is half the default value?\n",
    "\n",
    "\n",
    "## What's under the hood?\n",
    "\n",
    "The computer code implements a finite-difference (numerical) solution to the heat diffusion equation. The heat diffusion is what you get when you combine two pieces. The first piece is Fourier's law, which says that the flow of heat energy by conduction, $Q$ (in Watts per square meter), equals the product of thermal conductivity, $k$ (Watts per meter per Kelvin), and the temperature gradient, $dT/dz$:\n",
    "\n",
    "$Q = -k \\frac{\\partial T}{\\partial z}$,\n",
    "\n",
    "where $T$ is temperature and $z$ is depth.\n",
    "\n",
    "The second piece is conservation of heat energy, which in 1D (a depth profile) says:\n",
    "\n",
    "$\\frac{\\partial T}{\\partial t} = - \\frac{1}{\\rho C} \\frac{\\partial Q}{\\partial z}$,\n",
    "\n",
    "where $\\rho$ is the density of the soil/rock, and $C$ is the material's heat capacity (Joules per kilogram per Kelvin).\n",
    "\n",
    "Combining these two yields the 1D heat diffusion equation (with no sources/sinks):\n",
    "\n",
    "$\\boxed{\\frac{\\partial T}{\\partial t} = \\kappa \\frac{\\partial^2 T}{\\partial z^2}}$.\n",
    "\n",
    "Here $\\kappa$ is the *thermal diffusivity*, defined as:\n",
    "\n",
    "$\\kappa = \\frac{k}{\\rho C}$.\n",
    "\n",
    "A typical thermal diffusivity for soil or rock might be on the order of $10^{-6}$ m$^2$/s.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import rc\n",
    "from matplotlib.animation import FuncAnimation\n",
    "from IPython.display import HTML"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The block below defines the code for the model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "class HeatConductionSimulator(object):\n",
    "    \n",
    "    _ONE_DAY = 24.0 * 60.0 * 60.0\n",
    "    _ONE_YEAR = 365.25 * _ONE_DAY\n",
    "\n",
    "    def __init__(self,\n",
    "                 thermal_diffusivity=1.0e-6,\n",
    "                 surface_temp=-5.0,\n",
    "                 amplitude=10.0,\n",
    "                 period=24.0 * 3600.0,\n",
    "                 profile_depth=1.0,\n",
    "                 geotherm=0.025,\n",
    "                 num_nodes=101,\n",
    "                 diffusion_number=0.2,\n",
    "                ):\n",
    "        \"\"\"Initialize the HeatConductionSimulator.\n",
    "\n",
    "        Parameters\n",
    "        ----------\n",
    "        thermal_diffusivity: float\n",
    "            Thermal diffusivity, in m2/s (default 10^-6)\n",
    "        surface_temp: float\n",
    "            Mean surface temperature, degrees C (default -5)\n",
    "        amplitude: float\n",
    "            Amplitude of daily temperature excursions from the mean, degrees C (default 10)\n",
    "        period: float\n",
    "            Period of temperature variation, s (default 1 day = 86,400 s)\n",
    "        profile_depth: float\n",
    "            Depth of profile, m (default 1 m)\n",
    "        geotherm: float\n",
    "            Background geothermal gradient (default 0.025 deg C / m)\n",
    "        num_nodes: int\n",
    "            Number of nodes (default 101)\n",
    "        diffusion_number: float\n",
    "            Dimensionless diffusion number = D dt / dz^2 (default 0.2)\n",
    "        \"\"\"\n",
    "        self.profile_depth = profile_depth\n",
    "        self.surface_temp = surface_temp\n",
    "        self.amplitude = amplitude\n",
    "        self.period = period\n",
    "        self.alpha = diffusion_number  # must be <= 1/2 for numerical stability\n",
    "        self.dz = profile_depth / (num_nodes - 1)\n",
    "        self.timestep_duration = self.alpha * self.dz * self.dz / thermal_diffusivity\n",
    "\n",
    "        self.depth = np.linspace(0, profile_depth, num_nodes)\n",
    "        self.temp = surface_temp + geotherm * self.depth\n",
    "        self.current_time = 0.0\n",
    "        self.max_act_layer_depth = 0.0    \n",
    "        \n",
    "    def run_one_step(self):\n",
    "        \"\"\"Advance for one time step\"\"\"\n",
    "        self.temp[0] = (self.surface_temp + \n",
    "                        self.amplitude * np.sin(2.0 * np.pi * self.current_time / self.period))\n",
    "        self.temp[1:-1] += self.alpha * (self.temp[2:] - 2 * self.temp[1:-1] + self.temp[:-2])\n",
    "        self.current_time += self.timestep_duration\n",
    "\n",
    "    def run_n_steps(self, n):\n",
    "        for i in range(n):\n",
    "            self.run_one_step()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The function below initializes and runs the model, and generates and displays and animation of the result. Here, we are using default parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_model(num_days=2.0,\n",
    "              thermal_diffusivity=1.0e-6,\n",
    "              surface_temp=-5.0,\n",
    "              amplitude=10.0,\n",
    "              period=24.0 * 3600.0,\n",
    "              profile_depth=1.0,\n",
    "              geotherm=0.025,\n",
    "              num_nodes=101,\n",
    "              save_every=100,\n",
    "              diffusion_number=0.2,\n",
    "             ):\n",
    "    \"\"\"Initialize, run, and display output from model.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    num_days: float\n",
    "        Number of days to run\n",
    "    thermal_diffusivity: float\n",
    "        Thermal diffusivity, in m2/s (default 10^-6)\n",
    "    surface_temp: float\n",
    "        Mean surface temperature, degrees C (default -5)\n",
    "    amplitude: float\n",
    "        Amplitude of daily temperature excursions from the mean, degrees C (default 10)\n",
    "    period: float\n",
    "        Period of temperature variation, s (default 1 day = 86,400 s)\n",
    "    profile_depth: float\n",
    "        Depth of profile, m (default 1 m)\n",
    "    geotherm: float\n",
    "        Background geothermal gradient (default 0.025 deg C / m)\n",
    "    num_nodes: int\n",
    "        Number of nodes (default 101)\n",
    "    save_every: int\n",
    "        Interval to save an animation frame, in iterations (default 100)\n",
    "    diffusion_number: float\n",
    "        Dimensionless diffusion number = D dt / dz^2 (default 0.2)\n",
    "    \"\"\"\n",
    "\n",
    "    # Instantiate and initialize a simulator\n",
    "    hcs = HeatConductionSimulator(thermal_diffusivity=thermal_diffusivity,\n",
    "                                  profile_depth=profile_depth,\n",
    "                                  geotherm=geotherm,\n",
    "                                  surface_temp=surface_temp,\n",
    "                                  amplitude=amplitude,\n",
    "                                  period=period,\n",
    "                                  num_nodes=num_nodes,\n",
    "                                  diffusion_number=diffusion_number,\n",
    "                                 )\n",
    "    \n",
    "    # Calculate number of animation iterations\n",
    "    nsteps = int(num_days * 24.0 * 3600.0 / (hcs.timestep_duration * save_every))\n",
    "    \n",
    "    # Set up a blank figure with placeholder lists for data\n",
    "    fig, ax = plt.subplots()\n",
    "    xdata = []\n",
    "    ydata = []\n",
    "    obj = ax.plot([], [], color = 'k')\n",
    "    \n",
    "    # Then, set up an initialization function\n",
    "    def init():\n",
    "        ax.set_ylim(0, hcs.profile_depth)\n",
    "        ax.set_xlim(hcs.surface_temp - hcs.amplitude, hcs.surface_temp + hcs.amplitude)\n",
    "        ax.set_ylabel('Depth (m)')\n",
    "        ax.set_xlabel('Temperature (degrees C)')\n",
    "        return(obj)\n",
    "    \n",
    "    # Next, define the update function\n",
    "    def update(i):\n",
    "        ax.cla()\n",
    "        hcs.run_n_steps(save_every)\n",
    "        xdata = hcs.temp\n",
    "        ydata = hcs.depth\n",
    "        ax.set_ylim(0, hcs.profile_depth)\n",
    "        xmin = hcs.surface_temp - hcs.amplitude\n",
    "        xmax = hcs.surface_temp + hcs.amplitude\n",
    "        ax.set_xlim(xmin, xmax)\n",
    "        ax.set_ylabel('Depth (m)')\n",
    "        ax.set_xlabel('Temperature (degrees C)')\n",
    "        ax.invert_yaxis()\n",
    "        obj = ax.plot(xdata, ydata, color = 'r')\n",
    "        obj = ax.plot([0.0, 0.0], [0.0, hcs.profile_depth], 'k--')\n",
    "        if np.amax(hcs.temp) > 0.0:\n",
    "            above_zero = np.where(hcs.temp > 0.0)[0]\n",
    "            z_active_layer = hcs.depth[above_zero[-1]]\n",
    "            obj = ax.plot([xmin, xmax], [z_active_layer, z_active_layer], 'b')\n",
    "            if z_active_layer > hcs.max_act_layer_depth:\n",
    "                hcs.max_act_layer_depth = z_active_layer\n",
    "        return(obj)\n",
    "    \n",
    "    # Run the animation!\n",
    "    anim = FuncAnimation(fig, update, nsteps, init_func=init, blit = True)\n",
    "    \n",
    "    # Convert the animation to HTML\n",
    "    vid = HTML(anim.to_jshtml())\n",
    "    \n",
    "    # Report maximum active layer depth\n",
    "    print('Maximum active layer depth: ' + str(hcs.max_act_layer_depth) + ' meters')\n",
    "    \n",
    "    return vid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## An example model run with default parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_movie = run_model()\n",
    "my_movie"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## An example of changing parameters to simulate seasonal temperature cycles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "my_movie = run_model(num_days=365.25*2, amplitude=15., period=365.25*24*3600.0, profile_depth=12.,\n",
    "                     save_every=500)\n",
    "my_movie"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Your turn\n",
    "\n",
    "Use the cells below to experiment with running the model with different parameters and/or boundary conditions, as indicated in the lab exercise. To run the model and produce a movie, use `my_movie = run_model(`*keyword parameters*`)`. The main parameters that you might want to explore include:\n",
    "\n",
    "-    num_days: Number of days to run (default 2)\n",
    "-    thermal_diffusivity: Thermal diffusivity, in m2/s (default 10^-6)\n",
    "-    surface_temp: Mean surface temperature, degrees C (default -5)\n",
    "-    amplitude: Amplitude of daily temperature excursions from the mean, deg C (default 10)\n",
    "-    period: Period of temperature variation, s (default 1 day = 86,400 s)\n",
    "-    profile_depth: Depth of profile, m (default 1 m)\n",
    "\n",
    "For example, to run a model with a thermal diffusivity of $2\\times 10^6$ m$^2$/s and a mean surface temperature of -10 $^\\circ$C, you would create and run a new cell with the code:\n",
    "\n",
    "`my_movie = run_model(thermal_diffusivity=2.0e-6, surface_temp=-10.0)`\n",
    "\n",
    "`my_movie`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}