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).

For more details, see Recinos et al., (2019).

History:

Feb 2019: this file was created for the first round of reviews. This original version is still available here
Jun 2019: the file was moved to a new repository and modified following a new methodology as suggested by the reviewers in the second round of reviews
April 2020: updated for recent changes in code

Set-up#

These are the usual OGGM workflow commands.

import seaborn as sns
sns.set_context('notebook')

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.inversion import find_inversion_calving, calving_flux_from_depth

cfg.initialize()
cfg.PATHS['working_dir'] = utils.gettempdir(dirname='OGGM-Calving', reset=True)
cfg.PARAMS['border'] = 10

2023-03-07 12:35:24: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2023-03-07 12:35:24: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2023-03-07 12:35:24: oggm.cfg: Multiprocessing: using all available processors (N=2)
2023-03-07 12:35:24: oggm.cfg: PARAMS['border'] changed from `40` to `10`.

Pick a glacier#

Here we experiment with Leconte glacier (RGI60-01.03622). But you can try with other glaciers as well! See this list for nearly all glaciers simulated in Recinos et al., (2019) (we took out the glaciers that need special processing, i.e. because of incorrect outlines like Columbia).

gdir = workflow.init_glacier_directories(['RGI60-01.03622'], from_prepro_level=3)[0]

2023-03-07 12:35:24: oggm.workflow: init_glacier_directories from prepro level 3 on 1 glaciers.
2023-03-07 12:35:24: oggm.workflow: Execute entity tasks [gdir_from_prepro] on 1 glaciers
2023-03-07 12:35:24: oggm.utils: Downloading https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.4/L3-L5_files/CRU/centerlines/qc3/pcp2.5/no_match/RGI62/b_010/L3/RGI60-01/RGI60-01.03.tar to /github/home/OGGM/download_cache/cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.4/L3-L5_files/CRU/centerlines/qc3/pcp2.5/no_match/RGI62/b_010/L3/RGI60-01/RGI60-01.03.tar...
2023-03-07 12:35:29: oggm.utils: /github/home/OGGM/download_cache/cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.4/L3-L5_files/CRU/centerlines/qc3/pcp2.5/no_match/RGI62/b_010/L3/RGI60-01/RGI60-01.03.tar verified successfully.

graphics.plot_centerlines(gdir, use_flowlines=True);

/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
  out = pickle.load(f)

../_images/inversion_with_frontal_ablation_10_1.png

Let’s first see the results of the inversion without calving:

graphics.plot_inversion(gdir);

/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
  out = pickle.load(f)

../_images/inversion_with_frontal_ablation_12_1.png

The basic principles#

We use a simple calving law borrowed from Oerlemans and Nick (2005), 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 \qquad \qquad \qquad (eq. \: 1)\]

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., (2019), 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:

# Default calving flux:
calving_flux_from_depth(gdir, water_level=0)

{'flux': 0.001618211699679134,
 'width': 775.0815795377595,
 'thick': 113.79073958092376,
 'inversion_calving_k': 0.6,
 'water_depth': 30.579455682270137,
 'water_level': 0,
 'free_board': 83.21128389865362}

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 the surface elevation at the terminus is the only thing we know with “certainty”):

df = []
for thick in np.linspace(0, 500, 51):
    # This function simply computes the calving law
    df.append(calving_flux_from_depth(gdir, thick=thick, water_level=0))
df = pd.DataFrame(df).set_index('thick')

df[['flux']].plot();
plt.ylabel('Calving flux (km$^3$ yr$^{-1}$)');

../_images/inversion_with_frontal_ablation_18_0.png

The flux is zero as long as the ice isn’t thick enough to reach the water surface, 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 are two solutions for the frontal thickness which complies with the calving law and the ice thickness inversion model of OGGM, but only one solution that provides a realistic thickness.

To make our point, 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 or Recinos et al., (2019) for more info).

Let’s see how this works:

