Sermeq 1 ks

docs/install_pygem.md

@@ -6,22 +6,23 @@ The model is stored in two repositories [(see model structure)](model_structure_
 A conda environment is a directory that contains a specific collection of installed packages. The use of environments reduces issues caused by package dependencies. The model is designed to be compatible with OGGM. We therefore get started by following the [installation instructions from OGGM](https://docs.oggm.org/en/stable/installing-oggm.html).
 
 Once your conda environment is setup for OGGM, add the core of PyGEM using pip.
-
 ```
 pip install pygem
 ```
 
 This will provide you with a conda environment that has the basic functionality to run simple calibration options (e.g., 'HH2015', 'HH2015mod') and simulations. If you want to use the emulators and Bayesian inference, the advanced environment is required.
 
 ### Developing PyGEM
-Are you interested in developing PyGEM? If so, we recommend uninstalling pygem, forking the [PyGEM's github repository](https://github.com/drounce/PyGEM) and then [cloning the github repository](https://docs.github.com/en/repositories/creating-and-managing-repositories/cloning-a-repository) onto your local machine.
+Are you interested in developing PyGEM? If so, we recommend forking the [PyGEM's github repository](https://github.com/PyGEM-Community/PyGEM) and then [cloning the github repository](https://docs.github.com/en/repositories/creating-and-managing-repositories/cloning-a-repository) onto your local machine.
 
+Note, if PyGEM was already installed via PyPI, first uninstall:
 ```
 pip uninstall pygem
-git clone https://github.com/drounce/PyGEM.git
+````
+
+You can then use pip to install your locally cloned fork of PyGEM in 'editable' mode like so:
 ```
-```{warning}
-The directory structure here is important. The cloned <em>PyGEM</em> directory should be on the same level as the <em>PyGEM-scripts</em>.
+pip install -e /path/to/your/PyGEM/clone
 ```
 
 ### Advanced environment: GPyTorch (emulator only)
@@ -64,7 +65,7 @@ The only way to find out if your package dependencies work is to test it by runn
 
 
 ## Install PyGEM-Scripts
-The scripts that are used to run PyGEM are located in the [PyGEM-Scripts repository](https://github.com/drounce/PyGEM-scripts) on github. To run the model, you can either (i) clone the repository or (ii) fork the repository to develop/add your own scripts. For instructions, follow github’s instructions on [cloning](https://docs.github.com/en/repositories/creating-and-managing-repositories/cloning-a-repository) or [forking a repository](https://docs.github.com/en/get-started/quickstart/fork-a-repo). Once the repository is installed on your local machine, you can run the model from this directory.
+The scripts that are used to run PyGEM are located in the [PyGEM-Scripts repository](https://github.com/PyGEM-Community/PyGEM-scripts) on github. To run the model, you can either (i) clone the repository or (ii) fork the repository to develop/add your own scripts. For instructions, follow github’s instructions on [cloning](https://docs.github.com/en/repositories/creating-and-managing-repositories/cloning-a-repository) or [forking a repository](https://docs.github.com/en/get-started/quickstart/fork-a-repo). Once the repository is installed on your local machine, you can run the model from this directory.
 
 ```{note}
 Be sure that your [directory structure](directory_structure_target) is setup properly before you try running the model!

pygem/gcmbiasadj.py

@@ -321,8 +321,9 @@ def prec_biasadj_opt1(ref_prec, ref_elev, gcm_prec, dates_table_ref, dates_table
     gcm_prec_biasadj_subset = (
             gcm_prec_biasadj[:,gcm_subset_idx_start:gcm_subset_idx_end+1][:,gcm_spinupyears*12:])
     gcm_prec_biasadj_frac = gcm_prec_biasadj_subset.sum(axis=1) / ref_prec_nospinup.sum(axis=1)
-    assert np.min(gcm_prec_biasadj_frac) > 0.5 and np.max(gcm_prec_biasadj_frac) < 2, (
-            'Error with gcm precipitation bias adjustment: total ref and gcm prec differ by more than factor of 2')
+    assert np.min(gcm_prec_biasadj_frac) > 0.45 and np.max(gcm_prec_biasadj_frac) < 2, (
+            f'Error with gcm precipitation bias adjustment: total ref and gcm prec differ by more than factor of 2.'
+            f' min: {np.min(gcm_prec_biasadj_frac)}, max: {np.max(gcm_prec_biasadj_frac)}')
     assert gcm_prec_biasadj.max() <= 10, 'gcm_prec_adj (precipitation bias adjustment) too high, needs to be modified'
     assert gcm_prec_biasadj.min() >= 0, 'gcm_prec_adj is producing a negative precipitation value'   
 

pygem/glacierdynamics.py

@@ -21,6 +21,303 @@
 
 cfg.initialize()
 
+# calving law
+    #%% ----- Calving Law ----------
+    #%% ----- Calving Law ----------
+def fa_sermeq_speed_law(model,last_above_wl, v_scaling=1, verbose=False,
+                    tau0=1.5, yield_type='constant', mu=0.01,
+                    trim_profile=0):
+    """
+    This function is used to calculate frontal ablation given ice speed forcing,
+    for lake-terminating and tidewater glaciers
+
+    @author: Ruitang Yang & Lizz Ultee
+
+    Authors: Ruitang Yang & Lizz Ultee
+    Parameters
+    ----------
+
+    model : oggm.core.flowline.FlowlineModel
+        the model instance calling the function
+    flowline : oggm.core.flowline.Flowline
+        the instance of the flowline object on which the calving law is called
+    fl_id : float, optional
+        the index of the flowline in the fls array (might be ignored by some MB models)
+    last_above_wl : int
+        the index of the last pixel above water (in case you need to know
+        where it is).
+    v_scaling: float
+        velocity scaling factor, >0, default is 1
+    Terminus_mb : array
+        Mass balance along the flowline or nearest the terminus [m/a]. Default None...
+        TODO: set default behavior, check the unit meter of ice per year or m w.e. per year?
+    verbose: Boolean, optional
+        Whether to print component parts for inspection.  Default False.
+
+    tau0: float, optional
+        This glacier's yield strength [Pa]. Default is 150 kPa.
+    yield_type: str, optional
+        'constant' or 'variable' (Mohr-Coulomb) yielding. Default is constant.
+    mu: float, optional
+        Mohr-Coulomb cohesion, a coefficient between 0 and 1. Default is 0.01.
+        Only used if we have variable yield
+
+    trim_profile: int, optional
+        How many grid cells at the end of the profile to ignore.  Default is 1.
+        If the initial profile is set by k-calving (as in testing) there can be a
+        weird cliff shape with very thin final grid point and large velocity gradient
+
+    Returns
+    -------
+    fa_viscoplastic: float
+        Frontal ablation rate [m/a] based on viscoplastic assumptions
+    SQFA: dict
+        Frontal ablation rate [m/a] based on viscoplastic assumptions
+        serface elevation at the terminus [m a.s.l.] based on OGGM ## TODO CHANGE IT BY PyGEM Surface mass balance result
+        bed elevation at the terminus [m a.s.l.]
+        Terminus Thickness [m]
+        Yield terminus thickness [m]
+        Velocity at the terminus [m/a]
+        Surface mass balance at the terminus [m/a]
+        Length change at the terminus [m/a] based on viscoplastic assumptions
+        TODO: output the length change in case we
+        have the length change results from observations in-situ or remote sensing (Thomas Schellenberger has the machine learning products)
+        Frontal ablation rate [m/a] based on viscoplastic assumptions
+
+    Explanation of sign
+    -------
+    fa_viscoplastic: negative, mass loss
+    dLdt: length change rate, positive if advance; negative if retreat
+    terminus mass balance: negative if mass loss; positive if mass gain
+    """
+    # ---------------------------------------------------------------------------
+    # class NegativeValueError(Exception):
+    #     pass
+    # ---------------------------------------------------------------------------
+    ## Global constants
+    G = 9.8  # acceleration due to gravity in m/s^2
+    RHO_ICE = 920.0  # ice density kg/m^3
+    RHO_SEA = 1020.0  # seawater density kg/m^3
+
+    # ---------------------------------------------------------------------------
+    # the yield strength
+    def tau_y(tau0=1.5, yield_type='constant', bed_elev=None, thick=None, mu=0.01):
+        """
+        Functional form of yield strength.
+        Can do constant or Mohr-Coulomb yield strength.  Ideally, the glacier's yield type
+        ('constant' or 'variable') would be saved in a model instance.
+
+        Parameters
+        ----------
+        tau0: float, optional
+            Initial guess for yield strength [Pa]. Default is 150 kPa.
+        yield_type: str, optional
+            'constant' or 'variable' (Mohr-Coulomb) yielding. Default is constant.
+        bed_elev: float, optional
+            Bed elevation, dimensional [m]. The default is None.
+        thick: float, optional
+            Ice thickness, dimensional [m]. The default is None.
+        mu: float, optional
+            Mohr-Coulomb cohesion, a coefficient between 0 and 1. Default is 0.01.
+
+        Returns
+        -------
+        tau_y: float
+            The yield strength for these conditions.
+        """
+        tau1=tau0*1e5
+        if yield_type == 'variable':
+            try:
+                if bed_elev < 0:
+                    D = -1 * bed_elev  # Water depth D the nondim bed topography value when Z<0
+                else:
+                    D = 0
+            except:
+                print('You must set a bed elevation and ice thickness to use variable yield strength.')
+            N = RHO_ICE * G * thick - RHO_SEA * G * D  # Normal stress at bed
+            ty = tau1 + mu * N
+        else:  # assume constant if not set
+            ty = tau1
+        return ty
+
+
+    # ---------------------------------------------------------------------------
+    # calculate the yield ice thickness
+
+    def balance_thickness(yield_strength, bed_elev):
+        """
+        Ice thickness such that the stress matches the yield strength.
+
+        Parameters
+        ----------
+        yield_strength: float
+            The yield strength near the terminus.
+            If yield type is constant, this will of course be the same everywhere.  If yield type is
+            variable (Mohr-Coulomb), the yield strength at the terminus could differ from elsewhere.
+        bed_elev: float
+            Elevation of glacier bed at the terminus
+
+        Returns
+        -------
+        Hy: float
+            The ice thickness for stress balance at the terminus.
+        """
+        if bed_elev < 0:
+            D = -1 * bed_elev
+        else:
+            D = 0
+        return (2 * yield_strength / (RHO_ICE * G)) + np.sqrt(
+            (RHO_SEA * (D ** 2) / RHO_ICE) + ((2 * yield_strength / (RHO_ICE * G)) ** 2))
+        # TODO: Check on exponent on last term.  In Ultee & Bassis 2016, this is squared, but in Ultee & Bassis 2020 supplement, it isn't.
+
+    # ---------------------------------------------------------------------------
+    # calculate frontal ablation based on the ice thickness, speed at the terminus
+    fls=model.fls
+    flowline=fls[-1]
+    surface_m = flowline.surface_h
+    bed_m = flowline.bed_h
+    width_m = flowline.widths_m
+    velocity_m = model.u_stag[-1]*cfg.SEC_IN_YEAR
+    x_m = flowline.dis_on_line*flowline.map_dx/1000
+
+    # gdir : py:class:`oggm.GlacierDirectory`
+    #     the glacier directory to process
+    # fls = model.gdir.read_pickle('model_flowlines')
+    # mbmod_fl = massbalance.MultipleFlowlineMassBalance(model.gdir, fls=fls, use_inversion_flowlines=True,
+    #                                                    mb_model_class=MonthlyTIModel)
+    mb_annual=model.mb_model.get_annual_mb(heights=surface_m, fl_id=-1, year=model.yr, fls=model.fls)
+
+    Terminus_mb = mb_annual*cfg.SEC_IN_YEAR
+    # slice up to index+1 to include the last nonzero value
+    # profile: NDarray
+    #     The current profile (x, surface, bed,width) as calculated by the base model
+    #     Unlike core SERMeQ, these should be DIMENSIONAL [m].
+    profile=(x_m[:last_above_wl+1],
+                surface_m[:last_above_wl+1],
+                bed_m[:last_above_wl+1],width_m[:last_above_wl+1])
+    # model_velocity: array
+    #     Velocity along the flowline [m/a] as calculated by the base model
+    #     Should have values for the points nearest the terminus...otherwise
+    #     doesn't matter if this is the same shape as the profile array.
+    #     TODO: Check with the remote sensing products, or at least to validate the model products
+    model_velocity=velocity_m[:last_above_wl+1]
+    # remove lowest cells if needed
+    last_index = -1 * (trim_profile + 1)
+    ## TODO: Check the flowline model, the decrease the distance between two adjacent points along the flowline, and then calculate the averaged gradient for dhdx,dhydx,dudx
+    ##
+    if isinstance(Terminus_mb, (int, float)):
+        terminus_mb = Terminus_mb
+    elif isinstance(Terminus_mb, (list, np.ndarray)):
+        terminus_mb = Terminus_mb[last_index]
+    else:
+        print("please input the correct mass balance datatype")
+    #
+    if isinstance(model_velocity, (int, float)):
+        model_velocity = v_scaling * model_velocity
+    elif isinstance(model_velocity, list):
+        model_velocity = v_scaling * np.array(model_velocity)
+    elif isinstance(model_velocity, np.ndarray):
+        model_velocity = v_scaling * model_velocity
+    else:
+        print("please input the correct velocity datatype")
+    ## Ice thickness and yield thickness nearest the terminus
+    se_terminus = profile[1][last_index]
+    bed_terminus = profile[2][last_index]
+    h_terminus = se_terminus - bed_terminus
+    width_terminus = profile[3][last_index]
+    tau_y_terminus = tau_y(tau0=tau0, bed_elev=bed_terminus, thick=h_terminus, yield_type=yield_type)
+    Hy_terminus = balance_thickness(yield_strength=tau_y_terminus, bed_elev=bed_terminus)
+    if isinstance(model_velocity, (int, float)):
+        U_terminus = model_velocity
+        U_adj = model_velocity
+    else:
+        U_terminus = model_velocity[last_index]  ## velocity, assuming last point is terminus
+        U_adj = model_velocity[last_index - 1]
+    ## Ice thickness and yield thickness at adjacent point
+    se_adj = profile[1][last_index - 1]
+    bed_adj = profile[2][last_index - 1]
+    H_adj = se_adj - bed_adj
+    tau_y_adj = tau_y(tau0=tau0, bed_elev=bed_adj, thick=H_adj, yield_type=yield_type)
+    Hy_adj = balance_thickness(yield_strength=tau_y_adj, bed_elev=bed_adj)
+    # Gradients
+    dx_term = profile[0][last_index] - profile[0][last_index - 1]  ## check grid spacing close to terminus
+    dHdx = (h_terminus - H_adj) / dx_term
+    dHydx = (Hy_terminus - Hy_adj) / dx_term
+    if np.isnan(U_terminus) or np.isnan(U_adj):
+        dUdx = np.nan  ## velocity gradient
+        ## Group the terms
+        dLdt_numerator = np.nan
+        dLdt_denominator = np.nan  ## TODO: compute dHydx
+        dLdt_viscoplastic = np.nan
+        # fa_viscoplastic = dLdt_viscoplastic -U_terminus  ## frontal ablation rate
+        fa_viscoplastic = np.nan  ## frontal ablation rate
+    else:
+        # Gradients
+        # dx_term = profile[0][last_index] - profile[0][last_index - 1]  ## check grid spacing close to terminus
+        # dHdx = (h_terminus - H_adj) / dx_term
+        # dHydx = (Hy_terminus - Hy_adj) / dx_term
+        dUdx = (U_terminus - U_adj) / dx_term  ## velocity gradient
+        ## Group the terms
+        dLdt_numerator = terminus_mb - (h_terminus * dUdx) - (U_terminus * dHdx)
+        dLdt_denominator = dHydx - dHdx  ## TODO: compute dHydx
+        dLdt_viscoplastic = dLdt_numerator / dLdt_denominator
+        # fa_viscoplastic = dLdt_viscoplastic -U_terminus  ## frontal ablation rate
+
+        # try:
+        U_calving = U_terminus - dLdt_viscoplastic  ## frontal ablation rate
+        fa_viscoplastic=U_calving
+        # if U_calving<0:
+        #     print("The glacier is advancing, and the advancing rate is larger than ice flow speed at the terminus, please check ")
+        #     if U_calving>0 or U_calving==0:
+        #         fa_viscoplastic=U_calving
+        #     else:
+        #         fa_viscoplastic=U_calving
+        #         # fa_viscoplastic=np.nan
+        #         raise NegativeValueError("Something is wrong, right now the calving in negative, which should be positive or zero")
+        # except NegativeValueError as e:
+        #     print ("The glacier is advancing, and the advancing rate is larger than ice flow speed at the terminus, please check ")
+
+
+    SQFA = {'se_terminus': se_terminus,
+            'bed_terminus': bed_terminus,
+            'Thickness_termi': h_terminus,
+            'Width_termi':  width_terminus,
+            'Hy_thickness': Hy_terminus,
+            'Velocity_termi': U_terminus,
+            'Terminus_mb': terminus_mb,
+            'dLdt': dLdt_viscoplastic,
+            'Sermeq_fa': fa_viscoplastic}
+    if verbose:
+        print('For inspection on debugging - all should be DIMENSIONAL (m/a):')
+        #         print('profile_length={}'.format(profile_length))
+        print('last_index={}'.format(last_index))
+        print('se_terminus={}'.format(se_terminus))
+        print('bed_terminus={}'.format(bed_terminus))
+        print('se_adj={}'.format(se_adj))
+        print('bed_adj={}'.format(bed_adj))
+        print('Thicknesses: Hterm {}, Hadj {}'.format(h_terminus, H_adj))
+        print('Hy_terminus={}'.format(Hy_terminus))
+        print('Hy_adj={}'.format(Hy_adj))
+        print('U_terminus={}'.format(U_terminus))
+        print('U_adj={}'.format(U_adj))
+        print('dUdx={}'.format(dUdx))
+        print('dx_term={}'.format(dx_term))
+        print('Checking dLdt: terminus_mb = {}. \n H dUdx = {}. \n U dHdx = {}.'.format(terminus_mb, dUdx * h_terminus,
+                                                                                        U_terminus * dHdx))
+        print('Denom: dHydx = {} \n dHdx = {}'.format(dHydx, dHdx))
+        print('Viscoplastic dLdt={}'.format(dLdt_viscoplastic))
+        print('Terminus surface mass balance ma= {}'.format(terminus_mb))
+        print('Sermeq frontal ablation ma={}'.format(fa_viscoplastic))
+    else:
+        pass
+    return SQFA
+
+
+
+
+
+
+
 
 #%%
 class MassRedistributionCurveModel(FlowlineModel):
@@ -455,8 +752,11 @@ def updategeometry(self, year, debug=False):
             self.mb_model.glac_bin_width_annual[:,year+1] = fl.widths_m
             self.mb_model.glac_wide_area_annual[year+1] = glacier_area.sum()
             self.mb_model.glac_wide_volume_annual[year+1] = (fl.section * fl.dx_meter).sum()
-
-
+
+
+
+
+
     #%% ----- FRONTAL ABLATION -----
     def _get_annual_frontalablation(self, heights, year=None, fls=None, fl_id=None, calving_k=None, debug=False
                                     ):
@@ -520,7 +820,13 @@ def _get_annual_frontalablation(self, heights, year=None, fls=None, fl_id=None,
                     h = fl.thick[last_above_wl]
                     d = h - (fl.surface_h[last_above_wl] - self.water_level)
                     k = self.calving_k
-                    q_calving = k * d * h * fl.widths_m[last_above_wl]
+                    #q_calving = k * d * h * fl.widths_m[last_above_wl]
+                    # Sermeq calving law
+                    print('****************Sermeq claving law start********************')
+                    s_fa = fa_sermeq_speed_law(self,last_above_wl=last_above_wl,v_scaling = 1, verbose = False,tau0 = self.calving_k,
+                                        yield_type = 'constant', mu = 0.01,trim_profile = 0)
+                    q_calving = s_fa ['Sermeq_fa']*s_fa['Thickness_termi']*s_fa['Width_termi']/cfg.SEC_IN_YEAR
+                    print('****************Sermeq claving law end successfully********************')
 
                     # Max frontal ablation is removing all bins with bed below water level
                     glac_idx_bsl = np.where((fl.thick > 0) & (fl.bed_h < self.water_level))[0]