Skip to content

Store flowline diagnostics from OGGM dynamical model (#167) #168

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
btobers wants to merge 1 commit into from
Closed

Store flowline diagnostics from OGGM dynamical model (#167) #168

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
26 changes: 14 additions & 12 deletions pygem/bin/postproc/postproc_binned_subannual_thick.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -67,7 +67,12 @@ def getparser():


def get_binned_subannual(
bin_massbalclim, bin_mass_annual, bin_thick_annual, dates_subannual, dates_annual, debug=False
bin_massbalclim,
bin_mass_annual,
bin_thick_annual,
dates_subannual,
dates_annual,
debug=False,
):
"""
funciton to calculate the subannual binned ice thickness and mass
Expand All @@ -79,7 +84,7 @@ def get_binned_subannual(
thus, subannual thickness and mass is determined assuming
the flux divergence is constant throughout the year.

annual flux divergence is first estimated by combining the annual binned change in ice
annual flux divergence is first calculated by combining the annual binned change in ice
thickness and the annual binned mass balance. then, assume flux divergence is constant
throughout the year (divide annual by the number of steps in the binned climatic mass
balance to get subannual flux divergence).
Expand Down Expand Up @@ -118,30 +123,27 @@ def get_binned_subannual(
shape : [#glac, #elevbins, #steps]
"""

n_glac, n_bins, n_steps = bin_massbalclim.shape
years_annual = np.array([d.year for d in dates_annual])
years_subannual = np.array([d.year for d in dates_subannual])
yrs = np.unique(years_subannual)
nyrs = len(yrs)
assert nyrs > 1, 'Need at least two annual steps for flux divergence estimation'
assert len(yrs) > 1, 'Need at least two annual steps for flux divergence estimation'

# --- Step 1: convert mass balance from m w.e. to m ice ---
rho_w = pygem_prms['constants']['density_water']
rho_i = pygem_prms['constants']['density_ice']
bin_massbalclim_ice = bin_massbalclim * (rho_w / rho_i)

# --- Step 2: compute annual cumulative mass balance ---
# --- Step 2: compute annual thickness change ---
bin_delta_thick_annual = np.diff(bin_thick_annual, axis=-1) # [m ice]

# --- Step 3: compute annual cumulative mass balance ---
# Initialize annual cumulative mass balance (exclude last year for flux calculation)
bin_massbalclim_annual = np.zeros((n_glac, n_bins, nyrs))
bin_massbalclim_annual = np.zeros_like(bin_delta_thick_annual)
for i, year in enumerate(yrs):
idx = np.where(years_subannual == year)[0]
bin_massbalclim_annual[:, :, i] = bin_massbalclim_ice[:, :, idx].sum(axis=-1)

# --- Step 3: compute annual thickness change ---
bin_delta_thick_annual = np.diff(bin_thick_annual, axis=-1) # [m ice yr^-1]

# --- Step 4: compute annual flux divergence ---
bin_flux_divergence_annual = bin_massbalclim_annual - bin_delta_thick_annual # [m ice yr^-1]
bin_flux_divergence_annual = bin_massbalclim_annual - bin_delta_thick_annual # [m ice]

# --- Step 5: expand flux divergence to subannual steps ---
bin_flux_divergence_subannual = np.zeros_like(bin_massbalclim_ice)
Expand Down
49 changes: 34 additions & 15 deletions pygem/bin/run/run_calibration.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -317,26 +317,32 @@ def run_oggm_dynamics(gdir, modelprms, glacier_rgi_table, fls):
ev_model.mb_model.glac_wide_massbaltotal = (
ev_model.mb_model.glac_wide_massbaltotal + ev_model.mb_model.glac_wide_frontalablation
)
ev_model.mb_model.ensure_mass_conservation(diag)
# safely catch any errors with dynamical run
except Exception:
ds = None

return mbmod, ds


def calc_thick_change_1d(gdir, mbmod, ds):
def calc_thick_change_1d(gdir, fls, mbmod, ds):
"""
calculate binned change in ice thickness assuming constant annual flux divergence.
sub-annual ice thickness is differenced at timesteps coincident with observations.
"""
years_subannual = np.array([d.year for d in gdir.dates_table['date']])
yrs = np.unique(years_subannual)
nyrs = len(yrs)
# grab components of interest
bin_thick_annual = ds[0].thickness_m.values.T # glacier thickness [m ice], (nbins, nyears)

# set any < 0 thickness to nan
bin_thick_annual[bin_thick_annual <= 0] = np.nan
bin_massbalclim_ice_annual = ds[0].climatic_mb_myr.values.T[
:, 1:
] # annual climatic mass balance [m ice] (nbins, nyears - 1)
bin_delta_thick_annual = ds[0].dhdt_myr.values.T[
:, 1:
] # annual change in ice thickness [m ice] (nbins, nyears - 1)
bin_flux_divergence_annual = (
bin_massbalclim_ice_annual - bin_delta_thick_annual
) # annual flux divergence [m ice] (nbins, nyears - 1)

# --- Step 1: convert mass balance from m w.e. to m ice ---
bin_massbalclim = mbmod.glac_bin_massbalclim # climatic mass balance [m w.e.] per step
Expand All @@ -348,7 +354,6 @@ def calc_thick_change_1d(gdir, mbmod, ds):
# --- Step 2: expand flux divergence to subannual steps ---
# assume flux divergence is constant throughout the year
# (divide annual by the number of steps in the binned climatic mass balance to get subannual flux divergence)
bin_flux_divergence_annual = -ds[0].flux_divergence_myr.values.T[:, 1:]
bin_flux_divergence_subannual = np.zeros_like(bin_massbalclim_ice)
for i, year in enumerate(yrs):
idx = np.where(years_subannual == year)[0]
Expand All @@ -358,12 +363,16 @@ def calc_thick_change_1d(gdir, mbmod, ds):
bin_delta_thick_subannual = bin_massbalclim_ice - bin_flux_divergence_subannual

# --- Step 4: calculate subannual thickness = running thickness change + initial thickness---
running_bin_delta_thick_subannual = np.cumsum(bin_delta_thick_subannual, axis=-1)
running_bin_delta_thick_subannual = np.nancumsum(bin_delta_thick_subannual, axis=-1)
bin_thick_subannual = running_bin_delta_thick_subannual + bin_thick_annual[:, 0][:, np.newaxis]

# --- Step 5: rebin ---
# get surface height at the specified reference year
ref_surface_height = ds[0].bed_h.values + ds[0].thickness_m.sel(time=gdir.elev_change_1d['ref_dem_year']).values
ref_surface_height = (
fls[0].surface_h
- ds[0].thickness_m.isel(time=0).values
+ ds[0].thickness_m.sel(time=gdir.elev_change_1d['ref_dem_year']).values
)
# aggregate model bin thicknesses as desired
with warnings.catch_warnings():
warnings.filterwarnings('ignore')
Expand Down Expand Up @@ -416,7 +425,7 @@ def mcmc_model_eval(
if calib_elev_change_1d:
mbmod, ds = run_oggm_dynamics(gdir, modelprms, glacier_rgi_table, fls)
# note, the binned thickness change is scaled by modeled density in mcmc.mbPosterior.log_likelihood() to calculate modeled surface elevation change
results['elev_change_1d'] = calc_thick_change_1d(gdir, mbmod, ds) if ds else float('-inf')
results['elev_change_1d'] = calc_thick_change_1d(gdir, fls, mbmod, ds) if ds else float('-inf')

if mbfxn is not None:
# grab current values from modelprms for the emulator
Expand Down Expand Up @@ -859,6 +868,19 @@ def run(list_packed_vars):
if not isinstance(gdir.elev_change_1d['dh_sigma'], int)
else gdir.elev_change_1d['dh_sigma']
)

# optionally adjust ref_startyear based on earliest available elevation calibration data (must be <= 2000)
if args.spinup:
args.ref_startyear = min(
2000, *(int(date[:4]) for pair in gdir.elev_change_1d['dates'] for date in pair)
)
# adjust dates table and climate data
mask = gdir.dates_table['year'] >= args.ref_startyear
gdir.dates_table = gdir.dates_table.loc[mask].reset_index(drop=True)
nrows_rm = (~mask).sum() # number of dropped rows
for key in ('temp', 'tempstd', 'prec', 'lr'):
gdir.historical_climate[key] = gdir.historical_climate[key][nrows_rm:]

# get observation period indices in model date_table
# create lookup dict (timestamp → index)
date_to_index = {d: i for i, d in enumerate(gdir.dates_table['date'])}
Expand All @@ -869,11 +891,6 @@ def run(list_packed_vars):
)
for start, end in gdir.elev_change_1d['dates']
]
# optionally adjust ref_startyear based on earliest available elevation calibration data (must be <= 2000)
if args.spinup:
args.ref_startyear = min(
2000, *(int(date[:4]) for pair in gdir.elev_change_1d['dates'] for date in pair)
)