# Reference data for the recalibration of the mass-balance
ref_tstars_url = 'https://cluster.klima.uni-bremen.de/~oggm/ref_mb_params/oggm_v1.4/RGIV62/CRU/centerlines/qc3/pcp2.5'
workflow.download_ref_tstars(base_url=ref_tstars_url)

2023-03-07 12:35:37: oggm.utils: Downloading https://cluster.klima.uni-bremen.de/~oggm/ref_mb_params/oggm_v1.4/RGIV62/CRU/centerlines/qc3/pcp2.5/ref_tstars.csv to /github/home/OGGM/download_cache/cluster.klima.uni-bremen.de/~oggm/ref_mb_params/oggm_v1.4/RGIV62/CRU/centerlines/qc3/pcp2.5/ref_tstars.csv...
2023-03-07 12:35:37: oggm.utils: /github/home/OGGM/download_cache/cluster.klima.uni-bremen.de/~oggm/ref_mb_params/oggm_v1.4/RGIV62/CRU/centerlines/qc3/pcp2.5/ref_tstars.csv verified successfully.
2023-03-07 12:35:37: oggm.utils: Downloading https://cluster.klima.uni-bremen.de/~oggm/ref_mb_params/oggm_v1.4/RGIV62/CRU/centerlines/qc3/pcp2.5/ref_tstars_params.json to /github/home/OGGM/download_cache/cluster.klima.uni-bremen.de/~oggm/ref_mb_params/oggm_v1.4/RGIV62/CRU/centerlines/qc3/pcp2.5/ref_tstars_params.json...
2023-03-07 12:35:38: oggm.utils: /github/home/OGGM/download_cache/cluster.klima.uni-bremen.de/~oggm/ref_mb_params/oggm_v1.4/RGIV62/CRU/centerlines/qc3/pcp2.5/ref_tstars_params.json verified successfully.

with utils.DisableLogger():  # this scraps some output - to use with caution!!!
    # Ouput dataframe
    df = []
    # Range of frontal ice thickness from 0 to 500m
    for thick in np.linspace(0, 500, 51):

        # This function simply computes the flux from the calving law
        out = calving_flux_from_depth(gdir, thick=thick, water_level=0)
        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)'] = calving_flux_from_depth(gdir, water_level=0)['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)'] = calving_flux_from_depth(gdir, water_level=0)['thick']

        # Store
        df.append(out)

df = pd.DataFrame(df)

/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
  out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
  out = pickle.load(f)

We’ve got:

df[['flux', 'thick (prescribed)', 'thick (oggm)']].plot(x='flux');
plt.xlabel('Calving flux (km$^3$ yr$^{-1}$)');
plt.ylabel('Ice thickness (m)');

../_images/inversion_with_frontal_ablation_23_0.png

We already know that the calving law (eq. 1) 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). The graphic shows that, for this glacier, there is at least one solution to the problem of finding a calving flux which is compatible with both the calving law and the physics of ice deformation (where the two lines cross).

Note that adding sliding doesn’t change the problem (we still solve a polynomial of degree 5 in OGGM, with a new term in degree 3):

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)');   

../_images/inversion_with_frontal_ablation_25_0.png

This plot is a bit misleading. The orange and green curve cross the blue line in two places, as we are going to show.

Finding the “correct” flux#

What we are looking for are the frontal ice thickness (and corresponding calving flux) that satisfy the following equation:

\[q_{calving}(h) = q_{OGGM}(h) \qquad \qquad \qquad (eq. \: 2)\]

Where \(q_{calving}\) is the flux of ice from the calving law and \(q_{OGGM}\) is the flux of ice due to deformation. See Recinos et al., (2019) for the details on the equations.

In this case, eq. 2 is a polynomial and could therefore be solved analytically and with root finding methods. However, we prefer to solve it numerically, for two main reasons:

finding a numerical solution implied the least change (and copy-paste) in the existing and well-tested code base, therefore minimizing the risk of creating bugs
numerical solvers are very flexible and can virtually be applied to any other formulation of \(q_{calving}\) and \(q_{OGGM}\), i.e. the method will still be applicable if we use the lateral drag parameterization in OGGM or use another formulation for the calving law.

