Dynamic model initialization using observed thickness data#
This notebook demonstrates the translation of thickness observations into bed topography in ‘flowline-space’ and a dynamic model initialization for projections. Guidelines for adding thickness observations to your glacier directory are provided in the tutorials Ingest gridded products such as ice velocity into OGGM and OGGM-Shop and Glacier Directories in OGGM.
from oggm import utils, workflow, tasks, cfg, graphics
import pandas as pd
import matplotlib.pyplot as plt
import copy
import numpy as np
import xarray as xr
from scipy import optimize
Translate thickness observations into elevation band flowlines#
To convert thickness observations into bed topography, the initial step involves binning the data into elevation bands, as detailed in the tutorial Ingest gridded products such as ice velocity into OGGM. Fortunately, preprocessed directories are available, encompassing all data supported by the OGGM-Shop, already binned to elevation bands. This enables easy initialization of a dynamic flowline using tasks.init_present_time_glacier. Specify the data to be used with use_binned_thickness_data. Here, dynamic flowlines are defined using consensus thickness (Farinotti et al. 2019) and Millan thickness data (Millan et al. 2022).
# initialize our test glacier from the preprocessed directories including the shop data
cfg.initialize(logging_level='WARNING')
cfg.PATHS['working_dir'] = utils.gettempdir('OGGM_thk_dyn_spn', reset=True)
rgi_ids = ['RGI60-11.00897']
# preprocessed directories including shop data and the default oggm dynamic initialisation
base_url = 'https://cluster.klima.uni-bremen.de/~oggm/gdirs/oggm_v1.6/L3-L5_files/2023.3/elev_bands/W5E5_spinup_w_data/'
gdirs = workflow.init_glacier_directories(
rgi_ids, # which glaciers?
prepro_base_url=base_url, # where to fetch the data?
from_prepro_level=4, # what kind of data?
prepro_border=80 # how big of a map?
)
gdir = gdirs[0]
2024-08-25 21:31:26: oggm.cfg: Reading default parameters from the OGGM `params.cfg` configuration file.
2024-08-25 21:31:26: oggm.cfg: Multiprocessing switched OFF according to the parameter file.
2024-08-25 21:31:26: oggm.cfg: Multiprocessing: using all available processors (N=4)
2024-08-25 21:31:28: oggm.workflow: init_glacier_directories from prepro level 4 on 1 glaciers.
2024-08-25 21:31:28: oggm.workflow: Execute entity tasks [gdir_from_prepro] on 1 glaciers
# This cell is only needed due to a bug in oggm v1.6.1 which was fixed and can be deleted with updated preprocessed directories
bin_variables = ['consensus_ice_thickness', 'millan_ice_thickness']
workflow.execute_entity_task(tasks.elevation_band_flowline,
gdirs,
bin_variables=bin_variables)
workflow.execute_entity_task(tasks.fixed_dx_elevation_band_flowline,
gdirs,
bin_variables=bin_variables);
2024-08-25 21:31:28: oggm.workflow: Execute entity tasks [elevation_band_flowline] on 1 glaciers
2024-08-25 21:31:28: oggm.workflow: Execute entity tasks [fixed_dx_elevation_band_flowline] on 1 glaciers
# create the dynamic flowlines for consensus and millan thickness
tasks.init_present_time_glacier(gdir, filesuffix='_consensus',
use_binned_thickness_data='consensus_ice_thickness')
tasks.init_present_time_glacier(gdir, filesuffix='_millan',
use_binned_thickness_data='millan_ice_thickness')
# open the default oggm flowline and the two we just created
fls_oggm = gdir.read_pickle('model_flowlines')
fls_consensus = gdir.read_pickle('model_flowlines', filesuffix='_consensus')
fls_millan = gdir.read_pickle('model_flowlines', filesuffix='_millan')
# plot the created dynamic flowlines, including the oggm default
fig, axs = plt.subplots(2, 2, figsize=(10, 7))
fig.subplots_adjust(wspace=0.6, hspace=0.5)
graphics.plot_modeloutput_section(fls_oggm, ax=axs[0][0], title='OGGM default')
graphics.plot_modeloutput_section(fls_consensus, ax=axs[0][1], title='Consensus')
graphics.plot_modeloutput_section(fls_millan, ax=axs[1][0], title='Millan')
ax_compare = axs[1][1]
ax_compare.plot(fls_oggm[0].surface_h, label='surface height') # surface height is the same for all flowlines
ax_compare.plot(fls_oggm[0].bed_h, label='OGGM default')
ax_compare.plot(fls_consensus[0].bed_h, label='Consensus')
ax_compare.plot(fls_millan[0].bed_h, label='Millan')
ax_compare.set_title('Direct comparison')
ax_compare.legend();
The resulting flowlines exhibit identical surface and area height distributions, defined by the outline and the DEM. Comparing individual bed heights, there is minimal difference between OGGM default and consensus flowline bed heights. This is attributed to OGGM’s past contribution to the consensus estimate and the use of the consensus volume estimate in OGGM inversion. In contrast, the Millan bed height differs significantly.
It’s crucial to note that these newly created flowlines maintain the total glacier volume of the observations. While OGGM default’s glacier volume is regionally calibrated, differences are expected compared to the consensus volume estimate, even though the consensus was used in the inversion (see next section).
The key question now is how to dynamically calibrate and initialize the model with these newly generated flowlines.
Dynamically calibrate and initialise todays glacier state, using flowlines from thickness observations#
Now, we delve into the process of initializing the model using OGGM’s dynamic calibration methodology, elaborated in the tutorial Dynamic spinup and dynamic melt_f calibration for past simulations. In summary, the approach involves dynamically aligning the RGI area by identifying a suitable glacier state in the past and dynamically adjusting the melt factor of the mass balance to match observed geodetic mass balance. With the glacier bed defined by observations, there is an additional capability to dynamically match for the volume (further details below).
To commence, use the dynamic calibration function run_dynamic_melt_f_calibration in conjunction with the newly created flowlines, which can be passed to init_model_fls. Refer to the designated tutorial Dynamic spinup and dynamic melt_f calibration for past simulations for more details on the other parameters.
err_dmdtda_scaling_factor = 0.2 # the default value of oggm runs to account for correlated observation errors
# Consensus flowline
workflow.execute_entity_task(
tasks.run_dynamic_melt_f_calibration, gdirs,
init_model_fls=fls_consensus,
err_dmdtda_scaling_factor=err_dmdtda_scaling_factor,
ys=1979, ye=2020,
melt_f_max=cfg.PARAMS['melt_f_max'],
kwargs_run_function={'minimise_for': 'area',
'do_inversion': False},
ignore_errors=True,
kwargs_fallback_function={'minimise_for': 'area',
'do_inversion': False},
output_filesuffix='_spinup_historical_consensus',
)
# Millan flowlines
workflow.execute_entity_task(
tasks.run_dynamic_melt_f_calibration, gdirs,
init_model_fls=fls_millan,
err_dmdtda_scaling_factor=err_dmdtda_scaling_factor,
ys=1979, ye=2020,
melt_f_max=cfg.PARAMS['melt_f_max'],
kwargs_run_function={'minimise_for': 'area',
'do_inversion': False},
ignore_errors=True,
kwargs_fallback_function={'minimise_for': 'area',
'do_inversion': False},
output_filesuffix='_spinup_historical_millan',
)
2024-08-25 21:31:29: oggm.workflow: Execute entity tasks [run_dynamic_melt_f_calibration] on 1 glaciers
2024-08-25 21:31:31: oggm.workflow: Execute entity tasks [run_dynamic_melt_f_calibration] on 1 glaciers
[<oggm.core.flowline.SemiImplicitModel at 0x7f652728d550>]
# define a plotting function to compare the dynamic model runs with the observations
def compare_model_to_observations(gdir, filesuffixes):
area_ref = gdir.rgi_area_m2
area_ref_error = area_ref * 0.01 # 1% is the default of the dynamic spinup
df_dmdtda_ref = utils.get_geodetic_mb_dataframe().loc[gdir.rgi_id]
sel = df_dmdtda_ref.loc[df_dmdtda_ref['period'] == '2000-01-01_2020-01-01'].iloc[0]
# reference geodetic mass balance from Hugonnet 2021
dmdtda_ref = float(sel['dmdtda'])
# dmdtda: in meters water-equivalent per year -> we convert
dmdtda_ref *= 1000 # kg m-2 yr-1
# error of reference geodetic mass balance from Hugonnet 2021
dmdtda_ref_error = float(sel['err_dmdtda'])
dmdtda_ref_error *= 1000 # kg m-2 yr-1
dmdtda_ref_error *= err_dmdtda_scaling_factor
dmdtda_comparison = f'dmdtda reference (with reduced uncertainty): {dmdtda_ref:.1f} +/- {dmdtda_ref_error:.1f} kg m-2 yr-1'
def get_dmdtda_modelled(volume, area_ref, yr0_ref_mb=2000, yr1_ref_mb=2020):
return ((volume.loc[yr1_ref_mb].values -
volume.loc[yr0_ref_mb].values) /
area_ref /
(yr1_ref_mb - yr0_ref_mb) *
cfg.PARAMS['ice_density'])
# get the data of all model runs and plot it
fig, axs = plt.subplots(2, 1, figsize=(10, 6))
for i, filesuffix in enumerate(filesuffixes):
if filesuffix == '':
label = 'OGGM default'
else:
label = filesuffix[1:]
fl = gdir.read_pickle('model_flowlines', filesuffix=filesuffix)[0]
with xr.open_dataset(
gdir.get_filepath('model_diagnostics',
filesuffix=f'_spinup_historical{filesuffix}')) as ds:
ds_tmp = ds.load()
# plot volume with reference
ref_year = gdir.rgi_date + 1
axs[0].plot(ds_tmp.time.values, ds_tmp.volume_m3.values,
c=f'C{i}', label=label)
axs[0].plot(ref_year, fl.volume_m3,
'o', c=f'C{i}', label=f'ref volume {label}')
# plot area with reference
axs[1].plot(ds_tmp.time, ds_tmp.area_m2,
c=f'C{i}', label=label)
axs[1].plot(ref_year, area_ref,
'o', c=f'C{i}', label='ref area')
axs[1].plot([ref_year, ref_year],
[area_ref - area_ref_error,
area_ref + area_ref_error],
c=f'C{i}')
# add dmdtda as text
dmdtda_tmp = get_dmdtda_modelled(ds_tmp.volume_m3, area_ref)
dmdtda_comparison += (f'\n{f"dmdtda {label}:":<44} {dmdtda_tmp:.1f} kg m-2 yr-1, diff. to ref.: '
f'{f"{dmdtda_tmp - dmdtda_ref:.1f}":>5} kg m-2 yr-1')
axs[0].legend()
axs[0].set_ylabel('Volume [m³]')
axs[1].legend()
axs[1].set_ylabel('Area [m²]')
print(dmdtda_comparison)
compare_model_to_observations(gdir, ['', '_consensus', '_millan'])
dmdtda reference (with reduced uncertainty): -1100.3 +/- 34.4 kg m-2 yr-1
dmdtda OGGM default: -1130.1 kg m-2 yr-1, diff. to ref.: -29.8 kg m-2 yr-1
dmdtda consensus: -1078.0 kg m-2 yr-1, diff. to ref.: 22.3 kg m-2 yr-1
dmdtda millan: -1134.6 kg m-2 yr-1, diff. to ref.: -34.3 kg m-2 yr-1
In all three simulations, we observe that both the area and geodetic mass balance align within the target boundaries. Nevertheless, the volume is only coincidentally matched for the consensus estimate. This discrepancy arises because the current deformation parameter was calibrated during the inversion for OGGM default initialization, which incorporates an equilibrium assumption (see documention for more information). However, when defining the glacier bed from thickness observations, it becomes possible/necessary to calibrate the deformation parameter to match the observed volume during initialization. Although there is currently no implemented function for this, the following code should provide an idea of how it can be achieved:
# create a other flowline for this, for later comparision
tasks.init_present_time_glacier(gdir, filesuffix='_millan_adapted',
use_binned_thickness_data='millan_ice_thickness')
def to_minimize(glen_a):
ref_volume = fls_millan[0].volume_m3
workflow.execute_entity_task(
tasks.run_dynamic_melt_f_calibration, gdirs,
init_model_fls=gdir.read_pickle('model_flowlines',
filesuffix='_millan_adapted'),
err_dmdtda_scaling_factor=0.2,
ys=1979, ye=2020,
melt_f_max=cfg.PARAMS['melt_f_max'],
kwargs_run_function={'minimise_for': 'area',
'do_inversion': False,
'glen_a': glen_a},
ignore_errors=True,
kwargs_fallback_function={'minimise_for': 'area',
'do_inversion': False},
output_filesuffix='_spinup_historical_millan_adapted',
)
with xr.open_dataset(gdir.get_filepath('model_diagnostics',
filesuffix='_spinup_historical_millan_adapted')) as ds:
ds_millan_adapted = ds.load()
return ref_volume - ds_millan_adapted.loc[{'time': gdir.rgi_date + 1}].volume_m3.values.item()
glen_a_millan = optimize.brentq(to_minimize, 1e-22, 1e-24, xtol=1e-25)
# conduct the dynamic initialisation once more to be sure the right one is saved
workflow.execute_entity_task(
tasks.run_dynamic_melt_f_calibration, gdirs,
init_model_fls=gdir.read_pickle('model_flowlines', filesuffix='_millan_adapted'),
err_dmdtda_scaling_factor=0.2,
ys=1979, ye=2020,
melt_f_max=cfg.PARAMS['melt_f_max'],
kwargs_run_function={'minimise_for': 'area',
'do_inversion': False,
'glen_a': glen_a_millan},
ignore_errors=True,
kwargs_fallback_function={'minimise_for': 'area',
'do_inversion': False,
'glen_a': glen_a_millan},
output_filesuffix='_spinup_historical_millan_adapted',
)
compare_model_to_observations(gdir, ['', '_consensus', '_millan', '_millan_adapted'])
dmdtda reference (with reduced uncertainty): -1100.3 +/- 34.4 kg m-2 yr-1
dmdtda OGGM default: -1130.1 kg m-2 yr-1, diff. to ref.: -29.8 kg m-2 yr-1
dmdtda consensus: -1078.0 kg m-2 yr-1, diff. to ref.: 22.3 kg m-2 yr-1
dmdtda millan: -1134.6 kg m-2 yr-1, diff. to ref.: -34.3 kg m-2 yr-1
dmdtda millan_adapted: -1110.3 kg m-2 yr-1, diff. to ref.: -10.0 kg m-2 yr-1
2024-08-25 21:31:35: oggm.workflow: Execute entity tasks [run_dynamic_melt_f_calibration] on 1 glaciers
2024-08-25 21:31:39: oggm.workflow: Execute entity tasks [run_dynamic_melt_f_calibration] on 1 glaciers
2024-08-25 21:31:41: oggm.workflow: Execute entity tasks [run_dynamic_melt_f_calibration] on 1 glaciers
2024-08-25 21:31:47: oggm.workflow: Execute entity tasks [run_dynamic_melt_f_calibration] on 1 glaciers
2024-08-25 21:31:49: oggm.workflow: Execute entity tasks [run_dynamic_melt_f_calibration] on 1 glaciers
2024-08-25 21:31:50: oggm.workflow: Execute entity tasks [run_dynamic_melt_f_calibration] on 1 glaciers
2024-08-25 21:31:54: oggm.workflow: Execute entity tasks [run_dynamic_melt_f_calibration] on 1 glaciers
2024-08-25 21:31:56: oggm.workflow: Execute entity tasks [run_dynamic_melt_f_calibration] on 1 glaciers
2024-08-25 21:31:58: oggm.workflow: Execute entity tasks [run_dynamic_melt_f_calibration] on 1 glaciers
2024-08-25 21:31:59: oggm.workflow: Execute entity tasks [run_dynamic_melt_f_calibration] on 1 glaciers
2024-08-25 21:32:01: oggm.workflow: Execute entity tasks [run_dynamic_melt_f_calibration] on 1 glaciers
2024-08-25 21:32:03: oggm.workflow: Execute entity tasks [run_dynamic_melt_f_calibration] on 1 glaciers
2024-08-25 21:32:05: oggm.workflow: Execute entity tasks [run_dynamic_melt_f_calibration] on 1 glaciers
Now, we achieve dynamic consistency by calibrating the Glen A factor, ensuring a match in area, geodetic mass balance and volume.
glen_a_millan
1.7797147426841996e-23
Note that if the optimized Glen A factor falls outside expected boundaries, considering the use of a sliding factor fs is advisable.
Important: To retain the optimized Glen A (and fs) during projections, it must be set in the glacier directory.
# change glen_a in gdir for potential projection runs
gdir.add_to_diagnostics('inversion_glen_a', glen_a_millan)
# gdir.add_to_diagnostics('inversion_fs', fs_optimized)
With this configuration, you can conduct projections as outlined in the tutorial 10 minutes to… a glacier change projection with GCM data.
If you intend to undertake a comparative study using different thickness datasets, we highly recommend creating an individual glacier directory per thickness data to ensure consistent use of the correct dynamic parameters.
What’s next?#
Look at the more comprehensive tutorial Dynamic spinup and dynamic melt_f calibration for past simulations
return to the OGGM documentation
back to the table of contents