◆ bootstrap_ice_temperature()

void pism::energy::bootstrap_ice_temperature	(	const array::Scalar &	ice_thickness,
		const array::Scalar &	ice_surface_temp,
		const array::Scalar &	surface_mass_balance,
		const array::Scalar &	basal_heat_flux,
		array::Array3D &	result
	)

Create a temperature field within the ice from provided ice thickness, surface temperature, surface mass balance, and geothermal flux.

In bootstrapping we need to determine initial values for the temperature within the ice (and the bedrock). There are various data available at bootstrapping, but not the 3D temperature field needed as initial values for the temperature. Here we take a "guess" based on an assumption of steady state and a simple model of the vertical velocity in the column. The rule is certainly heuristic but it seems to work well anyway.

The result is not the temperature field which is in steady state with the ice dynamics. Spinup is most-definitely needed in many applications. Such spinup usually starts from the temperature field computed by this procedure and then runs for a long time (e.g. \(10^4\) to \(10^6\) years), possibly with fixed geometry, to get closer to thermomechanically-coupled equilibrium.

Consider a horizontal grid point. Suppose the surface temperature \(T_s\), surface mass balance \(m\), and geothermal flux \(g\) are given at that location. Within the column denote the temperature by \(T(z)\) at height \(z\) above the base of the ice. Suppose the column of ice has height \(H\), the ice thickness.

There are two alternative bootstrap methods determined by the configuration parameter config.get_number("bootstrapping.temperature_heuristic")). Allowed values are "smb" and "quartic_guess".

If the smb method is chosen, which is the default, and if \(m>0\), then the method sets the ice temperature to the solution of the steady problem [Paterson]
\[\rho_i c w \frac{\partial T}{\partial z} = k_i \frac{\partial^2 T}{\partial z^2} \qquad \text{with boundary conditions} \qquad T(H) = T_s \quad \text{and} \quad \frac{\partial T}{\partial z}(0) = - \frac{g}{k_i}, \]
where the vertical velocity is linear between the surface value \(w=-m\) and a velocity of zero at the base:
\[w(z) = - m z / H.\]
(Note that because \(m>0\), this vertical velocity is downward.) This is a two-point boundary value problem for a linear ODE. In fact, if \(K = k_i / (\rho_i c)\) then we can write the ODE as
\[K T'' + \frac{m z}{H} T' = 0.\]
Then let
\[C_0 = \frac{g \sqrt{\pi H K}}{k_i \sqrt{2 m}}, \qquad \gamma_0 = \sqrt{\frac{mH}{2K}}.\]
(Note \(\gamma_0\) is, up to a constant, the square root of the Peclet number [Paterson]; compare [vanderWeletal2013].) The solution to the two-point boundary value problem is then
\[T(z) = T_s + C_0 \left(\operatorname{erf}(\gamma_0) - \operatorname{erf}\left(\gamma_0 \frac{z}{H}\right)\right).\]
If usesmb is true and \(m \le 0\), then the velocity in the column, relative to the base, is taken to be zero. Thus the solution is
\[ T(z) = \frac{g}{k_i} \left( H - z \right) + T_s, \]
a straight line whose slope is determined by the geothermal flux and whose value at the ice surface is the surface temperature, \(T(H) = T_s\).
If the quartic_guess method is chosen, the "quartic guess" formula which was in older versions of PISM is used. Namely, within the ice we set
\[T(z) = T_s + \alpha (H-z)^2 + \beta (H-z)^4\]
where \(\alpha,\beta\) are chosen so that
\[\frac{\partial T}{\partial z}\Big|_{z=0} = - \frac{g}{k_i} \qquad \text{and} \qquad \frac{\partial T}{\partial z}\Big|_{z=H/4} = - \frac{g}{2 k_i}.\]
The purpose of the second condition is that when ice is advecting downward then the temperature gradient is much larger in roughly the bottom quarter of the ice column. However, without the surface mass balance, much less the solution of the stress balance equations, we cannot estimate the vertical velocity, so we make such a rough guess.

In either case the temperature within the ice is not allowed to exceed the pressure-melting temperature.

We set \(T(z)=T_s\) above the top of the ice.

This method determines \(T(0)\), the ice temperature at the ice base. This temperature is used by BedThermalUnit::bootstrap() to determine a bootstrap temperature profile in the bedrock.

Definition at line 327 of file utilities.cc.

References pism::ParallelSection::check(), pism::ParallelSection::failed(), pism::RuntimeError::formatted(), G, pism::array::Array3D::get_column(), pism::VariableMetadata::get_number(), pism::array::Array::grid(), ice_temperature_guess(), ice_temperature_guess_smb(), pism::hydrology::K(), pism::k, pism::array::Array::metadata(), pism::array::min(), PISM_ERROR_LOCATION, pism::diagnostics::SMB, and pism::array::Array::update_ghosts().

Referenced by bootstrap_ice_enthalpy(), pism::energy::TemperatureModel::bootstrap_impl(), and pism::energy::TemperatureModel::initialize_impl().