{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 1D simulation of a valley glacier\n", "\n", "This notebook creates and runs a one-dimensional simulation model of a hypothetical valley glacier." ] }, { "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 ValleyGlacierSimulator(object):\n", "\n", " def __init__(self,\n", " valley_length=10000.0,\n", " valley_slope=0.1,\n", " ela=900.0,\n", " gamma=0.01,\n", " flow_law_coefficient=6.8e-24,\n", " ice_density=850.0,\n", " sliding_factor=0.0,\n", " grav_accel=9.8,\n", " num_nodes=50,\n", " timestep_duration=0.04,\n", " ):\n", " \"\"\"Initialize the ValleyGlacierSimulator.\n", "\n", " Parameters\n", " ----------\n", " valley_length: float\n", " Length of valley, m (default 10 km)\n", " valley_slope: float\n", " Valley slope gradient, m/m\n", " ela: float\n", " Equilibrium line altitude, m (default 900)\n", " gamma: float\n", " Mass balance coefficient, 1/yr (default 0.01)\n", " flow_law_coefficient: float\n", " \"A\" in flow law, s^-1 Pa^-3 (default 6.8 x 10^-24)\n", " (Note: 6.8 x 10^-24 s^-1 Pa^-3 ~ 2 x 10^-16 yr^-1 Pa^3)\n", " ice_density: float\n", " Bulk density of glacial ice, kg/m3 (default 850)\n", " sliding_factor: float\n", " Sliding speed as proportion of deformation speed (default 0)\n", " grav_accel: float\n", " Gravitational acceleration (default 9.8)\n", " num_nodes: int\n", " Number of nodes (default 50)\n", " timestep_duration: float\n", " Duration of a time step, yr (default 0.04)\n", " \"\"\"\n", " self.valley_length = valley_length\n", " self.valley_slope = valley_slope\n", " self.ela = ela\n", " self.gamma = gamma\n", " self.flow_law_coefficient = flow_law_coefficient\n", " self.ice_density = ice_density\n", " self.sliding_factor = sliding_factor\n", " self.g = grav_accel\n", " self.timestep_duration = timestep_duration\n", " self.dx = valley_length / (num_nodes - 1)\n", " \n", " self.current_time = 0.0\n", "\n", " self.x = np.linspace(0.5 * self.dx, \n", " valley_length - 0.5 * self.dx,\n", " num_nodes)\n", " self.x_edges = np.linspace(0, valley_length, num_nodes + 1)\n", " zmax = self.valley_slope * self.valley_length\n", " self.z = zmax - self.valley_slope * self.x\n", " self.ice_thickness = np.zeros(num_nodes)\n", " self.q = np.zeros(num_nodes + 1)\n", " self.ice_slope = np.zeros(num_nodes)\n", " self.ice_slope[-1] = valley_slope\n", " self.b = gamma * (self.z - ela)\n", "\n", " self.Hstar = 0.01 # depth decay scale for melt rate, m\n", " \n", " # Calculate lumped parameter in flow law: U = k H^4, in m/(s m4)\n", " #\n", " # The units here are: \n", " # A => 1 / s Pa3 (= m3 s5 / kg3)\n", " # rho_ice => kg / m3\n", " # g => m / s2\n", " #\n", " # Therefore A (rho_ice g)^3 => m3 s5 kg3 m3 / kg3 m9 s6 => 1 / m3 s\n", " # and k_sec H^4 => m4 / m3 s => m / s\n", " # and k H^4 => m4 s / m3 s yr => m / yr\n", " #\n", " k_sec = (0.4 * flow_law_coefficient * (1 + sliding_factor)\n", " * (ice_density * grav_accel)**3) \n", " self.k = k_sec * 3600 * 24 * 365.25; # To per-year from per-sec\n", "\n", "\n", " def run_one_step(self):\n", " \"\"\"Advance for one time step\"\"\"\n", " #print('Time = ' + str(self.current_time))\n", " self.ice_slope[:-1] = -(np.diff(self.ice_thickness + self.z)\n", " / self.dx)\n", " self.q[1:] = self.k * self.ice_slope**3 * self.ice_thickness**5\n", " dqdx = np.diff(self.q) / self.dx # flux gradient, m/yr\n", " self.b[:] = self.gamma * ((self.z + self.ice_thickness) - self.ela)\n", " b_eff = self.b # effective net accum rate, m/yr\n", " below_ela = self.b < 0\n", " b_eff[below_ela] = (self.b[below_ela]\n", " * (1 - np.exp(-self.ice_thickness[below_ela]\n", " / self.Hstar)))\n", " dHdt = b_eff - dqdx # mass balance\n", " self.ice_thickness += dHdt * self.timestep_duration\n", " self.ice_thickness = np.maximum(self.ice_thickness, 0.0)\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()\n" ] }, { "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(run_duration=1000.0,\n", " valley_length=10000.0,\n", " valley_slope=0.1,\n", " ela=900.0,\n", " gamma=0.01,\n", " flow_law_coefficient=6.8e-24,\n", " ice_density=850.0,\n", " sliding_factor=0.0,\n", " grav_accel=9.8,\n", " num_nodes=50,\n", " timestep_duration=0.04,\n", " save_every=400,\n", " ):\n", " \"\"\"Initialize, run, and display output from model.\n", "\n", " Parameters\n", " ----------\n", " run_duration: float\n", " Number of years to run (default 1200)\n", " save_every: int\n", " Interval to save an animation frame, in iterations (default 5)\n", " (see ValleyGlacierSimulator for other parameters)\n", " \"\"\"\n", "\n", " # Instantiate and initialize a simulator\n", " model = ValleyGlacierSimulator(valley_length=valley_length,\n", " valley_slope=valley_slope,\n", " ela=ela,\n", " gamma=gamma,\n", " flow_law_coefficient=flow_law_coefficient,\n", " ice_density=ice_density,\n", " sliding_factor=sliding_factor,\n", " grav_accel=grav_accel,\n", " num_nodes=num_nodes,\n", " timestep_duration=timestep_duration\n", " )\n", "\n", " # Calculate number of animation iterations\n", " nsteps = int(run_duration / (model.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, 2 * np.amax(model.z))\n", " ax.set_xlim(0, model.valley_length)\n", " ax.set_ylabel('Height (m)')\n", " ax.set_xlabel('Distance (m)')\n", " return(obj)\n", "\n", " # Next, define the update function\n", " def update(i):\n", " ax.cla()\n", " model.run_n_steps(save_every)\n", " xdata = model.x\n", " ydata = model.z\n", " ax.set_ylim(0, 2 * np.amax(model.z))\n", " ax.set_xlim(0, model.valley_length)\n", " ax.set_ylabel('Height (m)')\n", " ax.set_xlabel('Distance (m)')\n", " ax.set_title('Time = ' + str(round(model.current_time)) + ' years')\n", " ax.text(0.05 * model.valley_length, 1.8 * np.amax(model.z), 'ACCUMULATION')\n", " ax.text(0.82 * model.valley_length, 0.25 * np.amax(model.z), 'ABLATION')\n", " ax.text(0.9 * model.valley_length, 1.02 * model.ela, 'ELA')\n", " obj = ax.plot(xdata, ydata + model.ice_thickness, color='c')\n", " obj = ax.plot(xdata, ydata, color='k')\n", " obj = ax.plot([0.0, model.valley_length], [model.ela, model.ela], 'm:')\n", " return(obj)\n", "\n", " def report_values(model):\n", " \"\"\"Report values of various quantities.\"\"\"\n", " print('Final glacier length: ' + str(model.dx \n", " * np.count_nonzero(model.ice_thickness > 0.0))\n", " + ' m')\n", " print('Final maximum thickness: ' + str(np.amax(model.ice_thickness)) + ' m')\n", " qq = model.q[1:]\n", " u = 0 * qq\n", " ice_present = model.ice_thickness > 0.0\n", " u[ice_present] = qq[ice_present] / model.ice_thickness[ice_present]\n", " print('Final maximum velocity: ' + str(np.amax(u)) + ' m/yr')\n", " print('Final maximum ice specific discharge: ' + str(np.amax(qq)) + ' m2/yr') # Clean up the model\n", " print('Final ice volume (per width): ' + str(model.dx * np.sum(model.ice_thickness))\n", " + ' m2')\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()) # or HTML(anim.to_jshtml())\n", "\n", " report_values(model)\n", "\n", " return vid, model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## An example model run with default parameters" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "my_movie, model = run_model()\n", "my_movie" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## An example model run with a different ELA" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "my_movie, model = run_model(ela=750.0)\n", "my_movie" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Example of plotting the ice thickness discharge\n", "\n", "The `model` data object contains arrays for ice thickness and discharge at the end of the run. To plot ice thickness, use the `model.x` array for the locations of the model cells and the `model.ice_thickness` for the ice thickness at those locations:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.plot(model.x, model.ice_thickness, 'c')\n", "plt.xlabel('Distance (m)')\n", "plt.ylabel('Thickness (m)')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "To plot ice discharge, use `model.x_edges` (the locations of the *edges* of the cells) and `model.q`, which is the ice discharge calculated *at the cell edges*." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.plot(model.x_edges, model.q, 'b')\n", "plt.xlabel('Distance (m)')\n", "plt.ylabel('Specific ice discharge (m2/yr)')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Your turn\n", "\n", "Use new code cells below to experiment with running the model with different parameters and/or boundary conditions, as indicated in the assignment. To run the model and produce a movie, use `my_movie = run_model(`*keyword parameters*`)`, as in the example above. See the header of the `ValleyGlacierSimulator` class for a list of parameters and their default values." ] } ], "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 }