.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "api/gallery/glaciology/plot_14_mahaffy_shallow_ice_model.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_api_gallery_glaciology_plot_14_mahaffy_shallow_ice_model.py>`
        to download the full example code or to run this example in your browser via JupyterLite.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_api_gallery_glaciology_plot_14_mahaffy_shallow_ice_model.py:


Mahaffy's shallow-ice model (1976)
==================================

M. W. Mahaffy built the first three-dimensional numerical ice-sheet
model, for the Barnes Ice Cap. Because ice sheets are thin, the flux at each
point follows from the local thickness and surface slope, and the thickness
obeys a nonlinear diffusion equation solved on a staggered grid. Its steady
state on a flat bed is Vialov's (1958) profile, flatter at the dome than
Nye's parabola. ``simulate_flowline_glacier`` applies the same scheme to a
mountain glacier, here grown on a bare bed and then shown responding to a
100 m rise of its equilibrium line.

.. GENERATED FROM PYTHON SOURCE LINES 14-59



.. image-sg:: /api/gallery/glaciology/images/sphx_glr_plot_14_mahaffy_shallow_ice_model_001.png
   :alt: Steady ice-sheet profiles, Flowline glacier, Mahaffy: retreat to a new steady state
   :srcset: /api/gallery/glaciology/images/sphx_glr_plot_14_mahaffy_shallow_ice_model_001.png
   :class: sphx-glr-single-img





.. code-block:: Python


    import matplotlib.pyplot as plt
    import numpy as np

    from earthsciencekit.constants import GRAVITY, ICE_DENSITY, SECONDS_PER_YEAR
    from earthsciencekit.glaciology import (
        glacier_response_time,
        nye_plastic_profile,
        simulate_flowline_glacier,
        vialov_dome_height,
        vialov_profile,
    )

    half_width = 500e3
    dome = vialov_dome_height(half_width, 0.3 / SECONDS_PER_YEAR, 2.4e-24, 3, ICE_DENSITY, GRAVITY)
    distance = np.linspace(0, half_width, 501)
    x = np.arange(0, 20e3 + 1, 200.0)
    bed = 3000 - 0.1 * x
    gradient = 0.007 / SECONDS_PER_YEAR
    spin_up = simulate_flowline_glacier(bed, 200.0, 2600.0, gradient, np.linspace(0, 1500, 76) * SECONDS_PER_YEAR)
    steady = spin_up.values[-1, 2:]
    years = np.linspace(0, 500, 251)
    step = simulate_flowline_glacier(bed, 200.0, 2700.0, gradient, years * SECONDS_PER_YEAR, initial_thickness=steady)
    terminus = steady[steady > 1][-1]
    tau = glacier_response_time(steady.max(), gradient * (2600 - bed[np.count_nonzero(steady > 1)])) / SECONDS_PER_YEAR
    area = step.values[:, 1]
    assert area[-1] < area[0] and abs(area[-1] - area[-2]) < 1e-6 * area[0] and terminus > 0

    fig, (top, middle, bottom) = plt.subplots(3, 1, figsize=(8, 10))
    top.plot(distance / 1e3, vialov_profile(distance, half_width, dome, 3), label="Vialov (shallow ice, n = 3)")
    top.plot(distance / 1e3, nye_plastic_profile(half_width - distance, 1e5, ICE_DENSITY, GRAVITY), "--", label="Nye (plastic, 100 kPa)")
    top.set(xlabel="distance from the divide [km]", ylabel="surface [m]", title="Steady ice-sheet profiles")
    top.legend()
    middle.plot(x / 1e3, bed, color="saddlebrown", label="bed")
    middle.plot(x / 1e3, bed + steady, label="steady, ELA 2600 m")
    middle.plot(x / 1e3, bed + step.values[-1, 2:], label="after 500 years, ELA 2700 m")
    middle.set(xlabel="distance along flow [km]", ylabel="elevation [m]", title="Flowline glacier")
    middle.legend()
    bottom.plot(years, step.values[:, 0] / 1e3)
    bottom.axvline(tau, color="0.6", ls=":", label=f"H / ablation = {tau:.0f} years")
    bottom.set(xlabel="years after the step", ylabel="length [km]", title="Mahaffy: retreat to a new steady state")
    bottom.legend()
    fig.tight_layout()

    plt.show()


.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 2.518 seconds)


.. _sphx_glr_download_api_gallery_glaciology_plot_14_mahaffy_shallow_ice_model.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: lite-badge

      .. image:: images/jupyterlite_badge_logo.svg
        :target: ../../../lite/lab/index.html?path=api/gallery/glaciology/plot_14_mahaffy_shallow_ice_model.ipynb
        :alt: Launch JupyterLite
        :width: 150 px

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: plot_14_mahaffy_shallow_ice_model.ipynb <plot_14_mahaffy_shallow_ice_model.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: plot_14_mahaffy_shallow_ice_model.py <plot_14_mahaffy_shallow_ice_model.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: plot_14_mahaffy_shallow_ice_model.zip <plot_14_mahaffy_shallow_ice_model.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_