idealized_sermeq_ablation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6d1abf99-31ba-4727-9e79-256b8abd5b9a",
   "metadata": {},
   "source": [
    "## Frontal ablation for global glacier models from SERMeQ"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad3e5f34-7ae5-4502-8a57-b411f3ad4e38",
   "metadata": {},
   "source": [
    "We would like to develop a frontal ablation function based on the physics of [SERMeQ](https://github.com/ehultee/SERMeQ) that can be used in global glacier models.  The function should become a module directly exchangeable with the common \"k-calving\" parametrization.\n",
    "\n",
    "This notebook is a test setting for developing that functionality.  Core developers are Lizz Ultee & Ruitang Yang.  Test cases are based on previous work by F. Maussion & L. Ultee for [chakra](https://github.com/ehultee/chakra), and on the OGGM [tutorial notebook covering k-calving](https://oggm.org/tutorials/stable/notebooks/kcalving_parameterization.html).\n",
    "\n",
    "15 Feb 2023 | EHU\n",
    "- Edited 16 Feb: flip sign convention and replace test case"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9ba1c5e-c5f1-4ac5-b0e4-eadb31dfd87f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# This allows changes in chakra.py to be automatically re-imported\n",
    "# (this is tricky with OOP though, to be used with care)\n",
    "%load_ext autoreload\n",
    "%autoreload 1\n",
    "%aimport chakra"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bd03f367-85f5-4f1e-bd4b-182fe40f7737",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea33cfe2-5520-4b95-ab4d-35963b6ecbcd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# OGGM imports\n",
    "from oggm.core.flowline import FluxBasedModel\n",
    "from oggm.core.massbalance import ScalarMassBalance\n",
    "from oggm.tests.funcs import bu_tidewater_bed\n",
    "from oggm import cfg\n",
    "cfg.initialize(logging_level='WORKFLOW')\n",
    "cfg.PARAMS['cfl_number'] = 0.01  # less numerical instabilities\n",
    "cfg.PARAMS['use_multiprocessing'] = False"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "997b04d3-56bb-48b0-99ba-a309753ae61e",
   "metadata": {},
   "source": [
    "## Velocity-forced viscoplastic calving"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1dd0c89-47cf-48ec-92c1-97a84360e811",
   "metadata": {},
   "source": [
    "First, we'll define some constants and a preliminary function for yield strength.  Then, we'll write a function that takes profiles and velocity from any model and outputs a viscoplastic frontal ablation rate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d1ae2c49-2628-4b75-81c8-3cb41fd87778",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Global constants\n",
    "G = 9.8 # acceleration due to gravity in m/s^2\n",
    "RHO_ICE = 920.0 #ice density kg/m^3\n",
    "RHO_SEA = 1020.0 #seawater density kg/m^3\n",
    "\n",
    "def tau_y(tau0=150e3, yield_type='constant', bed_elev=None, thick=None, mu=0.01):\n",
    "    \"\"\"Functional form of yield strength.\n",
    "    Can do constant or Mohr-Coulomb yield strength.  Ideally, the glacier's yield type\n",
    "    ('constant' or 'variable') would be saved in a model instance.\n",
    "\n",
    "    Parameters\n",
    "    ----------\n",
    "    tau0: float, optional\n",
    "        Initial guess for yield strength [Pa]. Default is 150 kPa.\n",
    "    yield_type: str, optional\n",
    "        'constant' or 'variable' (Mohr-Coulomb) yielding. Default is constant.\n",
    "    bed_elev: float, optional\n",
    "        Bed elevation, dimensional [m]. The default is None.\n",
    "    thick: float, optional\n",
    "        Ice thickness, dimensional [m]. The default is None.\n",
    "    mu: float, optional\n",
    "        Mohr-Coulomb cohesion, a coefficient between 0 and 1. Default is 0.01.\n",
    "\n",
    "    Returns\n",
    "    -------\n",
    "    tau_y: float\n",
    "        The yield strength for these conditions.\n",
    "    \"\"\"\n",
    "    if yield_type=='variable':\n",
    "        try:\n",
    "            if bed_elev<0:\n",
    "                D = -1*bed_elev #Water depth D the nondim bed topography value when Z<0\n",
    "            else:\n",
    "                D = 0\n",
    "        except:\n",
    "            print('You must set a bed elevation and ice thickness to use variable yield strength.')\n",
    "        N = RHO_ICE*G*thick - RHO_SEA*G*D #Normal stress at bed\n",
    "        ty = tau0 + mu*N\n",
    "    else: #assume constant if not set\n",
    "        ty = tau0\n",
    "    return ty\n",
    "\n",
    "\n",
    "def balance_thickness(yield_strength, bed_elev):\n",
    "    \"\"\"Ice thickness such that the stress matches the yield strength.\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    yield_strength: float\n",
    "        The yield strength near the terminus.  \n",
    "        If yield type is constant, this will of course be the same everywhere.  If yield type is \n",
    "        variable (Mohr-Coulomb), the yield strength at the terminus could differ from elsewhere.\n",
    "    bed_elev: float\n",
    "        Elevation of glacier bed at the terminus\n",
    "        \n",
    "    Returns\n",
    "    -------\n",
    "    Hy: float\n",
    "        The ice thickness for stress balance at the terminus.\n",
    "    \"\"\"\n",
    "    \n",
    "    if bed_elev<0:\n",
    "        D = -1*bed_elev\n",
    "    else:\n",
    "        D = 0\n",
    "    return (2*yield_strength/(RHO_ICE*G)) + np.sqrt((RHO_SEA*(D**2)/RHO_ICE)+(2*yield_strength/(RHO_ICE*G)))\n",
    "## TODO: Check on exponent on last term.  In Ultee & Bassis 2016, this is squared, but in Ultee & Bassis 2020 supplement, it isn't.\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09d44b74-e4ad-413e-971e-f0c386bc78a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def fa_from_velocity(profile, model_velocity, terminus_mb=None, verbose=False,\n",
    "                    tau0=150e3, yield_type='constant', mu=0.01,\n",
    "                    trim_profile=0):\n",
    "    \"\"\"Compute frontal ablation given velocity forcing\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    profile: NDarray\n",
    "        The current profile (x, surface, bed) as calculated by the base model\n",
    "        Unlike core SERMeQ, these should be DIMENSIONAL [m].\n",
    "    model_velocity: array\n",
    "        Velocity along the flowline [m/a] as calculated by the base model \n",
    "        Should have values for the points nearest the terminus...otherwise doesn't matter if this\n",
    "        is the same shape as the profile array.\n",
    "    terminus_mb : float, optional\n",
    "        Mass balance nearest the terminus [m/a]. Default None...TODO: set default behavior\n",
    "    verbose: Boolean, optional\n",
    "        Whether to print component parts for inspection.  Default False.\n",
    "        \n",
    "    tau0: float, optional\n",
    "        This glacier's yield strength [Pa]. Default is 150 kPa.\n",
    "    yield_type: str, optional\n",
    "        'constant' or 'variable' (Mohr-Coulomb) yielding. Default is constant.\n",
    "    mu: float, optional\n",
    "        Mohr-Coulomb cohesion, a coefficient between 0 and 1. Default is 0.01.\n",
    "        Only used if we have variable yield\n",
    "        \n",
    "    trim_profile: int, optional\n",
    "        How many grid cells at the end of the profile to ignore.  Default is 1.\n",
    "        If the initial profile is set by k-calving (as in testing) there can be a \n",
    "        weird cliff shape with very thin final grid point and large velocity gradient\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "    fa_viscoplastic: float\n",
    "        Frontal ablation rate [m/a] based on viscoplastic assumptions\n",
    "    \"\"\"\n",
    "    last_index=-1*(trim_profile+1) ## remove lowest cells if needed\n",
    "        \n",
    "    ## Ice thickness and yield thickness nearest the terminus\n",
    "    se_terminus = profile[1][last_index]\n",
    "    bed_terminus = profile[2][last_index]\n",
    "    H_terminus = se_terminus - bed_terminus\n",
    "    tau_y_terminus = tau_y(tau0=tau0, bed_elev=bed_terminus, thick=H_terminus)\n",
    "    Hy_terminus = balance_thickness(yield_strength=tau_y_terminus, bed_elev=bed_terminus)\n",
    "    U_terminus = model_velocity[last_index] ## velocity, assuming last point is terminus\n",
    "    \n",
    "    ## Ice thickness and yield thickness at adjacent point\n",
    "    se_adj = profile[1][last_index-1]\n",
    "    bed_adj = profile[2][last_index-1]\n",
    "    H_adj = se_adj - bed_adj\n",
    "    tau_y_adj = tau_y(tau0=tau0,bed_elev=bed_adj, thick=H_adj)\n",
    "    Hy_adj = balance_thickness(yield_strength=tau_y_adj, bed_elev=bed_adj)\n",
    "    U_adj = model_velocity[last_index-1]\n",
    "    \n",
    "    \n",
    "    # Gradients\n",
    "    dx_term = abs(profile[0][last_index-1] - profile[0][last_index]) ## check grid spacing close to terminus\n",
    "    dHdx = (H_terminus-H_adj)/dx_term\n",
    "    dHydx = (Hy_terminus-H_adj)/dx_term\n",
    "    dUdx = (U_terminus-U_adj)/ dx_term ## velocity gradient\n",
    "    \n",
    "    \n",
    "    ## Group the terms\n",
    "    dLdt_numerator = terminus_mb - (H_terminus * dUdx) - (U_terminus * dHdx)\n",
    "    dLdt_denominator = dHydx - dHdx ## TODO: compute dHydx\n",
    "    dLdt_viscoplastic = dLdt_numerator / dLdt_denominator\n",
    "    \n",
    "    fa_viscoplastic = U_terminus - dLdt_viscoplastic ## frontal ablation rate\n",
    "    \n",
    "    if verbose:\n",
    "        print('For inspection on debugging - all should be DIMENSIONAL (m/a):')\n",
    "#         print('profile_length={}'.format(profile_length))\n",
    "        print('se_terminus={}'.format(se_terminus))\n",
    "        print('bed_terminus={}'.format(bed_terminus))\n",
    "        print('Thicknesses: Hterm {}, Hadj {}'.format(H_terminus, H_adj))\n",
    "        print('Hy_terminus={}'.format(Hy_terminus))\n",
    "        print('U_terminus={}'.format(U_terminus))\n",
    "        print('dx_term={}'.format(dx_term))\n",
    "        print('Checking dLdt: terminus_mb = {}. \\n H dUdx = {}. \\n U dHdx = {}.'.format(terminus_mb, dUdx*H_terminus, U_terminus*dHdx)) \n",
    "        print('Denom: dHydx = {} \\n dHdx = {}'.format(dHydx, dHdx))\n",
    "        print('Viscoplastic dLdt={}'.format(dLdt_viscoplastic))\n",
    "    else:\n",
    "        pass\n",
    "    \n",
    "    return fa_viscoplastic\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f02f7b64-467f-481c-9468-499eb816058d",
   "metadata": {},
   "source": [
    "### Testing on an idealized profile"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0a65f50-a96d-47bd-92e5-a7e60c862bbc",
   "metadata": {},
   "source": [
    "#### BU bed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4df2ae2-8771-433a-a573-dd88e8caba06",
   "metadata": {},
   "outputs": [],
   "source": [
    "bu_fl = chakra.bu_tidewater_bed()[0]\n",
    "\n",
    "xc = bu_fl.dis_on_line * bu_fl.dx_meter / 1000\n",
    "f, ax = plt.subplots(1, 1, figsize=(12, 5))\n",
    "ax.plot(xc, bu_fl.bed_h, color='k')\n",
    "plt.hlines(0, *xc[[0, -1]], color='C0', linestyles=':')\n",
    "plt.ylim(-350, 1000); plt.ylabel('Altitude [m]'); plt.xlabel('Distance along flowline [km]');"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e02b754-019a-4c45-8f2b-057320031dcb",
   "metadata": {},
   "source": [
    "#### Put a k-calving profile on top of this and check what viscoplastic calving rates we'd get"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4502f2cc-3166-4533-a83e-9bffc1c5ea9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "mb_model = ScalarMassBalance()\n",
    "\n",
    "model = FluxBasedModel(bu_tidewater_bed(), mb_model=mb_model,\n",
    "                       is_tidewater=True, \n",
    "                       calving_use_limiter=True,  # default is True\n",
    "                       flux_gate=0.06,  # default is 0\n",
    "                       calving_k=0.2,  # default is 2.4\n",
    "                       do_kcalving=True\n",
    "                      )\n",
    "# long enough to reach approx. equilibrium \n",
    "ds = model.run_until_and_store(7000)\n",
    "df_diag = model.get_diagnostics()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fdb45b1-365b-48ce-b163-710fad34e1fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "cfg.PARAMS['calving_k']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cace65ad-9f0a-4d61-9f13-b7c15ce73677",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_diag"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "307eff69-04fd-4e04-bc34-328d9fb114d2",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_diag.ice_velocity.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3c1ba2d-a571-429c-ad3b-1875e1b124ca",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "np.argwhere(df_diag.surface_h.values>0)[-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7082df59-2716-4924-b934-a96537cfa64b",
   "metadata": {},
   "outputs": [],
   "source": [
    "scaled_x = 1e-3 * (df_diag.index.values)\n",
    "\n",
    "fig, ax = plt.subplots(1,1, figsize=(12,5))\n",
    "# df_diag.surface_h.plot(ax=ax, color='Gainsboro')\n",
    "# df_diag.bed_h.plot(ax=ax, color='k')\n",
    "ax.plot(scaled_x, df_diag.surface_h, color='Gainsboro')\n",
    "ax.plot(scaled_x, df_diag.bed_h, color='k')\n",
    "ax.axhline(0, color='C0', linestyle=':')\n",
    "plt.ylim(-350, 1000); plt.ylabel('Altitude [m]'); plt.xlabel('Distance along flowline [km]');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "68c2d428-dba7-4acb-a793-cd2175ecfd1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "model.calving_rate_myr"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c596f8d-0d28-4989-ae15-5146ea2dc654",
   "metadata": {},
   "source": [
    "Okay, this clean test case has a sensible calving rate of 18 m/yr, consistent with what's shown in the Oerlemans-Nick k-calving tutorial on the OGGM site."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9ddef5ac-9720-4b57-b465-28751f026812",
   "metadata": {},
   "outputs": [],
   "source": [
    "surface_profile = df_diag.surface_h\n",
    "bed_profile = df_diag.bed_h\n",
    "x = df_diag.index\n",
    "model_U = df_diag.ice_velocity * cfg.SEC_IN_YEAR ## convert ice velocity to m/a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05952086-e826-4f28-ace2-a1d59a233373",
   "metadata": {},
   "outputs": [],
   "source": [
    "max(model_U)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "acf485fd-cdcf-4866-9a70-e977f43a72cc",
   "metadata": {},
   "source": [
    "Okay, the maximum ice velocity is 21 m/a.  This is pretty slow.  Let's see how this goes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa79fe66-4da8-4d99-9b98-a502077ee426",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Find index of the terminus\n",
    "term_index = int(np.argwhere(surface_profile.values>0)[-1])\n",
    "\n",
    "input_profile = (x.values[:term_index+1], ## slice up to index+1 to include the last nonzero value\n",
    "                 surface_profile.values[:term_index+1],\n",
    "                 bed_profile.values[:term_index+1])\n",
    "input_velocity = model_U.values[:term_index+1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5783eebf-857b-4886-bb03-79d80067750f",
   "metadata": {},
   "outputs": [],
   "source": [
    "no_mb.get_annual_mb(heights=surface_profile.values[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52377094-d441-4be3-a4ed-65049b944b34",
   "metadata": {},
   "outputs": [],
   "source": [
    "testval = fa_from_velocity(profile=input_profile, model_velocity=input_velocity, terminus_mb=0, verbose=False,\n",
    "                    tau0=150e3, yield_type='constant', mu=0.01)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d107643c-5c72-4977-ac8a-6af6760e7eaf",
   "metadata": {},
   "outputs": [],
   "source": [
    "testval"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c09df6cc-3d19-4814-a8ab-c1954b0a5b91",
   "metadata": {},
   "source": [
    "Okay, we got a value that seems reasonable!  This is a pretty high ablation rate given the flow speed of the glacier.  Let's explore a bit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "082a4c62-f1f0-4052-938d-28f12bc3f959",
   "metadata": {},
   "outputs": [],
   "source": [
    "ty_tests = np.linspace(25e3, 500e3, num=50)\n",
    "fa_results = [fa_from_velocity(profile=input_profile, model_velocity=input_velocity, terminus_mb=0, verbose=False,\n",
    "                    tau0=ty, yield_type='constant', mu=0.01,trim_profile=1) for ty in ty_tests]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "449ff89c-cda6-44e3-8554-bdeb597d6b5a",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(1e-3*ty_tests, fa_results)\n",
    "ax.set(xlabel='Yield strength [kPa]',\n",
    "       ylabel='Frontal ablation rate')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "629d9c45-b17e-4cfc-9e68-2fba0cdcc301",
   "metadata": {},
   "outputs": [],
   "source": [
    "for ty in ty_tests[0:3]:\n",
    "    t = fa_from_velocity(profile=input_profile, model_velocity=input_velocity, terminus_mb=0, verbose=True,\n",
    "                    tau0=ty, yield_type='constant', mu=0.01, trim_profile=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80dbfffe-f412-4b83-8a9c-3f0652a43ddb",
   "metadata": {},
   "outputs": [],
   "source": [
    "mb_tests = np.linspace(-10, 10, num=50)\n",
    "fa_results_mb = [fa_from_velocity(profile=input_profile, model_velocity=input_velocity, terminus_mb=mb, verbose=False,\n",
    "                    tau0=150e3, yield_type='constant', mu=0.01, trim_profile=1) for mb in mb_tests]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cf634c19-d5ca-4e88-a240-ceb644693aad",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "ax.plot(mb_tests, fa_results_mb)\n",
    "ax.set(xlabel='Terminus mass balance [m/a]',\n",
    "       ylabel='Frontal ablation rate')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b4eb8b71-f7f2-4965-a2ef-f44ca05a6198",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2) = plt.subplots(1,2,sharey=True)\n",
    "ax1.plot(1e-3*ty_tests, fa_results)\n",
    "ax1.set(xlabel='Yield strength [kPa], mb=0',\n",
    "       ylabel='Frontal ablation rate')\n",
    "ax2.plot(mb_tests, fa_results_mb)\n",
    "ax2.set(xlabel=r'Terminus mass balance [m/a], $\\tau_y$=150 KPa')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "adb96f26-f6c7-4505-aa3f-5f19cf4b61f6",
   "metadata": {},
   "source": [
    "We know our glacier is moving slowly.  What does the frontal ablation rate look like if we scale up the velocity?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7e1e0c14-145b-445b-ab88-d0c2bea26240",
   "metadata": {},
   "outputs": [],
   "source": [
    "v_scalings = np.linspace(0.1,20,num=50)\n",
    "fa_results_v = [fa_from_velocity(profile=input_profile, model_velocity=v*np.asarray(input_velocity), terminus_mb=0, verbose=False,\n",
    "                    tau0=150e3, yield_type='constant', mu=0.01,trim_profile=1) for v in v_scalings]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ba04aa5c-49c9-4538-90a7-a2173ebe2286",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2,ax3) = plt.subplots(1,3,sharey=True, figsize=(12,5))\n",
    "ax1.plot(1e-3*ty_tests, fa_results, label='Constant mb=0 m/a, \\n Base velocity' )\n",
    "ax1.set(xlabel='Yield strength [kPa]',\n",
    "       ylabel='Frontal ablation rate')\n",
    "ax2.plot(mb_tests, fa_results_mb, label=r'Constant $\\tau_y$=150 kPa,'+ '\\n Base velocity')\n",
    "ax2.set(xlabel=r'Terminus mass balance [m/a]')\n",
    "ax3.plot(v_scalings, fa_results_v, label=r'Constant $\\tau_y$=150 kPa,'+ '\\n mb=0 m/a')\n",
    "ax3.set(xlabel=r'Velocity scaling [factor]')\n",
    "for ax in (ax1,ax2,ax3):\n",
    "    ax.legend(loc='upper left')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "440c8f7d-b111-445a-af6c-047a23535e09",
   "metadata": {},
   "source": [
    "This case appears to have an instability around $\\tau_y$=450 kPa.  We can investigate that later, but for now let's just trim it so we can see everything else more clearly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c7ea9891-bce6-4e91-9f5f-3b5a0a49e9cb",
   "metadata": {},
   "outputs": [],
   "source": [
    "fa_results_m = np.ma.masked_where(abs(np.asarray(fa_results))>=10000, fa_results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "012958b5-29a4-4bbb-aa39-5406df36d236",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, (ax1, ax2,ax3) = plt.subplots(1,3,sharey=True, figsize=(12,5))\n",
    "ax1.plot(1e-3*ty_tests, fa_results_m, label='Constant mb=0 m/a, \\n Base velocity' )\n",
    "ax1.set(xlabel='Yield strength [kPa]',\n",
    "       ylabel='Frontal ablation rate')\n",
    "ax2.plot(mb_tests, fa_results_mb, label=r'Constant $\\tau_y$=150 kPa,'+ '\\n Base velocity')\n",
    "ax2.set(xlabel=r'Terminus mass balance [m/a]')\n",
    "ax3.plot(v_scalings, fa_results_v, label=r'Constant $\\tau_y$=150 kPa,'+ '\\n mb=0 m/a')\n",
    "ax3.set(xlabel=r'Velocity scaling [factor]')\n",
    "for ax in (ax1,ax2,ax3):\n",
    "    ax.legend(loc='upper left')\n",
    "    ax.axhline(0, color='k', ls=':')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "97b7be5f-7166-4b53-9760-082a9182b1f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "input_profile[1][-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3381e5d8-9097-4fcb-a481-9b4c18ada606",
   "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.8.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}