Add atmosphere.elevation_change.precipitation.temp_lapse_rate

This new parameter make it possible to use "-atmosphere ...,lapse_rate" with a model that
already includes a temperature lapse rate. (For example: now "-atmosphere pik,lapse_rate"
can be used to adjust precipitation in response to changing ice geometry *without* messing
up the near-surface air temperature coming from "-atmosphere pik".)

This takes care of the issue #492 tried to address.

CHANGES.rst

@@ -97,6 +97,13 @@ Changes since v1.2
   contributed by Ronja Reese and Torsten Albrecht.)
 - Implement a mechanism for "optimizing" the till friction angle used by the Mohr-Coulomb
   yield stress model. The implementation is based on the code contributed by T. Albrecht.
+- Add `atmosphere.elevation_change.precipitation.temp_lapse_rate` to the `-atmosphere
+  ...,elevation_change` modifier. Now this parameter is used to compute the temperature
+  change `dT` used to scale precipitation by a factor `exp(C * dT)` with `C =
+  atmosphere.precip_exponential_factor_for_temperature`. We need a separate temperature
+  lapse rate parameter to be able to use this modifier with atmosphere models that include
+  an elevation-dependent near-surface air temperature parameterizations, e.g. `-atmosphere
+  pik,elevation_change`.
 
 Changes from v1.1 to v1.2
 =========================

doc/sphinx/climate_forcing/atmosphere.rst

@@ -275,6 +275,14 @@ The near-surface air temperature is modified using an elevation lapse rate `\gam
 .. math::
    \gamma_T = -\frac{dT}{dz}.
 
+.. warning::
+
+   Some atmosphere models (:ref:`sec-atmosphere-pik`, for example) use elevation-dependent
+   near-surface air temperature parameterizations that include an elevation lapse rate.
+
+   In most cases one should not combine such a temperature parameterization with an
+   additional elevation lapse rate for temperature.
+
 Two methods of adjusting precipitation are available:
 
 - Scaling using an exponential factor
@@ -284,8 +292,8 @@ Two methods of adjusting precipitation are available:
      \mathrm{P} = \mathrm{P_{input}} \cdot \exp(C \cdot \Delta T),
 
   where `C =` :config:`atmosphere.precip_exponential_factor_for_temperature` and `\Delta
-  T` is the temperature difference produced by applying
-  :config:`atmosphere.elevation_change.temperature_lapse_rate`.
+  T` is the temperature difference produced by applying the lapse rate
+  :config:`atmosphere.elevation_change.precipitation.temp_lapse_rate`.
 
   This mechanisms increases precipitation by `100(\exp(C) - 1)` percent for each degree of
   temperature increase.

src/coupler/atmosphere/ElevationChange.cc

@@ -37,6 +37,8 @@ ElevationChange::ElevationChange(IceGrid::ConstPtr grid, std::shared_ptr<Atmosph
   m_precip_lapse_rate = m_config->get_number("atmosphere.elevation_change.precipitation.lapse_rate",
                                              "(kg m-2 / s) / m");
 
+  m_precip_temp_lapse_rate = m_config->get_number("atmosphere.elevation_change.precipitation.temp_lapse_rate",
+                                                  "K / m");
   m_precip_exp_factor = m_config->get_number("atmosphere.precip_exponential_factor_for_temperature");
 
   m_temp_lapse_rate = m_config->get_number("atmosphere.elevation_change.temperature_lapse_rate",
@@ -87,8 +89,10 @@ void ElevationChange::init_impl(const Geometry &geometry) {
                    convert(m_sys, m_precip_lapse_rate, "(kg m-2 / s) / m", "(kg m-2 / year) / km"));
   } else {
     m_log->message(2,
-                   "   precipitation scaling factor with temperature: %3.3f Kelvin-1\n",
-                   m_precip_exp_factor);
+                   "   precipitation scaling factor with temperature: %3.3f Kelvin-1\n"
+                   "   temperature lapse rate: %3.3f K per km\n",
+                   m_precip_exp_factor,
+                   convert(m_sys, m_precip_temp_lapse_rate, "K / m", "K / km"));
   }
 
   ForcingOptions opt(*m_grid->ctx(), "atmosphere.elevation_change");
