Add option for longer dynamic spinups

docs/whats-new.rst

@@ -12,6 +12,11 @@ Enhancements
 - `run_prepro_levels` now allows to set a custom climate data processing
   module and custom geodetic data file to create glacier directories (:pull:`1910`).
   By `Fabien Maussion <https://github.com/fmaussion>`_
+- New kwarg `spinup_periods_to_try` in `run_dynamic_spinup` to be able to
+  provide a list of additional spinup periods, which are tried in order if both
+  the initially defined spinup period (`spinup_period_initial`) and the minimum
+  spinup period (`min_spinup_period`) fail (:pull:`1914`).
+  By `Patrick Schmitt <https://github.com/pat-schmitt>`_
 
 Bug fixes
 ~~~~~~~~~

oggm/cli/prepro_levels.py

@@ -116,6 +116,7 @@ def run_prepro_levels(rgi_version=None, rgi_reg=None, border=None,
                       logging_level='WORKFLOW',
                       dynamic_spinup=False, err_dmdtda_scaling_factor=0.2,
                       dynamic_spinup_start_year=1979,
+                      dynamic_spinup_periods_to_try=None,
                       continue_on_error=True, store_fl_diagnostics=False):
     """Generate the preprocessed OGGM glacier directories for this OGGM version
 
@@ -239,6 +240,11 @@ def run_prepro_levels(rgi_version=None, rgi_reg=None, border=None,
     dynamic_spinup_start_year : int
         if dynamic_spinup is set, define the starting year for the simulation.
         The default is 1979, unless the climate data starts later.
+    dynamic_spinup_periods_to_try : list or None
+        If the spinup_period defined by rgi_date - dynamic_spinup_start_yr was
+        not successful, you can provide here a list of spinup periods which
+        should be tried in order.
+        Default is None
     continue_on_error : bool
         if True the workflow continues if a task raises an error. For operational
         runs it should be set to True (the default).
@@ -361,7 +367,7 @@ def _time_log():
 
             # In AA almost all large ice bodies are actually ice caps
             if rgi_reg == '19':
-                rgidf.loc[rgidf.Area > 100, 'Form'] = '1'
+                rgidf.loc[rgidf.Area > 100, 'Form'] = 1
 
             # For greenland we omit connectivity level 2
             if rgi_reg == '05':
@@ -902,9 +908,15 @@ def _time_log():
                 err_dmdtda_scaling_factor=err_dmdtda_scaling_factor,
                 ys=dynamic_spinup_start_year, ye=ye,
                 melt_f_max=melt_f_max,
-                kwargs_run_function={'minimise_for': minimise_for},
+                kwargs_run_function={'minimise_for': minimise_for,
+                                     'spinup_periods_to_try':
+                                         dynamic_spinup_periods_to_try
+                                     },
                 ignore_errors=True,
-                kwargs_fallback_function={'minimise_for': minimise_for},
+                kwargs_fallback_function={'minimise_for': minimise_for,
+                                          'spinup_periods_to_try':
+                                              dynamic_spinup_periods_to_try
+                                          },
                 output_filesuffix='_spinup_historical',)
             # Now compile the output
             opath = os.path.join(sum_dir, f'spinup_historical_run_output_{rgi_reg}.nc')
@@ -1125,7 +1137,7 @@ def parse_args(args):
                         help="include a dynamic spinup for matching glacier area "
                              "('area/dmdtda') OR volume ('volume/dmdtda') at "
                              "the RGI-date, AND mass-change from Hugonnet "
-                             "in the period 2000-2020 (dynamic mu* "
+                             "in the period 2000-2020 (dynamic melt_f "
                              "calibration).")
     parser.add_argument('--err-dmdtda-scaling-factor', type=float, default=0.2,
                         help="scaling factor to account for correlated "
@@ -1136,6 +1148,14 @@ def parse_args(args):
                         help="if --dynamic-spinup is set, define the starting"
                              "year for the simulation. The default is 1979, "
                              "unless the climate data starts later.")
+    parser.add_argument('--dynamic-spinup-periods-to-try', nargs='*',
+                        default=[30, 40, 50, 60, 70, 80, 90, 100],
+                        help="if --dynamic-spinup is set, define additional "
+                             "spinup periods to try, if the spinup starting "
+                             "from --dynamic-spinup-year is not successful. If"
+                             "you do not want to use set"
+                             "'--dynamic-spinup-periods-to-try none' in the"
+                             "terminal.")
     parser.add_argument('--geodetic-mb-file-path', type=str, default=None,
                         help='optional path or URL to a custom geodetic MB '
                              'file passed to MB calibration.')
@@ -1180,6 +1200,9 @@ def parse_args(args):
 
     dynamic_spinup = False if args.dynamic_spinup == '' else args.dynamic_spinup
 
+    if args.dynamic_spinup_periods_to_try == ['none']:
+        args.dynamic_spinup_periods_to_try = None
+
     # All good
     return dict(rgi_version=rgi_version, rgi_reg=rgi_reg,
                 border=border, output_folder=output_folder,
@@ -1211,6 +1234,7 @@ def parse_args(args):
                 dynamic_spinup=dynamic_spinup,
                 err_dmdtda_scaling_factor=args.err_dmdtda_scaling_factor,
                 dynamic_spinup_start_year=args.dynamic_spinup_start_year,
+                dynamic_spinup_periods_to_try=args.dynamic_spinup_periods_to_try,
                 mb_calibration_strategy=args.mb_calibration_strategy,
                 geodetic_mb_file_path=args.geodetic_mb_file_path,
                 store_fl_diagnostics=args.store_fl_diagnostics,

oggm/core/dynamic_spinup.py

@@ -32,8 +32,9 @@ def run_dynamic_spinup(gdir, init_model_filesuffix=None, init_model_yr=None,
                        climate_input_filesuffix='',
                        evolution_model=None,
                        mb_model_historical=None, mb_model_spinup=None,
-                       spinup_period=20, spinup_start_yr=None,
+                       spinup_period_initial=20, spinup_start_yr=None,
                        min_spinup_period=10, spinup_start_yr_max=None,
+                       spinup_periods_to_try=None,
                        target_yr=None, target_value=None,
                        minimise_for='area', precision_percent=1,
                        precision_absolute=1, min_ice_thickness=None,
@@ -81,7 +82,7 @@ def run_dynamic_spinup(gdir, init_model_filesuffix=None, init_model_yr=None,
         historical run. Default is to use a ConstantMassBalance model together
         with the provided parameter climate_input_filesuffix and during the
         period of spinup_start_yr until rgi_year (e.g. 1979 - 2000).
-    spinup_period : int
+    spinup_period_initial : int
         The period how long the spinup should run. Start date of historical run
         is defined "target_yr - spinup_period". Minimum allowed value is 10. If
         the provided climate data starts at year later than
@@ -104,6 +105,10 @@ def run_dynamic_spinup(gdir, init_model_filesuffix=None, init_model_yr=None,
         start from at least. If set, this overrides the min_spinup_period if
         target_yr - spinup_start_yr_max > min_spinup_period.
         Default is None
+    spinup_periods_to_try : list or None
+        If spinup_period_initial and min_spinup_period were not successful, you
+        can provide here a list of spinup periods which should be tried in order.
+        Default is None
     target_yr : int or None
         The year at which we want to match area or volume.
         If None, gdir.rgi_date + 1 is used (the default).
@@ -806,15 +811,21 @@ def get_mismatch(t_spinup):
         spinup_period_initial = min(target_yr - spinup_start_yr,
                                     target_yr - yr_min)
     else:
-        spinup_period_initial = min(spinup_period, target_yr - yr_min)
+        spinup_period_initial = min(spinup_period_initial, target_yr - yr_min)
+
+    # define the spinup periods we try to run in order
     if spinup_period_initial <= min_spinup_period:
-        spinup_periods_to_try = [min_spinup_period]
+        spinup_periods_to_run = [min_spinup_period]
     else:
-        # try out a maximum of three different spinup_periods
-        spinup_periods_to_try = [spinup_period_initial,
-                                 int((spinup_period_initial +
-                                      min_spinup_period) / 2),
+        # add initial and minimum spinup period
+        spinup_periods_to_run = [spinup_period_initial,
+                                 (spinup_period_initial + min_spinup_period) / 2,
                                  min_spinup_period]
+    if spinup_periods_to_try is not None:
+        for spn_prd in spinup_periods_to_try:
+            if min_spinup_period < spn_prd < target_yr - yr_min:
+                spinup_periods_to_run.append(spn_prd)
+
     # after defining the initial spinup period we can define the year for the
     # fixed_geometry_spinup
     if add_fixed_geometry_spinup:
@@ -828,7 +839,7 @@ def get_mismatch(t_spinup):
     if mb_model_spinup is not None:
         provided_mb_model_spinup = True
 
-    for spinup_period in spinup_periods_to_try:
+    for spinup_period in spinup_periods_to_run:
         yr_spinup = target_yr - spinup_period
 
         if not provided_mb_model_spinup:
@@ -852,7 +863,7 @@ def get_mismatch(t_spinup):
         except RuntimeError as e:
             # if the last spinup period was min_spinup_period the dynamic
             # spinup failed
-            if spinup_period == min_spinup_period:
+            if spinup_period == spinup_periods_to_run[-1]:
                 log.warning('No dynamic spinup could be conducted and the '
                             'original model with no spinup is saved using the '
                             f'provided output_filesuffix "{output_filesuffix}". '
@@ -949,13 +960,14 @@ def dynamic_melt_f_run_with_dynamic_spinup(
         gdir, melt_f, yr0_ref_mb, yr1_ref_mb, fls_init, ys, ye,
         output_filesuffix='', evolution_model=None,
         mb_model_historical=None, mb_model_spinup=None,
-        minimise_for='area', climate_input_filesuffix='', spinup_period=20,
+        minimise_for='area', climate_input_filesuffix='', spinup_period_initial=20,
         min_spinup_period=10, target_yr=None, precision_percent=1,
         precision_absolute=1, min_ice_thickness=None,
         first_guess_t_spinup=-2, t_spinup_max_step_length=2, maxiter=30,
         store_model_geometry=True, store_fl_diagnostics=None,
         local_variables=None, set_local_variables=False, do_inversion=True,
         spinup_start_yr_max=None, add_fixed_geometry_spinup=True,
+        spinup_periods_to_try=None,
         **kwargs):
     """
     This function is one option for a 'run_function' for the