Now solve eq. 2 with a Bound-Constrained minimization method. We solve for water depth \(d\) (not \(h\)), because the two are equivalent if the free-board \(E_{t}\) is known (\(h\) = \(d\) + \(E_{t}\)).

Let’s have a look at the function(s) we are seeking roots for:

from oggm.core.inversion import sia_thickness

# Slope and with at the glacier front
cls = gdir.read_pickle('inversion_input')[-1]
slope = cls['slope_angle'][-1]
width = cls['width'][-1]

def to_solve_default(wd):
    # Flux from calving law at given water depth
    fl = calving_flux_from_depth(gdir, water_depth=wd, water_level=0)
    # OGGM thickness obtained from the prescribed flux
    oggm = sia_thickness(slope, width, fl['flux'] * 1e9 / cfg.SEC_IN_YEAR)
    # Difference in thicknesses
    return fl['thick'] - oggm


def to_solve_with_sliding(wd):
    # Flux from calving law at given water depth
    fl = calving_flux_from_depth(gdir, water_depth=wd, water_level=0)
    # OGGM thickness obtained from the prescribed flux
    oggm = sia_thickness(slope, width, fl['flux'] * 1e9 / cfg.SEC_IN_YEAR, fs=5.7e-20)
    # Difference in thicknesses
    return fl['thick'] - oggm

And plot them for a range of possible water depths:

# Ouput dataframe
df = []

# Range of frontal water depth from 0.1 to 400m
for wd in np.linspace(0.1, 400, 100):
    
    out = {
        'Water depth [m]': wd,
        'Calving law - OGGM': to_solve_default(wd),
        'Calving law - OGGM (with sliding)': to_solve_with_sliding(wd),
          }
    
    # Store
    df.append(out)
    
df = pd.DataFrame(df).set_index('Water depth [m]')
df.plot(color=['C1', 'C2']);
plt.hlines([0], 0, 400);
plt.xlabel('Water depth [m]');
plt.ylabel('Difference in thickness [m]');
plt.legend(loc='upper right');

../_images/inversion_with_frontal_ablation_34_0.png

In this representation of the same data as above (here as a function of the prescribed water depth), one sees more clearly that there exist two locations where the zero line is crossed. This pattern is consistent for all but one glaciers in Alaska (see below), and we choose to pick the only realistic solution (the larger of the two). Indeed, the smallest solution leads to unrealistically small water depths and calving fluxes. Note that in the first two versions of the manuscript (with a different way to solve the equations), it is always the largest solution that was found: therefore, there is no noticeable change in the results between the last two manuscript versions.

To solve our problem, we first seek the absolute minimum of the function(s), and then look for zeros above and below this point. It’s the larger solution which is retained.

from scipy import optimize
# Absolute minima
abs_min = optimize.minimize(to_solve_default, [1], bounds=((1e-4, 1e4), ))
abs_min['x'][0]

42.147279259968215

Now we find the roots of the functions between the interval [0.0001 - 10000]. Using the python tool: optimize.brentq().

As observed in the previous plot we have two minima:

print('Optimum water depth 1 (in m)', optimize.brentq(to_solve_default, 1e-4, abs_min['x'][0]))
print('Optimum water depth 2 (in m)', optimize.brentq(to_solve_default, abs_min['x'][0], 1e4))

Optimum water depth 1 (in m) 2.9899016853556635
Optimum water depth 2 (in m) 143.7924777099962

A frontal ablation flux produced with such small water depth will result in a non realistic thickness at the front. In the code we always pick the larger solution.

opt = optimize.brentq(to_solve_default, abs_min['x'][0], 1e4)

Now we give this value to the calving law and calculate a flux:

out = calving_flux_from_depth(gdir, water_depth=opt, water_level=0)
out

{'flux': 0.015179864222832417,
 'width': 775.0815795377595,
 'thick': 227.00376160864982,
 'inversion_calving_k': 0.6,
 'water_depth': 143.7924777099962,
 'water_level': 0,
 'free_board': 83.21128389865362}

