Glens A, warnings and outputs (#152)

* Warnings for bad parameter choices.

This adds warnings when the user
- Sets a too low or high creep factor (Glen A)
- Sets a too low or high basal sliding.
- If the ELA is either above or below the altitude extent of the
  glacier.
- If spatial resolution is below 50m.

* Glens A temperature conversion.

This makes it possible to set the creep factor (Glens A) by supplying a
temperature below zero. The creep setter interprets the value and
converts it. If the value is already a creep factor (>0), it it set
directly.

* Accumulation area ratio

Adds the new attribute accumulation area ratio, as calculated in the
GlacierExplorer app. Also adds it to the summary table.

* Update tests for new bed angles. We can now set bed angles between -90
and 90 degrees.

oggm_edu/funcs.py

@@ -1,6 +1,7 @@
 """oggm-edu package: useful functions diffult to place elsewhere"""
 import urllib
 from functools import wraps
+import numpy as np
 import matplotlib as mpl
 import matplotlib.pyplot as plt
 import seaborn as sns
@@ -130,3 +131,36 @@ def expression_parser(expression, numeric):
         }
 
         return expression_dict[operator]
+
+
+def cp_glen_a(t):
+    """This implements the conversion between temperature and Glen's A,
+       as eq. 3.35 in  The physics of glaciers by Cuffey and Paterson.
+
+    Parameters
+    ----------
+    t : int, float
+        Temperature in degrees Celsius. Should be below 0.
+    """
+
+    # Constants.
+    a_star = 3.5 * 1e-25
+    q_plus = 115_000
+    t_star = 263 + 7 * 1e-8
+    q_minus = 60_000
+    r = 8.314
+
+    # Check if temp is below 0.
+    if t > 0:
+        raise ValueError("Supplied temperature should be below 0.")
+
+    # Convert to Kelvin
+    t = t + 273
+    t_h = t + 7 * 1e-8
+
+    if t < t_star:
+        q_c = q_minus
+    else:
+        q_c = q_plus
+
+    return a_star * np.exp(-(q_c / r) * (1 / t_h - 1 / t_star))

oggm_edu/glacier.py

@@ -10,7 +10,7 @@
 # Internals
 from oggm_edu.glacier_bed import GlacierBed
 from oggm_edu.mass_balance import MassBalance
-from oggm_edu.funcs import edu_plotter
+from oggm_edu.funcs import edu_plotter, cp_glen_a
 
 # Other libraries
 import numpy as np
@@ -91,6 +91,15 @@ def __init__(self, bed, mass_balance):
 
         self._mass_balance = copy.deepcopy(mass_balance)
 
+        if self._mass_balance.ela_h > self.bed.top:
+            warnings.warn(
+                "The ELA is above the top of the glacier. This will prevent the glacier from gaining mass."
+            )
+        elif self._mass_balance.ela_h < self.bed.bottom:
+            warnings.warn(
+                "The ELA is below the bottom of the glacier. It is likely that the glacier will grow out of its domain."
+            )
+
         # Initilaise the flowline
         self._initial_state = self._init_flowline()
         # Set current and model state to None.
@@ -179,6 +188,7 @@ def _to_json(self):
             "Volume [km3]": state.volume_km3,
             "Max ice thickness [m]": state.thick.max(),
             "Max ice velocity [m/yr]": ice_vel,
+            "AAR [%]": self.accumulation_area_ratio * 100,
             "Response time [yrs]": self.response_time,
         }
         return json
@@ -284,7 +294,13 @@ def state_history(self, obj):
 
     @property
     def basal_sliding(self):
-        """Set the sliding parameter of the glacier"""
+        """Set the sliding parameter of the glacier
+
+        Parameters
+        ----------
+        value : float
+            Value setting the basal sliding.
+        """
         return self._basal_sliding
 
     @basal_sliding.setter
@@ -294,23 +310,58 @@ def basal_sliding(self, value):
         Parameters
         ----------
         value : float
