#!/usr/bin/env python # coding: utf-8 # # # Ice thickness inversion with frontal ablation: a case study # In this notebook we are illustrating how OGGM searches for a frontal calving flux which is compatible with mass conservation, ice thickness estimation based on Glen's flow law, and the calving parameterization of [Oerlemans and Nick (2005)](https://www.cambridge.org/core/journals/annals-of-glaciology/article/minimal-model-of-a-tidewater-glacier/C6B72F547D8C44CDAAAD337E1F2FC97F). # # For more details, see [Recinos et al. (in review)](https://www.the-cryosphere-discuss.net/tc-2018-254/). # ## Set-up # These are the usual OGGM workflow commands. # In[1]: get_ipython().run_line_magic('matplotlib', 'inline') import seaborn as sns sns.set_context('notebook') # In[2]: import numpy as np import pandas as pd import geopandas as gpd import matplotlib.pyplot as plt import oggm from oggm import cfg, utils, workflow, tasks, graphics from oggm.core import inversion # In[3]: cfg.initialize(logging_level='WORKFLOW') cfg.PATHS['working_dir'] = utils.gettempdir(dirname='OGGM-Calving', reset=True) cfg.PARAMS['border'] = 10 # ## Pick a glacier # Here we experiment with Leconte glacier (`RGI60-01.03622`). But you can try with other glaciers as well! # In[4]: gdir = workflow.init_glacier_regions(['RGI60-01.03622'], from_prepro_level=3)[0] # In[5]: graphics.plot_centerlines(gdir, use_flowlines=True); # **Let's first see the results of the inversion without calving:** # In[6]: graphics.plot_inversion(gdir); # ## The basic principles # We use a simple calving law borrowed from [Oerlemans and Nick (2005)](https://www.cambridge.org/core/journals/annals-of-glaciology/article/minimal-model-of-a-tidewater-glacier/C6B72F547D8C44CDAAAD337E1F2FC97F), which relates the frontal calving flux $q_{calving}$ to the frontal ice thickness $H_f$, the water depth $d$ and the terminus width $w$ # # $$q_{calving} = k H_f d w$$ # # With $q_{calving}$ in km$^3$ yr$^{-1}$, $k$ a calibration parameter (default 2.4 yr$^{-1}$) and $d =$ $H_{f}$ - $E_{t}$ ($E_{t}$ being the free board). # # As explained in [Recinos et al. (in review)](https://www.the-cryosphere-discuss.net/tc-2018-254/), ice conservation methods applied to tidewater glaciers *must* take into account this mass-flux at the terminus, otherwise the ice thickness is underestimated. In fact, the default OGGM ice thickness inversion procedure assumes an ice flux of zero at the terminus: # In[8]: # Default calving flux: inversion.calving_flux_from_depth(gdir) # Here, `flux` is the calving flux ($q_{calving}$) in km$^3$ yr$^{-1}$, `free_board` is the height of ice above water ($E_{t}$, i.e. above 0 m a.s.l), `thick` the frontal ice thickness ($H_{f}$) in m (equal to zero per construction in OGGM), `water_depth` the water depth in m (negative here because of the zero ice thickness and a terminus elevation equal to the freeboard imply a positive bedrock elevation), `width` the front width in m. # # Now let's see how this calving flux would change if we increase the ice thickness (while keeping the free board fixed, because this is the only thing we know "for sure": the surface elevation at the terminus.): # In[10]: df = [] for thick in np.linspace(0, 500, 51): # This function simply computes the calving law df.append(inversion.calving_flux_from_depth(gdir, thick=thick)) df = pd.DataFrame(df).set_index('thick') # In[11]: df[['flux']].plot(); plt.ylabel('Calving flux (km$^3$ yr$^{-1}$)'); # The flux is zero as long as the ice isn't thick enough to reach water, after which the water depth is positive and calving occurs. The calving flux varies with $H$ as a polynomial of degree 2. # # We don't know which is the **real** value of the calving flux at this glacier. From here, let's make some very coarse assumptions: # - the Oerlemans and Nick calving law is perfectly exact # - the tuning parameter $k$ is known # - our glacier is in equilibrium (a fundamental assumption necessary for mass-conservation inversion) # - ice deformation at the glacier terminus follows Glen's flow law # # Under these assumptions, we are now going to show, that there is **one and only one** value for the frontal thickness which complies with *both* the calving law and the ice thickness inversion model of OGGM. # # To test this hypothesis, we are going to compute a calving flux (from the calving law) for a range of frontal ice thicknesses, then give this flux back to the OGGM inversion model, which will use this flux to compute the frontal ice thickness according to the physics of ice flow (see the [documentation](https://docs.oggm.org/en/latest/inversion.html) or [Recinos et al. (in review)](https://www.the-cryosphere-discuss.net/tc-2018-254/) for more info). Let's see how this works: # In[12]: df = [] for thick in np.linspace(0, 500, 51): # This function simply computes the calving law out = inversion.calving_flux_from_depth(gdir, thick=thick) out['thick (prescribed)'] = out.pop('thick') # Now we feed it back to OGGM gdir.inversion_calving_rate = out['flux'] # The mass-balance has to adapt in order to create a flux tasks.local_t_star(gdir) tasks.mu_star_calibration(gdir) tasks.prepare_for_inversion(gdir) v_inv, _ = tasks.mass_conservation_inversion(gdir) # Now we get the OGGM ice thickness out['thick (oggm)'] = inversion.calving_flux_from_depth(gdir)['thick'] # Add sliding (the fs value is outdated, but still) v_inv, _ = tasks.mass_conservation_inversion(gdir, fs=5.7e-20) out['thick (oggm+sliding)'] = inversion.calving_flux_from_depth(gdir)['thick'] # Store df.append(out) df = pd.DataFrame(df) # We've got: # In[13]: df[['flux', 'thick (prescribed)', 'thick (oggm)']].plot(x='flux'); plt.xlabel('Calving flux (km$^3$ yr$^{-1}$)'); plt.ylabel('Ice thickness (m)'); # We already know that the calving law relates the thickness to the flux with a root of degree two (blue curve). Now, what explains the shape of **the orange curve**?. It is Glen's flow law, which relates ice thickness to the flux with a 5th degree root (assuming n=3). This is the reason why there is one (and only one) non-zero solution to the problem of finding a calving flux which is compatible with both the calving law and the physics of ice deformation (under our simplified framework of course). # # Note that adding sliding doesn't change the problem (we still solve a polynome of degree 5 in OGGM, [with a new term in degree 3](https://docs.oggm.org/en/latest/ice-dynamics.html)): # In[14]: df[['flux', 'thick (prescribed)', 'thick (oggm)', 'thick (oggm+sliding)']].plot(x='flux'); plt.xlabel('Calving flux (km$^3$ yr$^{-1}$)'); plt.ylabel('Ice thickness (m)'); # ## Finding the optimum frontal thickness # There are several ways to find this "optimal" thickness, where mass-conservation inversion and the calving law are compatible. In OGGM, we implement an iterative procedure converging to this value in a few iterations: # In[15]: df = inversion.find_inversion_calving(gdir) # In[16]: f, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5)) df[['calving_flux']].plot(ax=ax1); ax1.set_ylabel('Flux (km$^3$ yr$^{-1}$)'); df[['mu_star']].plot(ax=ax2, color='C1'); ax2.set_ylabel('$\mu^*$ (mm mth$^{-1}$ K$^{-1}$)'); # The procedure starts with an arbitrary initial water depth (1/3 the terminus elevation, or at least 10 m), then iteratively feeds the calving flux back to the mass-conservation inversion function of OGGM, which adapts the mass-balance model to cope with a non-zero ice flux at the front. In order to do that, the temperature sensitivity $\mu^*$ has to be reduced (per construction, the original $\mu^*$ is defined such that the flux at the front is zero). # # Convergence is reached when the OGGM flux equals the calving law flux within 0.1%: # In[17]: print('At convergence:') print(df.iloc[-1]) # ### Testing consistency # We will know put the final thickness estimated in the calving law to demontrate we obtained the same flux. # In[22]: thick = df['thick'][10] # In[23]: out = inversion.calving_flux_from_depth(gdir, thick=thick) out # ### Sensitivity to the initial guess # The method always converges, regardless of the starting water depth: # In[16]: df = pd.DataFrame() df['wd_001'] = utils.find_inversion_calving(gdir, water_depth=1)['calving_flux'] df['wd_100'] = utils.find_inversion_calving(gdir, water_depth=100)['calving_flux'] df['wd_200'] = utils.find_inversion_calving(gdir, water_depth=200)['calving_flux'] df['wd_300'] = utils.find_inversion_calving(gdir, water_depth=300)['calving_flux'] df['wd_400'] = utils.find_inversion_calving(gdir, water_depth=400)['calving_flux'] df['wd_500'] = utils.find_inversion_calving(gdir, water_depth=500)['calving_flux'] # In[17]: df.iloc[1:].plot(); # This converges because the solution is unique. # ### What if it "overshoots"? # For some glaciers, the calving flux given by the calving law is larger than can be explained by climate alone, i.e. even without melt, the computed flux is larger than total accumulation over our glacier. This can happen for several reasons: # - precipitation is underestimated # - the flux is overestimated because of uncertainties in k and the terminus geometry # - the equilibrium assumption is not valid # - ... # # Let's simulate this case on our glacier, here by setting an unrealistically large calving $k$: # In[18]: cfg.PARAMS['k_calving'] = 10 # default is 2.4 # In[19]: df = utils.find_inversion_calving(gdir) # In[20]: f, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5)) df[['calving_flux']].plot(ax=ax1); ax1.set_ylabel('Flux (km$^3$ yr$^{-1}$)'); df[['mu_star']].plot(ax=ax2, color='C1'); ax2.set_ylabel('$\mu^*$ (mm mth$^{-1}$ K$^{-1}$)'); # If this happens during the iteration, OGGM is going to set $\mu^*$ to zero and computes the corresponding flux (the maximal physically possible value). In this case, the calving law and OGGM disagree: # In[21]: print('At iteration stop:') print(df.iloc[-1]) # ## Discussion # We have shown that under certain assumptions, our method converges to an ice thickness consistent with the calving law. Of course, there are many uncertainties involved, and it is not our attempt to find the "true" calving flux here. See [Recinos et al. (in review)](https://www.the-cryosphere-discuss.net/tc-2018-254/) for our motivations and a discussion about why it is crucial to include frontal ablation in regional ice thickness estimates. # ## What's next? # You have several options from here: # - return to the [OGGM documentation](https://docs.oggm.org) # - explore other tutorials on the [OGGM-Edu](https://edu.oggm.org) platform.