The process described above is implemented in the OGGM codebase in the following task: find_inversion_calving(). This is what we use in operation, when use_kcalving_for_inversion is switched on, and it provides the same results:

# this is switched off in the default settings
cfg.PARAMS['use_kcalving_for_inversion'] = True
out = find_inversion_calving(gdir, water_level=0)
out

2023-03-07 12:36:12: oggm.cfg: PARAMS['use_kcalving_for_inversion'] changed from `False` to `True`.
2023-03-07 12:36:12: oggm.core.inversion: (RGI60-01.03622) find_inversion_calving
2023-03-07 12:36:13: oggm.core.inversion: (RGI60-01.03622) find_inversion_calving_from_any_mb: found calving flux of 0.015 km3 yr-1

{'calving_flux': 0.015179864222832573,
 'calving_rate_myr': 86.27548662599843,
 'calving_mu_star': 91.70437425591591,
 'calving_law_flux': 0.015179864222832498,
 'calving_water_level': 0,
 'calving_inversion_k': 0.6,
 'calving_front_slope': 0.11600123219190657,
 'calving_front_water_depth': 143.7924777099967,
 'calving_front_free_board': 83.21128389865362,
 'calving_front_thick': 227.0037616086503,
 'calving_front_width': 775.0815795377595}

Testing consistency#

Per definition, the two calving fluxes (OGGM and calving law) are now the same. It is a good consistency check to verify that this is the case:

thick = out['calving_front_thick'] 
df = calving_flux_from_depth(gdir, thick=thick, water_level=0)
np.testing.assert_almost_equal(df['flux'], out['calving_flux']) 
np.testing.assert_almost_equal(gdir.read_pickle('inversion_output')[-1]['thick'][-1], thick)

Same flux!!

Free board and water level#

In the example above, the glacier has a “free-board” (thickness of the calving front above water) of about 85 m. This value is given by the glacier DEM (by construction of the OGGM flowlines, it is the lowest point on the glacier outline). Sometimes, this value can be unrealistically high (> 200 m a.s.l) or low (0 m a.s.l), leading to unrealistic calving values.

Since this change, we are now artificially adapting the water level so that the free-board is within a reasonable range. This range is given by the model parameters:

cfg.PARAMS['free_board_marine_terminating']  # units: m

[10.0, 50.0]

This means that we do not allow free-board higher than 50 (conservative value, could be changed to 75 m). This has an influence on the results. See how the calving flux changes with these new defaults:

find_inversion_calving(gdir, water_level=None)  # use default values

2023-03-07 12:36:14: oggm.core.inversion: (RGI60-01.03622) find_inversion_calving
2023-03-07 12:36:15: oggm.core.inversion: (RGI60-01.03622) find_inversion_calving_from_any_mb: found calving flux of 0.022 km3 yr-1

{'calving_flux': 0.02220788362949491,
 'calving_rate_myr': 116.97117335284842,
 'calving_mu_star': 91.28106176613524,
 'calving_law_flux': 0.0222078836294945,
 'calving_water_level': 33.21128389865362,
 'calving_inversion_k': 0.6,
 'calving_front_slope': 0.11600123219190657,
 'calving_front_water_depth': 194.9519555880771,
 'calving_front_free_board': 50.0,
 'calving_front_thick': 244.9519555880771,
 'calving_front_width': 775.0815795377595}

In this case, the calving flux slightly increased as a result of an increase in the water depth and the frontal thickness.

Mass conservation and the temperature sensitivity \(\mu^*\)#

While the computed flux satisfies eq. 2, we still need to check with the mass balace module if the flux computed above complies with mass conservation. Per construction, OGGM calibrates the mass-balance model so that the flux at the glacier front is zero under equilibrium climate conditions. This is obviously not the case for calving glaciers, and this newly computed flux can be used to better calibrate the model:

# Feed a zero flux to OGGM
gdir.inversion_calving_rate = 0

# Calibrate the mass balance
tasks.local_t_star(gdir)    
tasks.mu_star_calibration(gdir)