# load calibrated calving_k values for tidewater glaciers
if gdir.is_tidewater and pygem_prms['setup']['include_frontalablation']:
Expand Down Expand Up @@ -1884,7 +1901,7 @@ def calc_mb_total_minelev(modelprms):
]
# Melt
# energy available for melt [degC day]
melt_energy_available = T_minelev * dates_table['days_in_step'].values
melt_energy_available = T_minelev * gdir.dates_table['days_in_step'].values
melt_energy_available[melt_energy_available < 0] = 0
# assume all snow melt because anything more would melt underlying ice in lowermost bin
# SNOW MELT [m w.e.]
Expand Down Expand Up @@ -2135,6 +2152,8 @@ def rho_constraints(**kwargs):
)
# prepare export modelprms dictionary
modelprms_export = {}
# store reference period
modelprms_export['ref_period'] = (args.ref_startyear, args.ref_endyear)
# store model parameters and priors
modelprms_export['precgrad'] = [pygem_prms['sim']['params']['precgrad']]
modelprms_export['tsnow_threshold'] = [pygem_prms['sim']['params']['tsnow_threshold']]
Expand Down
50 changes: 15 additions & 35 deletions pygem/bin/run/run_simulation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -976,7 +976,7 @@ def run(list_packed_vars):
nfls,
y0=args.sim_startyear,
mb_model=mbmod,
glen_a=cfg.PARAMS['glen_a'] * glen_a_multiplier,
glen_a=glen_a,
fs=fs,
is_tidewater=gdir.is_tidewater,
water_level=water_level,
Expand All @@ -987,7 +987,7 @@ def run(list_packed_vars):
nfls,
y0=args.sim_startyear,
mb_model=mbmod,
glen_a=cfg.PARAMS['glen_a'] * glen_a_multiplier,
glen_a=glen_a,
fs=fs,
)

