tBTU_Full.hh - pism - [fork] customized build of PISM, the parallel ice sheet model (tillflux branch)
HTML git clone git://src.adamsgaard.dk/pism
DIR Log
DIR Files
DIR Refs
DIR LICENSE
---
tBTU_Full.hh (5894B)
---
1 /* Copyright (C) 2016, 2017, 2019 PISM Authors
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
20 #ifndef BTU_FULL_H
21 #define BTU_FULL_H
22
23 #include "BedThermalUnit.hh"
24 #include "pism/util/Context.hh"
25
26 namespace pism {
27 namespace energy {
28
29 class BedrockColumn;
30
31 //! @brief Given the temperature of the top of the bedrock, for the duration of one time-step,
32 //! provides upward geothermal flux at that interface at the end of the time-step.
33 /*!
34 The geothermal flux actually applied to the base of an ice sheet is dependent, over time,
35 on the temperature of the basal ice itself. The purpose of a bedrock thermal layer
36 in an ice sheet model is to implement this dependency by using a physical model
37 for the temperature within that layer, the upper lithosphere. Because the
38 upper part of the lithosphere stores or releases energy into the ice,
39 the typical lithosphere geothermal flux rate is not the same thing as the
40 geothermal flux applied to the base of the ice. This issue has long been
41 recognized by ice sheet modelers [%e.g. \ref RitzFabreLetreguilly].
42
43 For instance, suppose the ice sheet is in a balanced state in which the geothermal
44 flux deep in the crust is equal to the heat flux into the ice base. If the
45 near-surface ice cools from this state then, because the ice temperature gradient
46 is now greater in magnitude, between the warm bedrock and the cooler ice, the ice
47 will for some period receive more than the deep geothermal flux rate. Similarly,
48 if the ice warms from the balanced state then the temperature difference with
49 the bedrock has become smaller and the magnitude of the ice basal heat flux will
50 be less than the deep geothermal rate.
51
52 We regard the lithosphere geothermal flux rate, which is applied in this model
53 to the base of the bedrock thermal layer, as a time-independent quantity. This
54 concept is the same as in all published ice sheet models, to our knowledge.
55
56 Because the relevant layer of bedrock below an ice sheet is typically shallow,
57 modeling the bedrock temperature is quite simple.
58 Let \f$T_b(t,x,y,z)\f$ be the temperature of the bedrock layer, for elevations
59 \f$-L_b \le z \le 0\f$. In this routine, \f$z=0\f$ refers to the top of the
60 bedrock, the ice/bedrock interface. (Note \f$z=0\f$ is the base of the ice in
61 IceModel, and thus a different location if ice is floating.)
62 Let \f$G\f$ be the lithosphere geothermal flux rate, namely the PISM input
63 variable `bheatflx`; see Related Page \ref std_names . Let \f$k_b\f$
64 = `bedrock_thermal_conductivity` in pism_config.cdl) be the constant thermal
65 conductivity of the upper lithosphere. In these terms the actual
66 upward heat flux into the ice/bedrock interface is the quantity,
67 \f[G_0 = -k_b \frac{\partial T_b}{\partial z}.\f]
68 This is the \e output of the method flux_through_top_surface() in this class.
69
70 The evolution equation solved in this class, for which a timestep is done by the
71 update() method, is the standard 1D heat equation
72 \f[\rho_b c_b \frac{\partial T_b}{\partial t} = k_b \frac{\partial^2 T_b}{\partial z^2}\f]
73 where \f$\rho_b\f$ = `bedrock_thermal_density` and \f$c_b\f$ =
74 `bedrock_thermal_specific_heat_capacity` in pism_config.cdl.
75
76 If 3 or more levels are used then everything is the general case. The lithospheric temperature in
77 `temp` is saved in files as `litho_temp`. The flux_through_top_surface() method uses second-order
78 differencing to compute the values of \f$G_0\f$.
79
80 If 2 levels are used then everything is the general case except that flux_through_top_surface()
81 method uses first-order differencing to compute the values of \f$G_0\f$.
82 */
83 class BTU_Full : public BedThermalUnit {
84 public:
85 BTU_Full(IceGrid::ConstPtr g, const BTUGrid &vertical_grid);
86 virtual ~BTU_Full();
87
88 //! Bedrock thermal layer temperature field.
89 const IceModelVec3Custom& temperature() const;
90
91 protected:
92 virtual void bootstrap(const IceModelVec2S &bedrock_top_temperature);
93
94 virtual void init_impl(const InputOptions &opts);
95
96 virtual double vertical_spacing_impl() const;
97 virtual double depth_impl() const;
98 virtual unsigned int Mz_impl() const;
99
100 virtual MaxTimestep max_timestep_impl(double my_t) const;
101
102 using BedThermalUnit::update_impl;
103 virtual void update_impl(const IceModelVec2S &bedrock_top_temperature,
104 double t, double dt);
105
106 virtual void define_model_state_impl(const File &output) const;
107 virtual void write_model_state_impl(const File &output) const;
108 protected:
109 //! bedrock thermal layer temperature, in degrees Kelvin; part of state; uses equally-spaced
110 //! layers.
111 IceModelVec3Custom::Ptr m_temp;
112
113 //! bedrock thermal conductivity
114 double m_k;
115 //! diffusivity of the heat flow within the bedrock layer
116 double m_D;
117
118 //! number of vertical levels within the bedrock
119 unsigned int m_Mbz;
120 //! thickness of the bedrock layer, in meters
121 double m_Lbz;
122
123 //! true if the model needs to "bootstrap" the temperature field during the first time step
124 bool m_bootstrapping_needed;
125
126 void update_flux_through_top_surface();
127
128 std::shared_ptr<BedrockColumn> m_column;
129 };
130
131 } // end of namespace energy
132 } // end of namespace pism
133
134
135 #endif /* BTU_FULL_H */