PISM, A Parallel Ice Sheet Model  stable v2.1-1-g6902d5502 committed by Ed Bueler on 2023-12-20 08:38:27 -0800
Go to the documentation of this file.
1 // Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022, 2023 Ed Bueler and Constantine Khroulev
2 //
3 // This file is part of PISM.
4 //
5 // PISM is free software; you can redistribute it and/or modify it under the
6 // terms of the GNU General Public License as published by the Free Software
7 // Foundation; either version 3 of the License, or (at your option) any later
8 // version.
9 //
10 // PISM is distributed in the hope that it will be useful, but WITHOUT ANY
11 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
12 // FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 // details.
14 //
15 // You should have received a copy of the GNU General Public License
16 // along with PISM; if not, write to the Free Software
17 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
22 #include "pism/util/Component.hh"
24 #include "pism/util/Diagnostic.hh"
25 #include <memory>
27 namespace pism {
29 class Context;
30 class Vars;
32 //! @brief Energy balance models and utilities.
33 namespace energy {
35 // Vertical grid information for BTU_Full.
36 struct BTUGrid {
37  BTUGrid(std::shared_ptr<const Context> ctx);
38  static BTUGrid FromOptions(std::shared_ptr<const Context> ctx);
40  unsigned int Mbz; // number of vertical levels
41  double Lbz; // depth of the bed thermal layer
42 };
44 //! Given the temperature of the top of the bedrock, for the duration of one time-step, provides upward geothermal flux at that interface at the end of the time-step.
45 /*!
46  The geothermal flux actually applied to the base of an ice sheet is dependent, over time,
47  on the temperature of the basal ice itself. The purpose of a bedrock thermal layer
48  in an ice sheet model is to implement this dependency by using a physical model
49  for the temperature within that layer, the upper lithosphere. Because the
50  upper part of the lithosphere stores or releases energy into the ice,
51  the typical lithosphere geothermal flux rate is not the same thing as the
52  geothermal flux applied to the base of the ice. This issue has long been
53  recognized by ice sheet modelers [%e.g. \ref RitzFabreLetreguilly].
55  For instance, suppose the ice sheet is in a balanced state in which the geothermal
56  flux deep in the crust is equal to the heat flux into the ice base. If the
57  near-surface ice cools from this state then, because the ice temperature gradient
58  is now greater in magnitude, between the warm bedrock and the cooler ice, the ice
59  will for some period receive more than the deep geothermal flux rate. Similarly,
60  if the ice warms from the balanced state then the temperature difference with
61  the bedrock has become smaller and the magnitude of the ice basal heat flux will
62  be less than the deep geothermal rate.
64  We regard the lithosphere geothermal flux rate, which is applied in this model
65  to the base of the bedrock thermal layer, as a time-independent quantity. This
66  concept is the same as in all published ice sheet models, to our knowledge.
68  Because the relevant layer of bedrock below an ice sheet is typically shallow,
69  modeling the bedrock temperature is quite simple.
70  Let \f$T_b(t,x,y,z)\f$ be the temperature of the bedrock layer, for elevations
71  \f$-L_b \le z \le 0\f$. In this routine, \f$z=0\f$ refers to the top of the
72  bedrock, the ice/bedrock interface. (Note \f$z=0\f$ is the base of the ice in
73  IceModel, and thus a different location if ice is floating.)
74  Let \f$G\f$ be the lithosphere geothermal flux rate, namely the PISM input
75  variable `bheatflx`; see Related Page \ref std_names . Let \f$k_b\f$
76  = `bedrock_thermal_conductivity` in pism_config.cdl) be the constant thermal
77  conductivity of the upper lithosphere. In these terms the actual
78  upward heat flux into the ice/bedrock interface is the quantity,
79  \f[G_0 = -k_b \frac{\partial T_b}{\partial z}.\f]
80  This is the \e output of the method top_heat_flux() in this class.
82  The evolution equation solved in this class, for which a timestep is done by the
83  update() method, is the standard 1D heat equation
84  \f[\rho_b c_b \frac{\partial T_b}{\partial t} = k_b \frac{\partial^2 T_b}{\partial z^2}\f]
85  where \f$\rho_b\f$ = `bedrock_thermal_density` and \f$c_b\f$ =
86  `bedrock_thermal_specific_heat_capacity` in pism_config.cdl.
88  If `n_levels` >= 3 then everything is the general case. The lithospheric temperature
89  in `temp` is saved in files as `litho_temp`. The top_heat_flux()
90  method uses second-order differencing to compute the values of \f$G_0\f$.
92  If `n_levels` <= 1 then this object becomes very simplified: there is no internal
93  state in array::Array3D temp. The update() and allocate() methods are null,
94  and the top_heat_flux() method does nothing other than to copy the
95  field \f$G\f$ = `bheatflx` into `result`.
97  If `n_levels` == 2 then everything is the general case except that
98  top_heat_flux() method uses first-order differencing to compute the
99  values of \f$G_0\f$.
100 */
101 class BedThermalUnit : public Component {
102 public:
104  static std::shared_ptr<BedThermalUnit> FromOptions(std::shared_ptr<const Grid> g,
105  std::shared_ptr<const Context> ctx);
107  BedThermalUnit(std::shared_ptr<const Grid> g);
109  virtual ~BedThermalUnit() = default;
111  typedef std::shared_ptr<BedThermalUnit> Ptr;
112  typedef std::shared_ptr<const BedThermalUnit> ConstPtr;
114  void init(const InputOptions &opts);
116  //! Return the upward heat flux through the top surface of the bedrock thermal layer.
119  //! Return the upward heat flux through the bottom surface of the bedrock thermal layer.
122  void update(const array::Scalar &bedrock_top_temperature,
123  double t, double dt);
125  double vertical_spacing() const;
126  double depth() const;
128  unsigned int Mz() const;
130 protected:
131  virtual void initialize_bottom_surface_flux();
133  virtual void init_impl(const InputOptions &opts);
135  virtual void update_impl(const array::Scalar &bedrock_top_temperature,
136  double t, double dt) = 0;
138  virtual double vertical_spacing_impl() const = 0;
139  virtual double depth_impl() const = 0;
140  virtual unsigned int Mz_impl() const = 0;
142  virtual void define_model_state_impl(const File &output) const;
143  virtual void write_model_state_impl(const File &output) const;
145  virtual DiagnosticList diagnostics_impl() const;
147  //! upward heat flux through the bottom surface of the bed thermal layer
150  //! upward heat flux through the top surface of the bed thermal layer
152 };
154 class BTU_geothermal_flux_at_ground_level : public Diag<BedThermalUnit> {
155 public:
157 protected:
158  virtual std::shared_ptr<array::Array> compute_impl() const;
159 };
161 } // end of namespace energy
162 } // end of namespace pism
164 #endif /* _PISMBEDTHERMALUNIT_H_ */
A class defining a common interface for most PISM sub-models.
Definition: Component.hh:118
A template derived from Diagnostic, adding a "Model".
Definition: Diagnostic.hh:166
High-level PISM I/O class.
Definition: File.hh:56
BTU_geothermal_flux_at_ground_level(const BedThermalUnit *m)
virtual std::shared_ptr< array::Array > compute_impl() const
BedThermalUnit(std::shared_ptr< const Grid > g)
virtual void initialize_bottom_surface_flux()
std::shared_ptr< const BedThermalUnit > ConstPtr
void update(const array::Scalar &bedrock_top_temperature, double t, double dt)
const array::Scalar & flux_through_top_surface() const
Return the upward heat flux through the top surface of the bedrock thermal layer.
virtual ~BedThermalUnit()=default
array::Scalar m_bottom_surface_flux
upward heat flux through the bottom surface of the bed thermal layer
virtual void init_impl(const InputOptions &opts)
Initialize the bedrock thermal unit.
void init(const InputOptions &opts)
virtual double depth_impl() const =0
static std::shared_ptr< BedThermalUnit > FromOptions(std::shared_ptr< const Grid > g, std::shared_ptr< const Context > ctx)
virtual void define_model_state_impl(const File &output) const
The default (empty implementation).
virtual void write_model_state_impl(const File &output) const
The default (empty implementation).
virtual double vertical_spacing_impl() const =0
virtual DiagnosticList diagnostics_impl() const
virtual unsigned int Mz_impl() const =0
const array::Scalar & flux_through_bottom_surface() const
Return the upward heat flux through the bottom surface of the bedrock thermal layer.
array::Scalar m_top_surface_flux
upward heat flux through the top surface of the bed thermal layer
std::shared_ptr< BedThermalUnit > Ptr
virtual void update_impl(const array::Scalar &bedrock_top_temperature, double t, double dt)=0
Given the temperature of the top of the bedrock, for the duration of one time-step,...
static const double g
Definition: exactTestP.cc:36
std::map< std::string, Diagnostic::Ptr > DiagnosticList
Definition: Diagnostic.hh:125
static BTUGrid FromOptions(std::shared_ptr< const Context > ctx)
BTUGrid(std::shared_ptr< const Context > ctx)