Idealized-Guide.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Notebook guide: linear viscoelastic ice cauldrons"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Welcome!  This notebook will guide you through some of the key functions included in `cauldron_funcs.py` to simulate linear elastic/viscoelastic collapse of ice cauldrons.  Here we will use an idealized example.  \n",
    "\n",
    "The subfolder `ESkafta-2015` applies these methods to the 2015 collapse of Eastern Skafta Cauldron, Vatnajokull, Iceland.  However, ArcticDEM and GPS data used to initialize and assess that example come from other sources and are not included in this repository."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import statements"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we import the packages and modules we will use.  Remember to have your `jokull` environment (defined in `environment.yml` in this repository) activated to use these packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.misc as scp\n",
    "from scipy import interpolate\n",
    "from scipy.ndimage import gaussian_filter\n",
    "# from osgeo import gdal\n",
    "from netCDF4 import Dataset\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import cm\n",
    "from matplotlib.patches import Rectangle\n",
    "import math\n",
    "import sys\n",
    "sys.path.insert(0, '/Users/lizz/Documents/GitHub/VE-cauldrons')\n",
    "from cauldron_funcs import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Setting up cauldron geometry"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you have surface elevation data from a real cauldron, you can read it in as done in `ESkafta-2015/Skafta-ArcticDEM-transecting.py`.  Here we will initialize an idealized cauldron."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "transect_length = 2000 # m, diameter of the cauldron\n",
    "npoints = 500\n",
    "x_cylcoords = np.linspace(-0.5*transect_length, 0.5*transect_length, num=npoints) # x-axis values in cylindrical coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "initial_elevation = 1200 # m, surface elevation of the cauldron before collapse\n",
    "ice_thickness = 300 # m\n",
    "cldrn = Cauldron(name='Example_Cauldron', \n",
    "                 initial_surface = lambda x: initial_elevation, # function relating initial elevation to x-coordinate\n",
    "                 thickness = ice_thickness,\n",
    "                 radius = 0.5*transect_length\n",
    "                )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cldrn.set_viscoelastic_bendingmod()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute collapse profiles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we will use functions of our example `cldrn`, an instance of the `Cauldron` class, to compute post-collapse profiles of an elastic and linear viscoelastic plate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Elastic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "elastic_profile = [cldrn.LL_profile(x) for x in x_cylcoords] # LL for 'Landau & Lifshitz', analytical solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "initial_surf = np.full(shape=npoints, fill_value=initial_elevation)\n",
    "\n",
    "fig, ax = plt.subplots(1, figsize=(7, 3))\n",
    "ax.plot(x_cylcoords, initial_surf, color='k', ls='-.', label='Pre-collapse surface')\n",
    "ax.plot(x_cylcoords, elastic_profile, color='DarkBlue', lw=2, label='Elastic plate')\n",
    "ax.legend(bbox_to_anchor=(1.05, 1), loc='upper left')\n",
    "ax.set(aspect=5,\n",
    "       xlim=(-0.5*transect_length, 0.5*transect_length),\n",
    "       ylim=(initial_elevation-100, initial_elevation+100),\n",
    "       xlabel='Distance from cauldron center [m]',\n",
    "       ylabel='Surface elevation [m a.s.l.]', \n",
    "       title='Elastic collapse with E={:.1E} Pa'.format(cldrn.youngmod))\n",
    "ax.fill_between(x_cylcoords, initial_surf, elastic_profile, color='Gainsboro', hatch='/', edgecolor='DimGray', linewidth=0, alpha=0.7)\n",
    "ax.fill_between(x_cylcoords, elastic_profile, (ax.get_ylim()[0]), color='Azure')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Viscoelastic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Within `cauldron-funcs.py`, we apply a Laplace transform to take our analytical solution for the elastic case to the linear viscoelastic case.  This will take a little longer to run.  We will select only a few time slices to view, to minimize computational demand."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "days_of_collapse = 4 # how many days of collapse you want to simulate\n",
    "nseconds = days_of_collapse*24*60*60\n",
    "timestep = 40000 # seconds, approximately half a day\n",
    "times = np.arange(0, nseconds, step=timestep)\n",
    "\n",
    "VE_profiles = {t: [] for t in times}\n",
    "for t in times:\n",
    "    VE_profiles[t] = [cldrn.viscoelastic_profile(x, t) for x in x_cylcoords]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "colors = cm.get_cmap('winter_r')(np.linspace(0.1, 0.9, num=len(times)+1))\n",
    "ylim_lower = float(np.min(VE_profiles[times[-1]])) - 100 # ensure display captures full collapse range\n",
    "\n",
    "fig, ax1 = plt.subplots(1, figsize=(7, 3))\n",
    "ax1.plot(x_cylcoords, initial_surf, color='k', ls='-.') #, label='15 Oct 2012'\n",
    "ax1.plot(x_cylcoords, elastic_profile, color=colors[0], lw=2, label='Elastic plate')\n",
    "for i,ti in enumerate(times[::2]):\n",
    "    labeltime = int(round(ti/86400)) #time in days\n",
    "    ax1.plot(x_cylcoords, VE_profiles[ti], ls='--', color=colors[i+1], lw=2, label='Viscoelastic, t = {} days'.format(labeltime))\n",
    "ax1.legend(bbox_to_anchor=(1.05,1), loc='upper left') #\n",
    "ax1.set(aspect=5,\n",
    "        ylim=(ylim_lower, initial_elevation+100),\n",
    "        xlim=(-0.5*transect_length, 0.5*transect_length),\n",
    "        xlabel='Distance from cauldron center [m]',\n",
    "        title='Viscoelastic collapse with E={:.1E} Pa, eta={:.1E} Pa s'.format(cldrn.youngmod, cldrn.dyn_viscos))\n",
    "ax1.fill_between(x_cylcoords, y1=initial_surf, y2=[float(ve) for ve in VE_profiles[times[-1]]], color='Gainsboro', hatch='/', edgecolor='DimGray', linewidth=0, alpha=0.7)\n",
    "ax1.fill_between(x_cylcoords, y1=[float(ve) for ve in VE_profiles[times[-1]]], y2=ylim_lower, color='Azure')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computing stress"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We compute an approximation to the tensile/compressive stress at the surface of the ice.  Note that our approach here makes several assumptions that cause this approximation to be an _overestimate_: it will return a maximum plausible stress that could have been attained during collapse with this rheology.  See Ultee et al 2020, Journal of Glaciology for more details."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "elas_plate_stress = [cldrn.elastic_stress(x, config='radial_plate') for x in x_cylcoords]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax2 = plt.subplots()\n",
    "ax2.axhline(0, color='k', ls=':')\n",
    "ax2.plot(x_cylcoords, 1e-6*np.asarray(elas_plate_stress))\n",
    "ax2.fill_between(x_cylcoords, 1e-6*np.asarray(elas_plate_stress), 0, where=[sigma>0 for sigma in elas_plate_stress], color='DarkBlue', alpha=0.5, label='Compressive')\n",
    "ax2.fill_between(x_cylcoords, 1e-6*np.asarray(elas_plate_stress), 0, where=[sigma<0 for sigma in elas_plate_stress], color='FireBrick', alpha=0.5, label='Tensile')\n",
    "ax2.legend(bbox_to_anchor=(1.05,1), loc='upper left')\n",
    "ax2.set(xlabel='Distance from cauldron center',\n",
    "       ylabel='Max surface stress [MPa]',\n",
    "       title='Elastic case')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ve_plate_stress_init = [float(cldrn.viscoelastic_stress(x, times[0])) for x in x_cylcoords]\n",
    "ve_plate_stress_final = [float(cldrn.viscoelastic_stress(x, times[-1])) for x in x_cylcoords]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax3 = plt.subplots()\n",
    "ax3.axhline(0, color='k', ls=':')\n",
    "ax3.plot(x_cylcoords, 1e-6*np.asarray(ve_plate_stress_init))\n",
    "ax3.fill_between(x_cylcoords, 1e-6*np.asarray(ve_plate_stress_init), 0, where=[sigma>0 for sigma in ve_plate_stress_init], color='DarkBlue', alpha=0.5, label='Compressive')\n",
    "ax3.fill_between(x_cylcoords, 1e-6*np.asarray(ve_plate_stress_init), 0, where=[sigma<0 for sigma in ve_plate_stress_init], color='FireBrick', alpha=0.5, label='Tensile')\n",
    "ax3.legend(bbox_to_anchor=(1.05,1), loc='upper left')\n",
    "ax3.set(xlabel='Distance from cauldron center',\n",
    "       ylabel='Max surface stress [MPa]',\n",
    "       title='Viscoelastic case, E={:.1E}, eta={:.1E}'.format(cldrn.youngmod, cldrn.dyn_viscos))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that in this simple linear viscoelastic case, there is no creep, crevassing, etc. that would concentrate or dissipate stress.  Unphysically high values of surface stress are possible toward the end of the collapse period, when such effects would be expected to substantially modify the stress in nature."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Modifying material properties"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `Cauldron` instance we set up above has inherited default material properties from the `Ice` class defined in `cauldron_funcs.py`.  \n",
    "\n",
    "You can examine the material properties (density `rho_ice`, Young's modulus `youngmod`, Poisson's ratio `poisson_nu`, and dynamic viscosity `dyn_viscos`) like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cldrn.youngmod"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also change the material constants by calling them like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "cldrn.youngmod = 2E9 # Pa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now re-run the analysis cells with an altered Young's modulus (or other material property) to explore the parameter space."
   ]
  }
 ],
 "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}