diff --git a/pygem/bin/postproc/postproc_binned_subannual_thick.py b/pygem/bin/postproc/postproc_binned_subannual_thick.py index 13479ae3..b7c37821 100644 --- a/pygem/bin/postproc/postproc_binned_subannual_thick.py +++ b/pygem/bin/postproc/postproc_binned_subannual_thick.py @@ -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 @@ -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). @@ -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) diff --git a/pygem/bin/run/run_calibration.py b/pygem/bin/run/run_calibration.py index 1899e272..20bae6f3 100755 --- a/pygem/bin/run/run_calibration.py +++ b/pygem/bin/run/run_calibration.py @@ -317,6 +317,7 @@ 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 @@ -324,19 +325,24 @@ def run_oggm_dynamics(gdir, modelprms, glacier_rgi_table, fls): 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 @@ -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] @@ -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') @@ -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 @@ -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'])} @@ -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']: @@ -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.] @@ -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']] diff --git a/pygem/bin/run/run_simulation.py b/pygem/bin/run/run_simulation.py index 3a3bc0fd..fafcaf26 100755 --- a/pygem/bin/run/run_simulation.py +++ b/pygem/bin/run/run_simulation.py @@ -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, @@ -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, ) @@ -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] @@ -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 @@ -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 @@ -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,