# Output
print('mu* for the glacier without calving (in mm w.e. mth-1): '
      '{:.2f}'.format(gdir.read_json('local_mustar')['mu_star_glacierwide']))

# Feed the computed flux now
gdir.inversion_calving_rate = out['calving_flux']

# Calibrate the mass balance
tasks.local_t_star(gdir)    
tasks.mu_star_calibration(gdir)

# Output
print('mu* for the glacier with calving (in mm w.e. mth-1): '
      '{:.2f}'.format(gdir.read_json('local_mustar')['mu_star_glacierwide']))

# Invert for fun
tasks.prepare_for_inversion(gdir)
v_inv = tasks.mass_conservation_inversion(gdir)

mu* for the glacier without calving (in mm w.e. mth-1): 92.62
mu* for the glacier with calving (in mm w.e. mth-1): 91.70

2023-03-07 12:36:15: oggm.core.climate: (RGI60-01.03622) local_t_star
2023-03-07 12:36:15: oggm.core.climate: (RGI60-01.03622) local mu* computation for t*=1964
2023-03-07 12:36:15: oggm.core.climate: (RGI60-01.03622) mu_star_calibration
2023-03-07 12:36:16: oggm.core.climate: (RGI60-01.03622) local_t_star
2023-03-07 12:36:16: oggm.core.climate: (RGI60-01.03622) local mu* computation for t*=1964
2023-03-07 12:36:16: oggm.core.climate: (RGI60-01.03622) mu_star_calibration
2023-03-07 12:36:16: oggm.core.inversion: (RGI60-01.03622) prepare_for_inversion
2023-03-07 12:36:16: oggm.core.inversion: (RGI60-01.03622) mass_conservation_inversion

This is consistent: in order to allow for a frontal ablation flux, melt has to be reduced along the glacier. Note that not only the frontal thickness is affected by this change, the entire glacier becomes thicker:

graphics.plot_inversion(gdir);

/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
  out = pickle.load(f)

../_images/inversion_with_frontal_ablation_60_1.png

Cases where the calving flux isn’t consistent with climate#

