{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Excercises on ice sheet modeling\n", "Adapted from © Excercises on ice sheet modeling at GIA summerschool 2023, Torsten Albrecht, albrecht@pik-potsdam.de, PIK, 2023\n", "\n", "- code snippets are based on [PISM](https://www.pism.io), pism/examples/mismip/mismip2d" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 3 Marine ice sheet instability" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There is a whole theory on MISI since 50 years ([Weertman 1974](https://doi.org/10.3189/S0022143000023327)), in particular for unbuttressed ice sheets and ice shelves. [Schoof 2007](https://doi.org/10.1029/2006JF000664) came up with a semianalytical solution for the variational problem, with the flux at the grounding line\n", "\n", "$$\n", "Q(x_g) \\propto H(x_g)^{4.75}\n", "$$\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Consider a marine ice sheet resting on an overdeepende bed (typical for West Antarctic Ice Sheet):" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [], "source": [ "### Import packages\n", "\n", "%matplotlib inline\n", "%config InlineBackend.figure_formats = {'svg',}\n", "\n", "import matplotlib.pylab as plt\n", "from matplotlib import cm, colors\n", "import numpy as np" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "### Defining ice sheet functions\n", "\n", "# following pism/examples/mismip/mismip2d\n", "\n", "# ice softness\n", "A3b4 = 5.0e-26\n", "#A3b4 = 1.0e-26\n", "\n", "#sliding coefficient\n", "Csl = 7.624e6\n", "\n", "rhoi = 900.0\n", "rhosw = 1000.0\n", "\n", "spa = 60*60*365.25*24 # s per year\n", "\n", "#mode 2\n", "N = 1500\n", "\n", "# length\n", "Lm = 1800e3\n", "x = np.linspace(-Lm, Lm, 2 * N + 1)\n", "\n", "# convert x to kilometers\n", "xkm = x/1e3\n", "\n", "#calving_front\n", "xcf = 1750e3\n", "\n", "def bed(x):\n", " #\"Bed depth below sea level. (-b(x) = topg(x))\"\n", " xx = x / 7.5e5\n", " return -(729. - 2184.8 * xx ** 2. + 1031.72 * xx ** 4. - 151.72 * xx ** 6.)\n", "\n", "def M():\n", " #\"Accumulation rate (m/s)\"\n", " return 0.3 / spa\n", "\n", "def b_slope(x):\n", " #\"\"\"The x-derivative of b(experiment, x).\"\"\"\n", " xx = x / 7.5e5\n", " return -(- 2184.8 * (2. / 7.5e5) * xx\n", " + 1031.72 * (4. / 7.5e5) * xx ** 3.\n", " - 151.72 * (6. / 7.5e5) * xx ** 5.)\n", "\n", "def g():\n", " #\"\"\"Acceleration due to gravity. (Table 2 on page 19 of mismip_4.pdf\n", " #uses this value, i.e. g = 9.8 m s-2.)\"\"\"\n", " return 9.81\n", "\n", "def m():\n", " #\"Sliding law exponent\"\n", " return 1 / 3.0\n", " #return 1.0\n", "\n", "def n():\n", " #\"Glen exponent\"\n", " return 3.0\n", "\n", "\n", "def cold_function(x, sl, A0 , theta=0.0):\n", " \"\"\"Evaluates function whose zeros define x_g in 'cold' steady marine sheet problem.\"\"\"\n", " r = rhoi / rhosw\n", " h_f = r ** (-1.) * (bed(x)+sl)\n", " b_x = b_slope(x)\n", " s = M() * x\n", " rho_g = rhoi * g()\n", " return (theta * M()\n", " + Csl * s ** (m() + 1.0) / (rho_g * h_f ** (m() + 2.))\n", " - theta * s * b_x / h_f\n", " - A0 * (rho_g * (1.0 - r) / 4.0) ** n() * h_f ** (n() + 1.0))\n", "\n", "\n", "def x_g(sl, xi, A0, theta=0.0):\n", " \"\"\"Computes the theoretical grounding line location using Newton's method.\"\"\"\n", "\n", " # set the initial guess\n", " x = xi*1e3 #m\n", "\n", " delta_x = 10. # Finite difference step size (metres) for gradient calculation\n", " tolf = 1.e-4 # Tolerance for finding zeros\n", " eps = np.finfo(float).eps\n", " normf = tolf + eps\n", " toldelta = 1.e1 # Newton step size tolerance\n", " dx = toldelta + 1.0\n", "\n", " # this is just a shortcut\n", " def F(x):\n", " return cold_function(x, sl, A0, theta)\n", "\n", " while (normf > tolf) or (abs(dx) > toldelta):\n", " f = F(x)\n", " normf = abs(f)\n", " grad = (F(x + delta_x) - f) / delta_x\n", " dx = -f / grad\n", " x = x + dx\n", "\n", " return x\n", "\n", "\n", "def thickness(sl, x, xi, A0, theta=0.0):\n", " #\"\"\"Compute ice thickness for x > 0.\n", " #\"\"\"\n", " # compute the grounding line position\n", " xg = x_g(sl, xi, A0, theta)\n", "\n", " def surface(h, x):\n", " b_x = b_slope(x)\n", " rho_g = rhoi * g()\n", " s = M() * np.abs(x)\n", " return b_x - (Csl / rho_g) * s ** m() / h ** (m() + 1)\n", "\n", " # extract the grounded part of the grid\n", " x_grounded = x[x < xg]\n", "\n", " # We will integrate from the grounding line inland. odeint requires that\n", " # the first point in x_grid be the one corresponding to the initial\n", " # condition; append it and reverse the order.\n", " x_grid = np.append(xg, x_grounded[::-1])\n", "\n", " # use thickness at the grounding line as the initial condition\n", " h_f = bed(xg) * rhosw / rhoi\n", "\n", " import scipy.integrate\n", " thk_grounded = scipy.integrate.odeint(surface, [h_f], x_grid, atol=1.e-9, rtol=1.e-9)\n", "\n", " # now 'result' contains thickness in reverse order, including the grounding\n", " # line point (which is not on the grid); discard it and reverse the order.\n", " thk_grounded = np.squeeze(thk_grounded)[:0:-1]\n", "\n", " # extract the floating part of the grid\n", " x_floating = x[x >= xg]\n", "\n", " # compute the flux through the grounding line\n", " q_0 = M() * xg\n", "\n", " # Calculate ice thickness for shelf from van der Veen (1986)\n", " r = rhoi / rhosw\n", " rho_g = rhoi * g()\n", " numer = h_f * (q_0 + M() * (x_floating - xg))\n", " base = q_0 ** (n() + 1) + h_f ** (n() + 1) * ((1 - r) * rho_g / 4) ** n() * A0 \\\n", " * ((q_0 + M() * (x_floating - xg)) ** (n() + 1) - q_0 ** (n() + 1)) / M()\n", " thk_floating = numer / (base ** (1.0 / (n() + 1)))\n", "\n", " x_nonnegative = x[x >= 0]\n", " thk_nonnegative = np.r_[thk_grounded, thk_floating]\n", " thk_nonnegative = thk_nonnegative[x >= 0]\n", " thk_nonnegative[x_nonnegative > xcf] = 0\n", " \n", " thk = np.zeros_like(x)\n", " thk[x >= 0] = thk_nonnegative\n", " thk[x < 0] = thk_nonnegative[:0:-1]\n", " \n", " return thk\n", "\n" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "### Define plotting function\n", "def plot_mismip2d(topg,thk):\n", "\n", " # grounded / floating mask\n", " mask = np.ones_like(topg)*2 #grounded\n", " mask[-rhoi/rhosw*thk>topg]=3 # floating\n", " mask[thk==0]=4 # ice free ocean\n", "\n", " # surface elevation\n", " usurf = topg + thk\n", " usurf[mask>=3] = (1.0-rhoi/rhosw) * thk[mask>=3]\n", " usurf = np.ma.array(usurf, mask=mask == 4)\n", "\n", " # sub surface elevation\n", " lsurf = topg.copy()\n", " lsurf[mask == 3] = -rhoi / rhosw * thk[mask == 3]\n", " lsurf = np.ma.array(lsurf, mask=mask == 4)\n", "\n", " #######################################\n", "\n", "\n", " fig = plt.figure(12,figsize=(7,2))\n", " ax = fig.add_subplot(1,1,1)\n", "\n", " # sea level 0\n", " ax.plot(xkm, np.zeros_like(xkm), ls='dotted', color='C9')\n", " ax.plot(xkm, usurf, color='C0')\n", "\n", " # bedrock\n", " #ax.plot(xkm, topg, color='black')\n", " ax.fill_between(xkm,-1500,topg,color='C7')\n", "\n", " # ocean water\n", " ax.plot(xkm, lsurf, color='C1')\n", " ax.fill_between(xkm,topg,lsurf,color='C0')\n", " ax.fill_between(xkm[ ((mask==4) & (xkm<0)) ],topg[ ((mask==4) & (xkm<0)) ],0,color='C0')\n", " ax.fill_between(xkm[ ((mask==4) & (xkm>0)) ],topg[ ((mask==4) & (xkm>0)) ],0,color='C0')\n", "\n", " ax.set_xlabel('distance from the divide, km')\n", " ax.set_ylabel('elevation, m')\n", "\n", " #_, _, ymin, ymax = ax.axis(xmin=xkm.min(), xmax=xkm.max())\n", " _, _, ymin, ymax = ax.axis(xmin=0.0, xmax=xkm.max())\n", " ax.set_ylim(-1500,5000)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercise\n", "\n", "1. Try out different initial guesses for the grounding line position (xi). Observe where the real grounding line position appears (x_g) compared to the guessed position and how the ice sheet thickness changes.\n", "\n", " Can you identify stable grounding line positions (when guessed and real positions are similar) and unstable positions?\n", "\n", " And how does the grounding line position relate ice sheet thickness?\n", "\n", "2. Play around with the global mean sea level (sl) and observe how it changes the grounding line position (x_g).\n", "\n", " Which sea level rise would destabilize the ice sheet?" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "x_g: 1412.2 km\n" ] }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", " \n", " 2026-06-18T15:13:08.609151\n", " image/svg+xml\n", " \n", " \n", " Matplotlib v3.10.8, https://matplotlib.org/\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# change these:\n", "sl = 0.0 # sea level\n", "xi = 1370.0 # initial gues for grounding line position\n", "\n", "# leave these as they are:\n", "topg = -bed(x)-sl\n", "thk = thickness(sl, x, xi, A3b4 , theta=0.0)\n", " \n", "print('x_g:',np.around(x_g(sl, xi, A3b4 , theta=0.0)*1e-3,decimals=1),'km')\n", "\n", "plot_mismip2d(topg,thk)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "x_g: 974.2 km\n" ] }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", " \n", " 2026-06-18T15:13:37.244201\n", " image/svg+xml\n", " \n", " \n", " Matplotlib v3.10.8, https://matplotlib.org/\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sl = 0.0\n", "xi = 1200#1370.0\n", "\n", "topg = -bed(x)-sl\n", "thk = thickness(sl, x, xi, A3b4 , theta=0.0)\n", " \n", "print('x_g:',np.around(x_g(sl, xi, A3b4 , theta=0.0)*1e-3,decimals=1),'km')\n", "\n", "plot_mismip2d(topg,thk)\n" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "x_g: 1296.7 km\n" ] }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", " \n", " 2026-06-18T15:38:44.907742\n", " image/svg+xml\n", " \n", " \n", " Matplotlib v3.10.8, https://matplotlib.org/\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "sl = 160\n", "xi = 1370.0\n", "\n", "topg = -bed(x)-sl\n", "thk = thickness(sl, x, xi, A3b4 , theta=0.0)\n", " \n", "print('x_g:',np.around(x_g(sl, xi, A3b4 , theta=0.0)*1e-3,decimals=1),'km')\n", "\n", "plot_mismip2d(topg,thk)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Optional: \n", "\n", "Ice softness can also influence the stability of the ice sheet. At which ice softness the ice sheet would tip from a smaller state to to more extended one, and for which other way around?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "
\n", "A3b: click to show solution\n", "\n", "If ice cools (becoming harder, e.g $2.5 \\times 10^{-26} Pa^{-3} s^{-1}$) this specific (unbuttressed) marine ice sheet may tip to to more extended state.\n", "If the ice warms (becoming an order of magnitude softer) (e.g $2.5 \\times 10^{-25} Pa^{-3} s^{-1}$) it may tip to a less extended state (MISI).\n", " \n", "
\n", "
" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/tmp/ipykernel_303/2793940455.py:32: RuntimeWarning: invalid value encountered in power\n", " qschoof *= h_gr**((m()+n()+3.0)/(m()+1.0)) * spa*1e-6\n", "/tmp/ipykernel_303/2922672300.py:67: RuntimeWarning: invalid value encountered in scalar power\n", " + Csl * s ** (m() + 1.0) / (rho_g * h_f ** (m() + 2.))\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/svg+xml": [ "\n", "\n", "\n", " \n", " \n", " \n", " \n", " 2026-06-18T15:05:13.832503\n", " image/svg+xml\n", " \n", " \n", " Matplotlib v3.10.8, https://matplotlib.org/\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n" ], "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig = plt.figure(13,figsize=(10,4))\n", "ax = fig.add_subplot(1,1,1)\n", "\n", "# ice flux uH = ax\n", "qflux = M() * x #* spa m2/s\n", "qflux_scale = qflux * spa*1e-6 #km2/yr\n", "\n", "#thk = thickness(experiment, step, x(mode))[1]\n", "#uvel = qflux / thk #m/s\n", "\n", "ax.plot(xkm, qflux_scale, color='k')\n", "\n", "r = rhoi / rhosw\n", "h_gr = r ** (-1.) * (-topg-sl)\n", "rho_g = rhoi * g()\n", "\n", "# Values of A to be used in MISMIP experiment 3a.\n", "A3a = np.array([3.0e-25, 2.5e-25, 2.0e-25,\n", " 1.5e-25, 1.0e-25, 5.0e-26,\n", " 2.5e-26, 5.0e-26, 1.0e-25,\n", " 1.5e-25, 2.0e-25, 2.5e-25,\n", " 3.0e-25])\n", "\n", "\n", "xi1=900.0\n", "xi2=1370.0\n", "sl=0.0\n", "\n", "for s in range(0,7): \n", " \n", " qschoof = (A3a[s] * (rho_g)**(n()+1.0) * (1.0 - r)**n() / (4.0**n() * Csl) )**(1.0/(m()+1.0)) \n", " qschoof *= h_gr**((m()+n()+3.0)/(m()+1.0)) * spa*1e-6\n", "\n", " #######################################################\n", "\n", " ax.plot(xkm, qschoof,color=\"C\"+str(s),label=str(A3a[s])) #color='C1'\n", " ax.axvline(x_g(sl, xi1, A3a[s])*1e-3,color=\"C\"+str(s),ls='dashed')\n", " ax.axvline(x_g(sl, xi2, A3a[s])*1e-3,color=\"C\"+str(s),ls='dotted')\n", " \n", "ax.set_xlabel('distance from the divide, km')\n", "ax.set_ylabel('ice flux, km2/yr')\n", "_, _, ymin, ymax = ax.axis(xmin=0, xmax=xkm.max())\n", "ax.set_ylim(0,3) #km2/yr\n", "plt.legend()" ] }, { "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.13.13" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }