Add mb params perturbation utility func

oggm/cfg.py

@@ -312,14 +312,16 @@ def set_logging_config(logging_level='INFO'):
         Print confirmation that things are working as expected, e.g. when
         each task is run correctly (this is the default).
     WARNING
-        Indication that something unexpected happened on a glacier,
-        but that OGGM is still working on this glacier.
+        Do not print INFO or DEBUG but print WARNING (which indicate that
+        something unexpected happened on a glacier but that OGGM is
+        still working on this glacier).
     ERROR
         Print workflow messages and errors only, e.g. when a glacier cannot
         run properly.
     WORKFLOW
         Print only high level, workflow information (typically, one message
-        per task). Errors and warnings will NOT be printed.
+        per task). Errors and warnings will NOT be printed. This is the level
+        we recommend for operational large-scale runs.
     CRITICAL
         Print nothing but fatal errors.
 

oggm/core/flowline.py

@@ -1328,8 +1328,7 @@ def run_until_and_store(self, y1,
                             surface_h_previous[fl_id] = fl.surface_h
                         if 'flux_divergence' in ovars_fl and (yr > self.y0):
                             # calculated after the formula dhdt = mb + flux_div
-                            val = ds['dhdt_myr'].data[j, :] -\
-                                  ds['climatic_mb_myr'].data[j, :]
+                            val = ds['dhdt_myr'].data[j, :] - ds['climatic_mb_myr'].data[j, :]
                             # special treatment for retreating: If the glacier
                             # terminus is getting ice free it means the
                             # climatic mass balance is more negative than the
@@ -3173,6 +3172,7 @@ def decide_evolution_model(evolution_model=None):
 def flowline_model_run(gdir, output_filesuffix=None, mb_model=None,
                        ys=None, ye=None, zero_initial_glacier=False,
                        init_model_fls=None, store_monthly_step=False,
+                       glen_a_fac=None, fs_fac=None,
                        fixed_geometry_spinup_yr=None,
                        store_model_geometry=None,
                        store_fl_diagnostics=None,
@@ -3209,6 +3209,13 @@ def flowline_model_run(gdir, output_filesuffix=None, mb_model=None,
     store_monthly_step : bool
         whether to store the diagnostic data at a monthly time step or not
         (default is yearly)
+    glen_a_fac : float
+        perturbate the default glen_a (either from the inversion if
+        PARAMS['use_inversion_params_for_run'] is True or PARAMS['glen_a'] if
+        False) with a chosen multiplicative parameter. Useful for
+        parameter perturbation experiments.
+    fs_fac : float
+        same as glen_a_fac but for sliding
     store_model_geometry : bool
         whether to store the full model geometry run file to disk or not.
         (new in OGGM v1.4.1: default is to follow
@@ -3270,6 +3277,11 @@ def flowline_model_run(gdir, output_filesuffix=None, mb_model=None,
         fs = cfg.PARAMS['fs']
         glen_a = cfg.PARAMS['glen_a']
 
+    if glen_a_fac is not None:
+        glen_a *= glen_a_fac
+    if fs_fac is not None:
+        fs *= fs_fac
+
     kwargs.setdefault('fs', fs)
     kwargs.setdefault('glen_a', glen_a)
 
@@ -3320,11 +3332,11 @@ def flowline_model_run(gdir, output_filesuffix=None, mb_model=None,
                                        'to prevent this error.')
 
     model = evolution_model(fls, mb_model=mb_model, y0=ys,
-                           inplace=True,
-                           is_tidewater=gdir.is_tidewater,
-                           is_lake_terminating=gdir.is_lake_terminating,
-                           water_level=water_level,
-                           **kwargs)
+                            inplace=True,
+                            is_tidewater=gdir.is_tidewater,
+                            is_lake_terminating=gdir.is_lake_terminating,
+                            water_level=water_level,
+                            **kwargs)
 
     with warnings.catch_warnings():
         # For operational runs we ignore the warnings
@@ -3429,7 +3441,7 @@ def run_random_climate(gdir, nyears=1000, y0=None, halfsize=15,
         if false, every model year will be chosen from the random climate
         period with the same probability
     kwargs : dict
-        kwargs to pass to the FluxBasedModel instance
+        kwargs to pass to the flowline_model_run task
     """
 
     mb_model = MultipleFlowlineMassBalance(gdir,
@@ -3546,7 +3558,7 @@ def run_constant_climate(gdir, nyears=1000, y0=None, halfsize=15,
         list of flowlines to use to initialize the model (the default is the
         present_time_glacier file from the glacier directory)
     kwargs : dict
-        kwargs to pass to the FluxBasedModel instance
+        kwargs to pass to the flowline_model_run task
     """
 
     mb_model = MultipleFlowlineMassBalance(gdir,
@@ -3661,14 +3673,14 @@ def run_from_climate_data(gdir, ys=None, ye=None, min_ys=None, max_ys=None,
         multiply a factor to the precipitation time series (note that
         this factor is multiplied to any factor that was decided during
         calibration or by global parameters)
-    kwargs : dict
-        kwargs to pass to the FluxBasedModel instance
     fixed_geometry_spinup_yr : int
         if set to an integer, the model will artificially prolongate
         all outputs of run_until_and_store to encompass all time stamps
         starting from the chosen year. The only output affected are the
         glacier wide diagnostic files - all other outputs are set
         to constants during "spinup"
+    kwargs : dict
+        kwargs to pass to the flowline_model_run task
     """
 
     if init_model_filesuffix is not None:

oggm/core/massbalance.py

@@ -320,6 +320,7 @@ class MonthlyTIModel(MassBalanceModel):
     def __init__(self, gdir,
                  filename='climate_historical',
                  input_filesuffix='',
+                 mb_params_filesuffix='',
                  fl_id=None,
                  melt_f=None,
                  temp_bias=None,
@@ -340,7 +341,10 @@ def __init__(self, gdir,
             set to a different BASENAME if you want to use alternative climate
             data. Default is 'climate_historical'
         input_filesuffix : str, optional
-            append a suffix to the filename (useful for GCM runs).
+            append a suffix to the climate input filename (useful for GCM runs).
+        mb_params_filesuffix : str, optional
+            append a suffix to the mb params calibration file (useful for
+            sensitivity runs).
         fl_id : int, optional
             if this flowline has been calibrated alone and has specific
             model parameters.
@@ -378,6 +382,7 @@ def __init__(self, gdir,
         super(MonthlyTIModel, self).__init__()
         self.valid_bounds = [-1e4, 2e4]  # in m
         self.fl_id = fl_id  # which flowline are we the model of?
+        self._mb_params_filesuffix = mb_params_filesuffix  # which mb params?
         self.gdir = gdir
 
         if melt_f is None:
@@ -543,13 +548,16 @@ def temp_bias(self, new_temp_bias):
     @lazy_property
     def calib_params(self):
         if self.fl_id is None:
-            return self.gdir.read_json('mb_calib')
+            return self.gdir.read_json('mb_calib', self._mb_params_filesuffix)
         else:
             try:
-                return self.gdir.read_json('mb_calib',
-                                           filesuffix=f'_fl{self.fl_id}')
+                out = self.gdir.read_json('mb_calib', filesuffix=f'_fl{self.fl_id}')
+                if self._mb_params_filesuffix:
+                    raise InvalidWorkflowError('mb_params_filesuffix cannot be '
+                                               'used with multiple flowlines')
+                return out
             except FileNotFoundError:
-                return self.gdir.read_json('mb_calib')
+                return self.gdir.read_json('mb_calib', self._mb_params_filesuffix)
 
     def is_year_valid(self, year):
         return self.ys <= year <= self.ye
@@ -1403,6 +1411,7 @@ def mb_calibration_from_wgms_mb(gdir, **kwargs):
     Parameters
     ----------
     **kwargs : any kwarg accepted by mb_calibration_from_scalar_mb
+    except `ref_mb` and `ref_mb_years`
     """
 
     # Note that this currently does not work for hydro years (WGMS uses hydro)
@@ -1428,12 +1437,15 @@ def mb_calibration_from_geodetic_mb(gdir, *,
                                     calibrate_param2=None,
                                     calibrate_param3=None,
                                     mb_model_class=MonthlyTIModel,
+                                    filesuffix='',
                                     ):
     """Calibrate for geodetic MB data from Hugonnet et al., 2021.
 
     The data table can be obtained with utils.get_geodetic_mb_dataframe().
     It is equivalent to the original data from Hugonnet, but has some outlier
-    values filtered. See `this notebook` for more details.
+    values filtered. See [this notebook](https://nbviewer.org/urls/
+    cluster.klima.uni-bremen.de/~oggm/geodetic_ref_mb/convert_vold1.ipynb)
+    for more details.
 
     The problem of calibrating many unknown parameters on geodetic data is
     currently unsolved. This is OGGM's current take, based on trial and
@@ -1480,6 +1492,10 @@ def mb_calibration_from_geodetic_mb(gdir, *,
         the MassBalanceModel to use for the calibration. Needs to use the
         same parameters as MonthlyTIModel (the default): melt_f,
         temp_bias, prcp_fac.
+    filesuffix: str
+        add a filesuffix to mb_params.json. This could be useful for sensitivity
+        analyses with MB models, if they need to fetch other sets of params for
+        example.
 
     Returns
     -------
@@ -1551,6 +1567,7 @@ def mb_calibration_from_geodetic_mb(gdir, *,
                                              prcp_fac_max=prcp_fac_max,
                                              temp_bias=temp_bias,
                                              mb_model_class=mb_model_class,
+                                             filesuffix=filesuffix,
                                              )
 
     else:
@@ -1565,6 +1582,7 @@ def mb_calibration_from_geodetic_mb(gdir, *,
                                              calibrate_param3=calibrate_param3,
                                              temp_bias=temp_bias,
                                              mb_model_class=mb_model_class,
+                                             filesuffix=filesuffix,
                                              )
 
 
@@ -1895,6 +1913,68 @@ def to_minimize(x, model_attr):
     return df
 
 
+@entity_task(log, writes=['mb_calib'])
+def perturbate_mb_params(gdir, perturbation=None, reset_default=False, filesuffix=''):
+    """Replaces pre-calibrated MB params with perturbed ones for this glacier.
+
+    It simply replaces the existing `mb_calib.json` file with an
+    updated one with perturbed parameters. The original ones
+    are stored in the file for re-use after perturbation.
+
+    Users can change the following 4 parameters:
+    - 'melt_f': unit [kg m-2 day-1 K-1], the melt factor
+    - 'prcp_fac': unit [-], the precipitation factor
+    - 'temp_bias': unit [K], the temperature correction applied to the timeseries
+    - 'bias': unit [mm we yr-1], *substracted* from the computed MB. Rarely used.
+
+    All parameter perturbations are additive, i.e. the value
+    provided by the user is added to the *precalibrated* value.
+    For example, `temp_bias=1` means that the temp_bias used by the
+    model will be the precalibrated one, plus 1 Kelvin.
+
+    The only exception is prpc_fac, which is multiplicative.
+    For example prcp_fac=1 will leave the precalibrated prcp_fac unchanged,
+    while 2 will double it.
+
+    Parameters
+    ----------
+    perturbation : dict
+        the parameters to change and the associated value (see doc above)
+    reset_default : bool
+        reset the parameters to their original value. This might be
+        unnecessary if using the filesuffix mechanism.
+    filesuffix : str
+        write the modified parameters in a separate mb_calib.json file
+        with the filesuffix appended. This can then be read by the
+        MassBalanceModel for example instead of the default one.
+        Note that it's always the default, precalibrated params
+        file which is read to start with.
+    """
+    df = gdir.read_json('mb_calib')
+
+    # Save original params if not there
+    if 'bias_orig' not in df:
+        for k in ['bias', 'melt_f', 'prcp_fac', 'temp_bias']:
+            df[k + '_orig'] = df[k]
+
+    if reset_default:
+        for k in ['bias', 'melt_f', 'prcp_fac', 'temp_bias']:
+            df[k] = df[k + '_orig']
+        gdir.write_json(df, 'mb_calib', filesuffix=filesuffix)
+        return df
+
+    for k, v in perturbation.items():
+        if k == 'prcp_fac':
+            df[k] = df[k + '_orig'] * v
+        elif k in ['bias', 'melt_f', 'temp_bias']:
+            df[k] = df[k + '_orig'] + v
+        else:
+            raise InvalidParamsError(f'Perturbation not valid: {k}')
+
+    gdir.write_json(df, 'mb_calib', filesuffix=filesuffix)
+    return df
+
+
 def _check_terminus_mass_flux(gdir, fls):
     # Check that we have done this correctly
     rho = cfg.PARAMS['ice_density']

oggm/tasks.py

@@ -32,6 +32,7 @@
 from oggm.core.massbalance import mb_calibration_from_wgms_mb
 from oggm.core.massbalance import apparent_mb_from_linear_mb
 from oggm.core.massbalance import apparent_mb_from_any_mb
+from oggm.core.massbalance import perturbate_mb_params
 from oggm.core.massbalance import fixed_geometry_mass_balance
 from oggm.core.massbalance import compute_ela
 from oggm.shop.w5e5 import process_w5e5_data