"""Forward model for a (melt_f, prcp_fac, temp_bias) calibration of one glacier, on the v7 `sw_pf3.5_mf5_lam12` gdirs of the lmr_2025 project. One call = one run from the RGI-date flowlines in YS (1500) to the end of the climate file (2024 -> state at 2025), with the climate already in the gdir (HIST ens010: LMR before 1850, shared HIST T + ens148 precip after). It returns the three quantities to match: - dmdtda 2000-2020 (kg m-2 yr-1): Hugonnet, from observations.yml. Computed as OGGM's dynamic calibration does: mass change / RGI area. Also returned per *modelled* area (the model area, which is what a too small glacier hides). - area at the RGI date (km2): the RGI outline area. - volume at the RGI date (km3): the consensus-calibrated inversion volume (observations.yml, ref_volume_m3, year 2004 = start of 2004, i.e. the end of the RGI year 2003 in OGGM's time convention). Usage (any OGGM with incremental gdirs, i.e. master >= 5e9f1332): python aletsch_forward.py # reference checks below python aletsch_forward.py 6.87 2.8 -0.085 # one run: melt_f prcp_fac temp_bias Reference output (OGGM master 5e9f1332, 2026-10-30), state at time 2004: static mf 6.87 pf 2.80 tb -0.085 -> dmdtda -1658, A 84.9 km2, V 13.40 km3 dyn mf 14.15 (round5 calibration) -> dmdtda -1224, A 56.8 km2, V 4.00 km3 targets: dmdtda -1210.7 +- 131.5, A 82.2 km2, V 13.21 km3 """ import os import re import sys import tempfile import numpy as np import oggm.cfg as cfg from oggm import utils, workflow, tasks RGI_ID = os.environ.get('RGI_ID', 'RGI60-11.01450') # Great Aletsch # Gdir tree and experiment: ADA_EXP=sw_pf2.5_mf5_lam12 for the v7 prcp_fac 2.5 # prior; ADA_TREE=lmr_gdirs_v8_202610 ADA_EXP=pf3.5_mf5_lam12 for LAPrec precip TREE = os.environ.get('ADA_TREE', 'lmr_gdirs_v7_202610') EXP = os.environ.get('ADA_EXP', 'sw_pf3.5_mf5_lam12') GDIR_URL = ('https://cluster.klima.uni-bremen.de/~fmaussion/runs/lmr_2025/' f'{TREE}/exps/{EXP}/') PRCP_FAC = float(re.search(r'pf([0-9.]+)_', EXP).group(1)) MELT_F = float(re.search(r'mf([0-9.]+)_', EXP).group(1)) LAMBDA = float(re.search(r'lam([0-9.]+)', EXP).group(1)) BORDER = 240 YS = int(os.environ.get('YS', 1500)) REF_PERIOD = (2000, 2020) def setup(working_dir=None): # Same PARAMS as the gdirs were made with (params.cfg of the experiment). # The calibrated values (melt_f, prcp_fac, temp_bias, glen_a) live in each # gdir's settings.yml and take precedence anyway. cfg.initialize(logging_level='WARNING') cfg.PARAMS['map_proj'] = 'utm' cfg.PARAMS['prcp_fac'] = PRCP_FAC cfg.PARAMS['prcp_fac_min'] = 0.1 cfg.PARAMS['prcp_fac_max'] = 10 cfg.PARAMS['melt_f'] = MELT_F cfg.PARAMS['trapezoid_lambdas'] = LAMBDA cfg.PARAMS['border'] = BORDER cfg.PARAMS['use_multiprocessing'] = False cfg.PARAMS['store_diagnostic_variables'] = ['volume', 'volume_bsl', 'area', 'area_min_h', 'length', 'mass'] cfg.PARAMS['min_ice_thick_for_area'] = 2. cfg.PARAMS['min_ice_thick_for_length'] = 2. cfg.PATHS['working_dir'] = working_dir or tempfile.mkdtemp(prefix='adapbs_') # The L3 is a delta on the OGGM L2 (elev_bands_lambda12): OGGM fetches # the parents named in L3.manifest.json itself gdir = workflow.init_glacier_directories([RGI_ID], from_prepro_level=3, prepro_base_url=GDIR_URL, prepro_rgi_version='62', prepro_border=BORDER)[0] return gdir def targets(gdir): obs = gdir.observations v = obs['ref_volume_m3'] return {'dmdtda': float(obs['ref_mb']['value']), 'dmdtda_err': float(obs['ref_mb']['err']), 'area_km2': gdir.rgi_area_km2, 'vol_km3': float(v['value']) * 1e-9, 'state_year': int(v['year'])} def _v(ds, name): """Run output variable, whatever the OGGM naming (volume / volume_m3...).""" for n in (name, {'volume': 'volume_m3', 'area': 'area_m2', 'mass': 'mass_kg'}[name]): if n in ds: return ds[n] raise KeyError(f'{name} not in the run output: {list(ds.data_vars)}') def forward(gdir, melt_f, prcp_fac, temp_bias, ys=YS, suffix='_fwd'): """One run; returns the matched quantities and the run dataset.""" gdir.settings['melt_f'] = float(melt_f) gdir.settings['prcp_fac'] = float(prcp_fac) gdir.settings['temp_bias'] = float(temp_bias) tasks.run_from_climate_data(gdir, ys=ys, climate_filename='climate_historical', output_filesuffix=suffix) ds = utils.compile_run_output(gdir, input_filesuffix=suffix, path=False) ds = ds.isel(rgi_id=0) ds['time'] = ds.time.astype(int) y0, y1 = REF_PERIOD mass, area, vol = _v(ds, 'mass'), _v(ds, 'area'), _v(ds, 'volume') dm = float(mass.sel(time=y1) - mass.sel(time=y0)) a_mod = float(area.sel(time=slice(y0, y1 - 1)).mean()) ty = targets(gdir)['state_year'] out = {'melt_f': melt_f, 'prcp_fac': prcp_fac, 'temp_bias': temp_bias, 'dmdtda_rgi_area': dm / gdir.rgi_area_m2 / (y1 - y0), # = OGGM dyn calib 'dmdtda_model_area': dm / a_mod / (y1 - y0), 'area_km2': float(area.sel(time=ty)) * 1e-6, 'vol_km3': float(vol.sel(time=ty)) * 1e-9, 'area_1850_km2': float(area.sel(time=1850)) * 1e-6, 'vol_1850_km3': float(vol.sel(time=1850)) * 1e-9} return out, ds def fmt(o): return (f"mf {o['melt_f']:6.2f} pf {o['prcp_fac']:4.2f} tb {o['temp_bias']:+5.2f} | " f"dmdtda {o['dmdtda_rgi_area']:7.0f} (per model area {o['dmdtda_model_area']:7.0f}) | " f"A {o['area_km2']:6.1f} km2 V {o['vol_km3']:6.2f} km3 | " f"1850: A {o['area_1850_km2']:6.1f} V {o['vol_1850_km3']:6.2f}") if __name__ == '__main__': gdir = setup(os.environ.get('WORKDIR')) t = targets(gdir) s = gdir.settings print(f"{RGI_ID}: targets dmdtda {t['dmdtda']} +- {t['dmdtda_err']}, area {t['area_km2']:.1f} km2, " f"volume {t['vol_km3']:.2f} km3 (state year {t['state_year']})") print(f"static calibration in the gdir: melt_f {s['melt_f']:.3f}, prcp_fac {s['prcp_fac']:.3f}, " f"temp_bias {s['temp_bias']:+.3f}") if len(sys.argv) == 4: print(fmt(forward(gdir, *map(float, sys.argv[1:]))[0])) else: mf0, pf0, tb0 = s['melt_f'], s['prcp_fac'], s['temp_bias'] print(fmt(forward(gdir, mf0, pf0, tb0)[0]), ' <- static calibration') print(fmt(forward(gdir, 14.15146, pf0, tb0)[0]), ' <- round5 dyn melt_f (reference)')