Expand All @@ -997,7 +997,7 @@ def run(list_packed_vars):
ev_model, ax=ax, lnlabel=f'Glacier year {args.sim_startyear}'
)

diag = ev_model.run_until_and_store(args.sim_endyear + 1)
diag, ds = ev_model.run_until_and_store(args.sim_endyear + 1, fl_diag_path=True)
ev_model.mb_model.glac_wide_volume_annual[-1] = diag.volume_m3[-1]
ev_model.mb_model.glac_wide_area_annual[-1] = diag.area_m2[-1]

Expand Down Expand Up @@ -1231,8 +1231,16 @@ def run(list_packed_vars):
output_offglac_melt_steps[:, n_iter] = mbmod.offglac_wide_melt
output_offglac_snowpack_steps[:, n_iter] = mbmod.offglac_wide_snowpack
output_offglac_runoff_steps[:, n_iter] = mbmod.offglac_wide_runoff

if output_glac_bin_icethickness_annual is None:
# binned ouputs
if args.option_dynamics == 'OGGM':
# grab binned outputs from oggm flowline diagnostics
# note, transpose restructures (time, dis_along_flowline) -> (dis_along_flowline, time)
output_glac_bin_area_annual_sim = ds[0].area_m2.values.T[:, :, np.newaxis]
output_glac_bin_icethickness_annual_sim = ds[0].thickness_m.values.T[:, :, np.newaxis]
output_glac_bin_mass_annual_sim = (
ds[0].volume_m3.values.T[:, :, np.newaxis] * pygem_prms['constants']['density_ice']
)
else:
output_glac_bin_area_annual_sim = mbmod.glac_bin_area_annual[:, :, np.newaxis]
output_glac_bin_mass_annual_sim = (
mbmod.glac_bin_area_annual
Expand Down Expand Up @@ -1263,6 +1271,8 @@ def run(list_packed_vars):
output_glac_bin_mass_annual_sim[:, -1, 0] = (
glacier_vol_t0 * pygem_prms['constants']['density_ice']
)
# append individual simulation outputs together
if output_glac_bin_icethickness_annual is None:
output_glac_bin_area_annual = output_glac_bin_area_annual_sim
output_glac_bin_mass_annual = output_glac_bin_mass_annual_sim
output_glac_bin_icethickness_annual = output_glac_bin_icethickness_annual_sim
Expand All @@ -1288,36 +1298,6 @@ def run(list_packed_vars):
output_glac_bin_melt_steps = output_glac_bin_melt_steps_sim[:, :, np.newaxis]

else:
# Update the latest thickness and volume
output_glac_bin_area_annual_sim = mbmod.glac_bin_area_annual[:, :, np.newaxis]
output_glac_bin_mass_annual_sim = (
mbmod.glac_bin_area_annual
* mbmod.glac_bin_icethickness_annual
* pygem_prms['constants']['density_ice']
)[:, :, np.newaxis]
output_glac_bin_icethickness_annual_sim = (mbmod.glac_bin_icethickness_annual)[
:, :, np.newaxis
]
if ev_model is not None:
fl_dx_meter = getattr(ev_model.fls[0], 'dx_meter', None)
fl_widths_m = getattr(ev_model.fls[0], 'widths_m', None)
fl_section = getattr(ev_model.fls[0], 'section', None)
else:
fl_dx_meter = getattr(nfls[0], 'dx_meter', None)
fl_widths_m = getattr(nfls[0], 'widths_m', None)
fl_section = getattr(nfls[0], 'section', None)
if fl_section is not None and fl_widths_m is not None:
# thickness
icethickness_t0 = np.zeros(fl_section.shape)
icethickness_t0[fl_widths_m > 0] = (
fl_section[fl_widths_m > 0] / fl_widths_m[fl_widths_m > 0]
)
output_glac_bin_icethickness_annual_sim[:, -1, 0] = icethickness_t0
# mass
glacier_vol_t0 = fl_widths_m * fl_dx_meter * icethickness_t0
output_glac_bin_mass_annual_sim[:, -1, 0] = (
glacier_vol_t0 * pygem_prms['constants']['density_ice']
)
output_glac_bin_area_annual = np.append(
output_glac_bin_area_annual,
output_glac_bin_area_annual_sim,
Expand Down