@@ -107,11 +111,13 @@ void ElevationChange::update_impl(const Geometry &geometry, double t, double dt)
   // temperature and precipitation time series
   m_surface.copy_from(geometry.ice_surface_elevation);
 
+  const auto &reference_surface = *m_reference_surface;
+
   // temperature
   {
     m_temperature->copy_from(m_input_model->mean_annual_temp());
 
-    lapse_rate_correction(m_surface, *m_reference_surface,
+    lapse_rate_correction(m_surface, reference_surface,
                           m_temp_lapse_rate, *m_temperature);
   }
 
@@ -122,16 +128,12 @@ void ElevationChange::update_impl(const Geometry &geometry, double t, double dt)
     switch (m_precip_method) {
     case SCALE:
       {
-        const IceModelVec2S &input_temperature = m_input_model->mean_annual_temp();
-
-        IceModelVec::AccessList list{m_temperature.get(),
-                                     &input_temperature,
-                                     m_precipitation.get()};
+        IceModelVec::AccessList list{&m_surface, &reference_surface, m_precipitation.get()};
 
         for (Points p(*m_grid); p; p.next()) {
           const int i = p.i(), j = p.j();
 
-          double dT = (*m_temperature)(i, j) - input_temperature(i, j);
+          double dT = -m_precip_temp_lapse_rate * (m_surface(i, j) - reference_surface(i, j));
 
           (*m_precipitation)(i, j) *= std::exp(m_precip_exp_factor * dT);
         }
@@ -178,37 +180,37 @@ void ElevationChange::init_timeseries_impl(const std::vector<double> &ts) const
 }
 
 void ElevationChange::temp_time_series_impl(int i, int j, std::vector<double> &result) const {
-  std::vector<double> usurf(m_ts_times.size());
+  std::vector<double> reference_surface(m_ts_times.size());
 
   m_input_model->temp_time_series(i, j, result);
 
-  m_reference_surface->interp(i, j, usurf);
+  m_reference_surface->interp(i, j, reference_surface);
 
   for (unsigned int m = 0; m < m_ts_times.size(); ++m) {
-    result[m] -= m_temp_lapse_rate * (m_surface(i, j) - usurf[m]);
+    result[m] -= m_temp_lapse_rate * (m_surface(i, j) - reference_surface[m]);
   }
 }
 
 void ElevationChange::precip_time_series_impl(int i, int j, std::vector<double> &result) const {
   auto N = m_ts_times.size();
-  std::vector<double> usurf(N);
+  std::vector<double> reference_surface(N);
 
   m_input_model->precip_time_series(i, j, result);
 
-  m_reference_surface->interp(i, j, usurf);
+  m_reference_surface->interp(i, j, reference_surface);
 
   switch (m_precip_method) {
   case SCALE:
     {
       for (unsigned int m = 0; m < N; ++m) {
-        double dT = -m_temp_lapse_rate * (m_surface(i, j) - usurf[m]);
+        double dT = -m_precip_temp_lapse_rate * (m_surface(i, j) - reference_surface[m]);
         result[m] *= std::exp(m_precip_exp_factor * dT);
       }
     }
     break;
   case SHIFT:
     for (unsigned int m = 0; m < N; ++m) {
-      result[m] -= m_precip_lapse_rate * (m_surface(i, j) - usurf[m]);
+      result[m] -= m_precip_lapse_rate * (m_surface(i, j) - reference_surface[m]);
     }
     break;
   }

src/coupler/atmosphere/ElevationChange.hh

@@ -50,8 +50,15 @@ protected:
   enum Method {SCALE, SHIFT};
 
   Method m_precip_method;
+
+  // parameter used when m_precip_method == SHIFT
   double m_precip_lapse_rate;
+
+  // parameters used when m_precip_method == SCALE
+  double m_precip_temp_lapse_rate;
   double m_precip_exp_factor;
+
+  // surface temperature lapse rate
   double m_temp_lapse_rate;
 
   std::shared_ptr<IceModelVec2T> m_reference_surface;

src/pism_config.cdl

@@ -49,7 +49,7 @@ netcdf pism_config {
     pism_config:atmosphere.elevation_change.periodic_type = "flag";
 
     pism_config:atmosphere.elevation_change.precipitation.lapse_rate = 0.0;
-    pism_config:atmosphere.elevation_change.precipitation.lapse_rate_doc = "Elevation lapse rate for the surface mass balance";
+    pism_config:atmosphere.elevation_change.precipitation.lapse_rate_doc = "Elevation lapse rate for the precipitation";
     pism_config:atmosphere.elevation_change.precipitation.lapse_rate_option = "precip_lapse_rate";
     pism_config:atmosphere.elevation_change.precipitation.lapse_rate_type = "number";
     pism_config:atmosphere.elevation_change.precipitation.lapse_rate_units = "(kg m-2 / year) / km";
@@ -60,6 +60,12 @@ netcdf pism_config {
     pism_config:atmosphere.elevation_change.precipitation.method_option = "precip_adjustment";
     pism_config:atmosphere.elevation_change.precipitation.method_type = "keyword";
 
+    pism_config:atmosphere.elevation_change.precipitation.temp_lapse_rate = 0.0;
+    pism_config:atmosphere.elevation_change.precipitation.temp_lapse_rate_doc = "Elevation lapse rate for the surface temperature used to compute `\\Delta T` in the precipitation scaling factor `\\exp(C \\cdot \\Delta T)`";
+    pism_config:atmosphere.elevation_change.precipitation.temp_lapse_rate_option = "precip_temp_lapse_rate";
+    pism_config:atmosphere.elevation_change.precipitation.temp_lapse_rate_type = "number";
+    pism_config:atmosphere.elevation_change.precipitation.temp_lapse_rate_units = "Kelvin / km";
+
     pism_config:atmosphere.elevation_change.temperature_lapse_rate = 0.0;
     pism_config:atmosphere.elevation_change.temperature_lapse_rate_doc = "Elevation lapse rate for the surface temperature";
     pism_config:atmosphere.elevation_change.temperature_lapse_rate_option = "temp_lapse_rate";

test/regression/atmosphere_models.py

@@ -543,6 +543,7 @@ def setUp(self):
         self.dz = 1000.0        # m
         self.dT = -self.dTdz * self.dz / 1000.0
         self.dP = -PISM.util.convert(self.dPdz * self.dz / 1000.0, "kg m-2 year-1", "kg m-2 s-1")
+        self.precip_dTdz = 2.0  # Kelvin per km
 
         self.geometry = PISM.Geometry(self.grid)
 
@@ -553,6 +554,8 @@ def setUp(self):
 
         config.set_number("atmosphere.elevation_change.precipitation.lapse_rate", self.dPdz)
 
+        config.set_number("atmosphere.elevation_change.precipitation.temp_lapse_rate", self.precip_dTdz)
+
         config.set_number("atmosphere.elevation_change.temperature_lapse_rate", self.dTdz)
 
     def tearDown(self):
@@ -589,12 +592,13 @@ def test_atmosphere_elevation_change_scale(self):
         # change surface elevation
         self.geometry.ice_surface_elevation.shift(self.dz)
 
-        # check that the temperature changed accordingly
+        # check that the temperature and precipitation changed accordingly
         modifier.update(self.geometry, 0, 1)
 
         C = config.get_number("atmosphere.precip_exponential_factor_for_temperature")
+        dT = -self.precip_dTdz * self.dz / 1000.0
         P = sample(model.mean_precipitation())
-        dP = np.exp(C * self.dT) * P - P
+        dP = np.exp(C * dT) * P - P
 
         check_modifier(model, modifier, T=self.dT, P=dP,
                        ts=[0.5], Ts=[self.dT], Ps=[dP])