@@ -1008,7 +1020,7 @@ def dynamic_melt_f_run_with_dynamic_spinup(
     climate_input_filesuffix : str
         filesuffix for the input climate file
         Default is ''
-    spinup_period : int
+    spinup_period_initial : int
         The period how long the spinup should run. Start date of historical run
         is defined "target_yr - spinup_period". Minimum allowed value is
         defined with 'min_spinup_period'. If the provided climate data starts
@@ -1095,6 +1107,10 @@ def dynamic_melt_f_run_with_dynamic_spinup(
         fixed-geometry-spinup is added at the beginning so that the resulting
         model run always starts from ys.
         Default is True
+    spinup_periods_to_try : list or None
+        If spinup_period_initial and min_spinup_period were not successful, you
+        can provide here a list of spinup periods which should be tried in order.
+        Default is None
     kwargs : dict
         kwargs to pass to the evolution_model instance
 
@@ -1132,7 +1148,7 @@ def dynamic_melt_f_run_with_dynamic_spinup(
                  'adapted and it is the only period which is tried by the '
                  'dynamic spinup function!')
         min_spinup_period = target_yr - ys
-        spinup_period = target_yr - ys
+        spinup_period_initial = target_yr - ys
 
     if spinup_start_yr_max is None:
         spinup_start_yr_max = yr0_ref_mb
@@ -1183,6 +1199,15 @@ def dynamic_melt_f_run_with_dynamic_spinup(
     init_present_time_glacier(gdir, filesuffix=model_flowline_filesuffix,
                               add_to_log_file=False)
 
+    # check if there was already a successful run
+    if 'dynamic_spinup_period' in gdir.get_diagnostics():
+        last_spinup_period = gdir.get_diagnostics()['dynamic_spinup_period']
+        spinup_periods_to_try_forward = [last_spinup_period]
+        if spinup_periods_to_try is not None:
+            spinup_periods_to_try_forward.extend(spinup_periods_to_try)
+    else:
+        spinup_periods_to_try_forward = spinup_periods_to_try
+
     # Now do a dynamic spinup to match area
     # do not ignore errors in dynamic spinup, so all 'bad'/intermediate files
     # are deleted in run_dynamic_spinup function
@@ -1196,10 +1221,11 @@ def dynamic_melt_f_run_with_dynamic_spinup(
             evolution_model=evolution_model,
             mb_model_historical=mb_model_historical,
             mb_model_spinup=mb_model_spinup,
-            spinup_period=spinup_period,
+            spinup_period_initial=spinup_period_initial,
             spinup_start_yr=ys,
             spinup_start_yr_max=spinup_start_yr_max,
             min_spinup_period=min_spinup_period, target_yr=target_yr,
+            spinup_periods_to_try=spinup_periods_to_try_forward,
             precision_percent=precision_percent,
             precision_absolute=precision_absolute,
             min_ice_thickness=min_ice_thickness,
@@ -1240,13 +1266,14 @@ def dynamic_melt_f_run_with_dynamic_spinup_fallback(
         gdir, melt_f, fls_init, ys, ye, local_variables, output_filesuffix='',
         evolution_model=None, minimise_for='area',
         mb_model_historical=None, mb_model_spinup=None,
-        climate_input_filesuffix='', spinup_period=20, min_spinup_period=10,
-        target_yr=None, precision_percent=1,
+        climate_input_filesuffix='', spinup_period_initial=20,
+        min_spinup_period=10, target_yr=None, precision_percent=1,
         precision_absolute=1, min_ice_thickness=None,
         first_guess_t_spinup=-2, t_spinup_max_step_length=2, maxiter=30,
         store_model_geometry=True, store_fl_diagnostics=None,
         do_inversion=True, spinup_start_yr_max=None,
-        add_fixed_geometry_spinup=True, **kwargs):
+        add_fixed_geometry_spinup=True, spinup_periods_to_try=None,
+        **kwargs):
     """
     This is the fallback function corresponding to the function
     'dynamic_melt_f_run_with_dynamic_spinup', which are provided
@@ -1293,7 +1320,7 @@ def dynamic_melt_f_run_with_dynamic_spinup_fallback(
     climate_input_filesuffix : str
         filesuffix for the input climate file
         Default is ''
-    spinup_period : int
+    spinup_period_initial : int
         The period how long the spinup should run. Start date of historical run
         is defined "target_yr - spinup_period". Minimum allowed value is
         defined with 'min_spinup_period'. If the provided climate data starts
@@ -1367,6 +1394,10 @@ def dynamic_melt_f_run_with_dynamic_spinup_fallback(
         fixed-geometry-spinup is added at the beginning so that the resulting
         model run always starts from ys.
         Default is True
+    spinup_periods_to_try : list or None
+        If spinup_period_initial and min_spinup_period were not successful, you
+        can provide here a list of spinup periods which should be tried in order.
+        Default is None
     kwargs : dict
         kwargs to pass to the evolution_model instance
 
@@ -1409,7 +1440,7 @@ def dynamic_melt_f_run_with_dynamic_spinup_fallback(
                  'adapted and it is the only period which is tried by the '
                  'dynamic spinup function!')
         min_spinup_period = target_yr - ys
-        spinup_period = target_yr - ys
+        spinup_period_initial = target_yr - ys
 
     yr_clim_min = gdir.get_climate_info()['baseline_yr_0']
     try:
@@ -1422,10 +1453,11 @@ def dynamic_melt_f_run_with_dynamic_spinup_fallback(
             evolution_model=evolution_model,
             mb_model_historical=mb_model_historical,
             mb_model_spinup=mb_model_spinup,
-            spinup_period=spinup_period,
+            spinup_period_initial=spinup_period_initial,
             spinup_start_yr=ys,
             min_spinup_period=min_spinup_period,
             spinup_start_yr_max=spinup_start_yr_max,
+            spinup_periods_to_try=spinup_periods_to_try,
             target_yr=target_yr,
             minimise_for=minimise_for,
             precision_percent=precision_percent,

oggm/core/flowline.py

@@ -1319,17 +1319,17 @@ def run_until_and_store(self, y1,
                             ds['volume_bsl_m3'].data[i, :] = fl.volume_bsl_m3
                         if 'volume_bwl' in ovars_fl:
                             ds['volume_bwl_m3'].data[i, :] = fl.volume_bwl_m3
-                        if 'ice_velocity' in ovars_fl and (yr > self.y0):
+                        if 'ice_velocity' in ovars_fl and (yr > self.y0) and (i != 0):
                             # Velocity can only be computed with dynamics
                             var = self.u_stag[fl_id]
                             val = (var[1:fl.nx + 1] + var[:fl.nx]) / 2 * self._surf_vel_fac
                             ds['ice_velocity_myr'].data[i, :] = val * cfg.SEC_IN_YEAR
-                        if 'dhdt' in ovars_fl and (yr > self.y0):
+                        if 'dhdt' in ovars_fl and (yr > self.y0) and (i != 0):
                             # dhdt can only be computed after one step
                             val = fl.thick - thickness_previous_dhdt[fl_id]
                             ds['dhdt'].data[i, :] = val
                             thickness_previous_dhdt[fl_id] = fl.thick
-                        if 'climatic_mb' in ovars_fl and (yr > self.y0):
+                        if 'climatic_mb' in ovars_fl and (yr > self.y0) and (i != 0):
                             # yr - 1 to use the mb which lead to the current
                             # state, also using previous surface height
                             val = self.get_mb(surface_h_previous[fl_id],
@@ -1345,7 +1345,7 @@ def run_until_and_store(self, y1,
                                 0.,
                                 val * fac_sec)
                             surface_h_previous[fl_id] = fl.surface_h
-                        if 'flux_divergence' in ovars_fl and (yr > self.y0):
+                        if 'flux_divergence' in ovars_fl and (yr > self.y0) and (i != 0):
                             # calculated after the formula dhdt = mb + flux_div
                             val = ds['dhdt'].data[i, :] - ds['climatic_mb'].data[i, :]
                             # special treatment for retreating: If the glacier