-            Value fo the glacier sliding
+            Value setting the basal sliding.
         """
+        # Old default fs of oggm.
+        fs_def = 5.7e-20
+
+        if value > fs_def * 10:
+            warnings.warn("Basal sliding is very high.")
+        # Note that an edu glacier by default has zero basal sliding.
+        elif value < fs_def / 10:
+            warnings.warn("Basal sliding is very low.")
+
         self._basal_sliding = value
 
     @property
     def creep(self):
-        """Set the value for glen_a creep"""
+        """Set the value for the creep parameter A in Glens equation.
+        A value in the range 0 to -50 will be interpreted as a temperature
+        and converted to a creep factor using equation 3.35 in The physics of Glaciers (Cuffey & Paterson).
+        A value above 0 will be interpreted as a creep factor and assigned directly.
+
+        Parameters
+        ----------
+        value : float
+            Temperature in degrees or creep parameter.
+        """
         return self._creep
 
     @creep.setter
     def creep(self, value):
-        """Set the value for glen_a creep
+        """Set the value for the creep parameter A in Glens equation.
+        A value in the range 0 to -50 will be interpreted as a temperature
+        and converted to a creep factor using equation 3.35 in The physics of Glaciers (Cuffey & Paterson).
+        A value above 0 will be interpreted as a creep factor and assigned directly.
 
         Parameters
         ----------
         value : float
+            Temperature in degrees or creep parameter.
         """
+
+        if value > 0:
+            print("Value interpreted as a creep factor.")
+        # Convert temp to A.
+        else:
+            print("Value interpreted as a temperature and converted to a creep factor.")
+            value = cp_glen_a(value)
+        # Use 2.6e-27 as a rec. min for Glen A, from Physics of Glacier: Table 3.4
+        if value < 2.6e-27:
+            warnings.warn("Creep parameter (Glen A) is very small.")
+        elif value > cfg.PARAMS["glen_a"] * 10:
+            warnings.warn("Creep parameter (Glen A) is very large.")
+
         self._creep = value
 
     @history.setter
@@ -370,6 +421,29 @@ def response_time(self):
             response_time = all_vol_diff.time.isel(time=idx) - year_initial
             return response_time.values.item()
 
+    @property
+    def accumulation_area_ratio(self):
+        """The accumulation area ratio. As calculated in the GlacierExplorer app."""
+        if not self.current_state:
+            return np.nan
+        # What is the area of the glacier?
+        total_area = self.current_state.area_m2
+        # Can calculate AAR if there is an accumulation area
+        if np.any(self.current_state.surface_h > self.mass_balance.ela_h):
+            # How large is the accumulation area.
+            accumulation_area = np.sum(
+                np.where(
+                    self.current_state.surface_h > self.mass_balance.ela_h,
+                    self.current_state.bin_area_m2,
+                    0,
+                )
+            )
+
+            return accumulation_area / total_area
+        # if there is no accumulation area we return nan
+        else:
+            return np.nan
+
     def add_temperature_bias(self, bias, duration, noise=None):
         """Add a gradual temperature bias to the future mass balance of the
         glacier. Note that the method is cumulative, i.e. calling it multiple

oggm_edu/glacier_bed.py

@@ -8,6 +8,7 @@
 # Other libraries.
 import numpy as np
 import pandas as pd
+import warnings
 
 # Plotting
 from matplotlib import pyplot as plt
@@ -129,6 +130,10 @@ def __init__(
                 + " Bottom also has to be above 0"
             )
         # Set the resolution.
+        if map_dx <= 50:
+            warnings.warn(
+                "Setting the map resolution below 50 meters might lead to numerical instabilities."
+            )
         self.map_dx = map_dx
         # Do we have a specified slope?
         if slopes:

oggm_edu/tests/test_glacier_bed.py

@@ -117,7 +117,7 @@ def test_non_linear_bed_constructor():
 
     # Make sure we cant provide strange slop angles.
     with pytest.raises(Exception) as e_info:
-        bed = GlacierBed(top=3600, bottom=3000, width=300, slopes=[90])
+        bed = GlacierBed(top=3600, bottom=3000, width=300, slopes=[110])
 
     with pytest.raises(Exception) as e_info:
         _ = GlacierBed(