PISM fails with negative geothermal flux about pism

[mariazeitz]

https://cloud.pik-potsdam.de/index.php/s/3fqGyYRi2bwgaSX

[talbrecht]

I had this issue a few times before and used ncap2 to set negative values to 0 before running PISM (as an artefact of the regridding). But I wonder whether negative values could make physical sense, or could we just set a valid_min?

[ckhroulev]

In summary: don't regrid litho_temp or regrid both litho_temp and enthalpy.

Here's what's going on.

When you start PISM with -bootstrap it reads in 2D fields from an input file (in this case: longitude, latitude, ice thickness, bedrock elevation, bedrock uplift rate, till water amount, geothermal flux, basal melt rate) and fills 3D fields (temperature in the bedrock thermal layer and ice enthalpy) using heuristics.

When the bedrock thermal layer is active, its temperature (litho_temp) is set to the steady-state solution of the heat equation in each individual vertical column with a Neumann BC corresponding to the geothermal flux at the bottom boundary and a Dirichlet ensuring continuity of temperature at the bed-ice interface at the top boundary. (This is a straight line with the slope determined by the geothermal flux.)

The ice temperature in each ice column is set to the solution of a steady one-dimensional differential equation in which conduction and vertical advection are in balance, and the vertical velocity linearly-interpolates between the surface mass balance rate at the top and zero at the bottom. This uses the heat flux through the bottom ice surface (computed using litho_temp or geothermal flux if litho_temp is not available), ice thickness, surface mass balance, and top surface temperature. Once this is done, ice temperature is converted to ice enthalpy assuming that all ice is cold (i.e. has a zero water fraction).

If you add -regrid_file ... -regrid_vars litho_temp PISM will skip the litho_temp heuristic. The same is true for enthalpy.

There is one tricky case if we need to fill both the temperature of the bedrock thermal layer and the ice temperature:

• the bed thermal layer heuristic needs the temperature at the bottom of the ice column
• the ice column heuristic needs the flux through the top surface of the bed thermal layer column.
  To work around that we use the geothermal flux as the heat flux through the bottom of an ice column and then fill bed thermal layer columns using resulting ice temperatures for the top surface BC.

The error message PISM ERROR: geothermal flux G(23,77) = -0.037027 < 0.0 W m-2 is admittedly confusing, because it is really about the heat flux through the bottom surface of the ice, not geothermal flux.

I need to fix that.

Here's the underlying issue: the heuristic (as implemented) cannot handle the case in which heat is flowing from the ice into the bedrock as in the input file from this issue (see input.nc):

This is a limitation, of course, but more thinking is needed to decide if it is worth removing. One could add a some code so that if

• the bed below the ice is at the pressure-melting temperature of the ice
• the top surface temperature is at or above melting
• the SMB is negative

then the whole column is set to the pressure melting temperature (as a function of depth).

However, this configuration (together with the bed thermal layer column) does not seem to be close to a steady state to me.

[ckhroulev]

The error message is fixed in the dev branch, so the fix will be in the next release.

Now instead of complaining about geothermal flux it prints something like this:

PISM ERROR: Negative upward heat flux (-0.000144 W m-2) through the bottom of the ice column
            is not allowed by PISM's ice temperature bootstrapping method.
            Please check the temperature of the bedrock thermal layer at i=25, j=77.