tBedThermalUnit.hh - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
(HTM) git clone git://src.adamsgaard.dk/pism
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
tBedThermalUnit.hh (6978B)
---
1 // Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019 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
18
19 #ifndef _PISMBEDTHERMALUNIT_H_
20 #define _PISMBEDTHERMALUNIT_H_
21
22 #include "pism/util/Component.hh"
23 #include "pism/util/iceModelVec3Custom.hh"
24
25 #include "pism/util/Diagnostic.hh"
26
27 namespace pism {
28
29 class Vars;
30
31 //! @brief Energy balance models and utilities.
32 namespace energy {
33
34 // Vertical grid information for BTU_Full.
35 struct BTUGrid {
36 BTUGrid(Context::ConstPtr ctx);
37 static BTUGrid FromOptions(Context::ConstPtr ctx);
38
39 unsigned int Mbz; // number of vertical levels
40 double Lbz; // depth of the bed thermal layer
41 };
42
43 //! 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.
44 /*!
45 The geothermal flux actually applied to the base of an ice sheet is dependent, over time,
46 on the temperature of the basal ice itself. The purpose of a bedrock thermal layer
47 in an ice sheet model is to implement this dependency by using a physical model
48 for the temperature within that layer, the upper lithosphere. Because the
49 upper part of the lithosphere stores or releases energy into the ice,
50 the typical lithosphere geothermal flux rate is not the same thing as the
51 geothermal flux applied to the base of the ice. This issue has long been
52 recognized by ice sheet modelers [%e.g. \ref RitzFabreLetreguilly].
53
54 For instance, suppose the ice sheet is in a balanced state in which the geothermal
55 flux deep in the crust is equal to the heat flux into the ice base. If the
56 near-surface ice cools from this state then, because the ice temperature gradient
57 is now greater in magnitude, between the warm bedrock and the cooler ice, the ice
58 will for some period receive more than the deep geothermal flux rate. Similarly,
59 if the ice warms from the balanced state then the temperature difference with
60 the bedrock has become smaller and the magnitude of the ice basal heat flux will
61 be less than the deep geothermal rate.
62
63 We regard the lithosphere geothermal flux rate, which is applied in this model
64 to the base of the bedrock thermal layer, as a time-independent quantity. This
65 concept is the same as in all published ice sheet models, to our knowledge.
66
67 Because the relevant layer of bedrock below an ice sheet is typically shallow,
68 modeling the bedrock temperature is quite simple.
69 Let \f$T_b(t,x,y,z)\f$ be the temperature of the bedrock layer, for elevations
70 \f$-L_b \le z \le 0\f$. In this routine, \f$z=0\f$ refers to the top of the
71 bedrock, the ice/bedrock interface. (Note \f$z=0\f$ is the base of the ice in
72 IceModel, and thus a different location if ice is floating.)
73 Let \f$G\f$ be the lithosphere geothermal flux rate, namely the PISM input
74 variable `bheatflx`; see Related Page \ref std_names . Let \f$k_b\f$
75 = `bedrock_thermal_conductivity` in pism_config.cdl) be the constant thermal
76 conductivity of the upper lithosphere. In these terms the actual
77 upward heat flux into the ice/bedrock interface is the quantity,
78 \f[G_0 = -k_b \frac{\partial T_b}{\partial z}.\f]
79 This is the \e output of the method top_heat_flux() in this class.
80
81 The evolution equation solved in this class, for which a timestep is done by the
82 update() method, is the standard 1D heat equation
83 \f[\rho_b c_b \frac{\partial T_b}{\partial t} = k_b \frac{\partial^2 T_b}{\partial z^2}\f]
84 where \f$\rho_b\f$ = `bedrock_thermal_density` and \f$c_b\f$ =
85 `bedrock_thermal_specific_heat_capacity` in pism_config.cdl.
86
87 If `n_levels` >= 3 then everything is the general case. The lithospheric temperature
88 in `temp` is saved in files as `litho_temp`. The top_heat_flux()
89 method uses second-order differencing to compute the values of \f$G_0\f$.
90
91 If `n_levels` <= 1 then this object becomes very simplified: there is no internal
92 state in IceModelVec3 temp. The update() and allocate() methods are null,
93 and the top_heat_flux() method does nothing other than to copy the
94 field \f$G\f$ = `bheatflx` into `result`.
95
96 If `n_levels` == 2 then everything is the general case except that
97 top_heat_flux() method uses first-order differencing to compute the
98 values of \f$G_0\f$.
99 */
100 class BedThermalUnit : public Component {
101 public:
102
103 static BedThermalUnit* FromOptions(IceGrid::ConstPtr g,
104 Context::ConstPtr ctx);
105
106 BedThermalUnit(IceGrid::ConstPtr g);
107
108 virtual ~BedThermalUnit();
109
110 typedef std::shared_ptr<BedThermalUnit> Ptr;
111 typedef std::shared_ptr<const BedThermalUnit> ConstPtr;
112
113 void init(const InputOptions &opts);
114
115 //! Return the upward heat flux through the top surface of the bedrock thermal layer.
116 const IceModelVec2S& flux_through_top_surface() const;
117
118 //! Return the upward heat flux through the bottom surface of the bedrock thermal layer.
119 const IceModelVec2S& flux_through_bottom_surface() const;
120
121 void update(const IceModelVec2S &bedrock_top_temperature,
122 double t, double dt);
123
124 double vertical_spacing() const;
125 double depth() const;
126
127 unsigned int Mz() const;
128
129 protected:
130 virtual void initialize_bottom_surface_flux();
131
132 virtual void init_impl(const InputOptions &opts);
133
134 virtual void update_impl(const IceModelVec2S &bedrock_top_temperature,
135 double t, double dt) = 0;
136
137 virtual double vertical_spacing_impl() const = 0;
138 virtual double depth_impl() const = 0;
139 virtual unsigned int Mz_impl() const = 0;
140
141 virtual void define_model_state_impl(const File &output) const;
142 virtual void write_model_state_impl(const File &output) const;
143
144 virtual DiagnosticList diagnostics_impl() const;
145 protected:
146 //! upward heat flux through the bottom surface of the bed thermal layer
147 IceModelVec2S m_bottom_surface_flux;
148
149 //! upward heat flux through the top surface of the bed thermal layer
150 IceModelVec2S m_top_surface_flux;
151 };
152
153 class BTU_geothermal_flux_at_ground_level : public Diag<BedThermalUnit> {
154 public:
155 BTU_geothermal_flux_at_ground_level(const BedThermalUnit *m);
156 protected:
157 virtual IceModelVec::Ptr compute_impl() const;
158 };
159
160 } // end of namespace energy
161 } // end of namespace pism
162
163 #endif /* _PISMBEDTHERMALUNIT_H_ */
164