There is a physical lower bound for \(\mu^*\) of zero (which means no melt at all). In this case, all the mass accumulated on the glacier via solid precipitation is transported and lost at the calving front. In our experiments, this happens when the calving flux is overestimated (likely given all the uncertainties involved in the calving law and the boundary conditions (see Recinos et al., (2019)) or the climate data is biased (also likely). In these cases, we don’t accept the solution provided by the calving law but clip the temperature sensitivity \(\mu^*\) to zero, compute the associated frontal flux, and keep the later as the final solution.

We try a big number for the \(k\) parameter to illustrate this event:

# Will provide a flux too large for OGGM to cope with
cfg.PARAMS['inversion_calving_k'] = 2.4 * 10 
# Change the default clip mechanism to zero
cfg.PARAMS['calving_min_mu_star_frac'] = 0

out = find_inversion_calving(gdir)
out

2023-03-07 12:36:21: oggm.cfg: PARAMS['inversion_calving_k'] changed from `0.6` to `24.0`.
2023-03-07 12:36:21: oggm.cfg: PARAMS['calving_min_mu_star_frac'] changed from `0.7` to `0`.
2023-03-07 12:36:21: oggm.core.inversion: (RGI60-01.03622) find_inversion_calving
2023-03-07 12:36:22: oggm.core.inversion: (RGI60-01.03622) find_inversion_calving_from_any_mb: found calving flux of 1.538 km3 yr-1

{'calving_flux': 1.5376960608989283,
 'calving_rate_myr': 3470.3198902418503,
 'calving_mu_star': 0.0,
 'calving_law_flux': 5.547752204120104,
 'calving_water_level': 33.21128389865362,
 'calving_inversion_k': 24.0,
 'calving_front_slope': 0.11600123219190657,
 'calving_front_water_depth': 521.6807607403077,
 'calving_front_free_board': 50.0,
 'calving_front_thick': 571.6807607403077,
 'calving_front_width': 775.0815795377595}

After clipping \(\mu^*\), we keep the thickness and flux that still statisfy mass-conservation and ice deformation (\(q_{OGGM}\)) instead of the calving law (\(q_{calving}\)).

New as for version 1.4: we now have a mechanism to clip mu to a value higher than 0, given as a fraction of the non-calving mu:**

# Will provide a flux too large for OGGM to cope with
cfg.PARAMS['inversion_calving_k'] = 2.4 * 10 
# Change the default clip mechanism to 70 of non-calving mu*
cfg.PARAMS['calving_min_mu_star_frac'] = 0.7

out = find_inversion_calving(gdir)
out

2023-03-07 12:36:22: oggm.cfg: PARAMS['calving_min_mu_star_frac'] changed from `0` to `0.7`.
2023-03-07 12:36:22: oggm.core.inversion: (RGI60-01.03622) find_inversion_calving
2023-03-07 12:36:23: oggm.core.inversion: (RGI60-01.03622) find_inversion_calving_from_any_mb: found calving flux of 0.930 km3 yr-1

{'calving_flux': 0.9300311740488388,
 'calving_rate_myr': 2320.978475131525,
 'calving_mu_star': 36.60094278638051,
 'calving_law_flux': 4.4909861092712395,
 'calving_water_level': 33.21128389865362,
 'calving_inversion_k': 24.0,
 'calving_front_slope': 0.11600123219190657,
 'calving_front_water_depth': 466.98623222645926,
 'calving_front_free_board': 50.0,
 'calving_front_thick': 516.9862322264593,
 'calving_front_width': 775.0815795377595}

Check the robustness of our method#

We apply our algorithm to a handful of glaciers and check the solutions. This is not a formal test, see the original paper for regional results:

cfg.initialize(logging_level='CRITICAL')
cfg.PATHS['working_dir'] = utils.gettempdir(dirname='OGGM-Calving', reset=True)
cfg.PARAMS['border'] = 10

cfg.PARAMS['inversion_calving_k'] = 2.4  # We put it higher for this example

# Randomly selected IDs
ids = [
'RGI60-01.04375',
'RGI60-01.03377',
'RGI60-01.01390',
'RGI60-01.14391',
'RGI60-01.20891',
'RGI60-01.10612',
'RGI60-01.17807'
]

for rid in ids:

    gdir = workflow.init_glacier_directories([rid], from_prepro_level=3)[0]

    cls = gdir.read_pickle('inversion_input')[-1]
    slope = cls['slope_angle'][-1]
    width = cls['width'][-1]

    def to_minimize(wd):
        fl = calving_flux_from_depth(gdir, water_depth=wd, water_level=0)
        oggm = sia_thickness(slope, width, np.array([fl['flux'] * 1e9 / cfg.SEC_IN_YEAR]))
        return fl['thick'] - oggm
    
    abs_min = optimize.minimize(to_minimize, [1], bounds=((1e-4, 1e4), ))
    print(rid)
    print('Optimum small (in m)', optimize.brentq(to_minimize, 1e-4, abs_min['x'][0]))
    print('Optimum large (in m)', optimize.brentq(to_minimize, abs_min['x'][0], 1e4))
# this can take a few minutes

RGI60-01.04375
Optimum small (in m) 3.8791150180855616
Optimum large (in m) 766.7562471578146
RGI60-01.03377
Optimum small (in m) 0.8767092368563039
Optimum large (in m) 525.2806591959095
RGI60-01.01390
Optimum small (in m) 0.12725555240651668
Optimum large (in m) 521.2668860494764
RGI60-01.14391
Optimum small (in m) 0.2399488654573076
Optimum large (in m) 921.3178100109675
RGI60-01.20891
Optimum small (in m) 0.08746271387787898
Optimum large (in m) 656.9118156598653
RGI60-01.10612
Optimum small (in m) 2.3932170290279915
Optimum large (in m) 306.7325531822173
RGI60-01.17807
Optimum small (in m) 7.371858617233925
Optimum large (in m) 130.69421545268608

/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
  out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
  out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
  out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
  out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
  out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
  out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
  out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
  out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
  out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
  out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
  out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
  out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
  out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
  out = pickle.load(f)

The algorithm implemented in OGGM always seeks for the larger solution.

What about glaciers where no optimal thickness exists?#

For 6 glaciers in Alaska, there is no solution which complies with the calving law and ice deformation with the previous default parameters. Here for example glacier RGI60-01.23642 (Tsaa Glacier):

# Reference data for the recalibration of the mass-balance
ref_tstars_url = 'https://cluster.klima.uni-bremen.de/~oggm/ref_mb_params/oggm_v1.4/RGIV62/CRU/centerlines/qc3/pcp2.5'
workflow.download_ref_tstars(base_url=ref_tstars_url)

gdir = workflow.init_glacier_directories(['RGI60-01.23642'], from_prepro_level=3)[0]

cls = gdir.read_pickle('inversion_input')[-1]
slope = cls['slope_angle'][-1]
width = cls['width'][-1]

def to_minimize(wd):
    fl = calving_flux_from_depth(gdir, water_depth=wd, water_level=0)
    oggm = sia_thickness(slope, width, np.array([fl['flux'] * 1e9 / cfg.SEC_IN_YEAR]))
    return fl['thick'] - oggm

wd = np.linspace(0.1, 400)
out = []
for w in wd:
    out.append(to_minimize(w)) 
    
plt.plot(wd, out);
plt.hlines([0], 0, 400);
plt.xlabel('Water depth');
plt.ylabel('Difference Calving Law - OGGM');

/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
  out = pickle.load(f)
/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
  out = pickle.load(f)

../_images/inversion_with_frontal_ablation_73_1.png

Put another way, the two laws are simply incompatible. Both curves will never meet in a domain between \(q_{calving}\) \(>\) 0 and a \(q_{calving}\) that results in a non negative \(\mu*\). This could be due to several reasons:

a difference > 0 between both fluxes in eq. 2 corresponds to a glacier state outside of the equilibrium condition under this specific geometry
the calving law and the ice deformation equation do not represent the physical process affecting this glacier in particulary. Both are parameterised and therefore an approximation of reality
uncertain boundary conditions

Note that the number of cases where this happens is reduced to only one case when we correct for observed water depths and frontal width, indicating that uncertainties in the boundary conditions are the main problem here.

With the new defaults, let’s see if OGGM can find a suitable solution to this glacier:

cfg.PARAMS['inversion_calving_k'] = 0.6
cfg.PARAMS['use_kcalving_for_inversion'] = True
out = tasks.find_inversion_calving(gdir)
out

/usr/local/pyenv/versions/3.10.10/lib/python3.10/site-packages/oggm/utils/_workflow.py:2939: UserWarning: Unpickling a shapely <2.0 geometry object. Please save the pickle again; shapely 2.1 will not have this compatibility.
  out = pickle.load(f)

{'calving_flux': 0.02810564881034136,
 'calving_rate_myr': 57.953483737361545,
 'calving_mu_star': 87.96935557164645,
 'calving_law_flux': 0.028105648810341295,
 'calving_water_level': 83.39508240583478,
 'calving_inversion_k': 0.6,
 'calving_front_slope': 0.18201268515440294,
 'calving_front_water_depth': 96.58913956226903,
 'calving_front_free_board': 50.0,
 'calving_front_thick': 146.58913956226903,
 'calving_front_width': 3308.356272260822}

The glacier now calves.

Discussion#

We have shown that under certain assumptions, our method is able to estimate a ice thickness consistent with the calving law and ice deformation. Of course, there are many uncertainties involved, and it is not our attempt to find the “true” calving flux here. See Recinos et al., (2019) for our motivations and a discussion about why it is crucial to include frontal ablation in regional ice thickness estimates.

What’s next?#

return to the OGGM documentation
back to the table of contents

OGGM tutorials

Ice thickness inversion with frontal ablation: a case study

Contents