Surface mass and energy process model components ¶

Contents

Surface mass and energy process model components

The “invisible” model ¶

Options: -surface simple
Variables: none
C++ class: pism::surface::Simple

This is the simplest “surface model” available in PISM, enabled using -surface simple. Its job is to re-interpret precipitation as climatic mass balance, and to re-interpret mean annual near-surface (2m) air temperature as the temperature of the ice at the depth at which firn processes cease to change the temperature of the ice. (I.e. the temperature below the firn.) This implies that there is no melt. Though primitive, this model component may be desired in cold environments (e.g. East Antarctic ice sheet) in which melt is negligible and heat from firn processes is ignored.

Reading top-surface boundary conditions from a file ¶

Options: -surface given
Variables: ice_surface_temp, climatic_mass_balance \(kg / (m^{2} s)\)
C++ class: pism::surface::Given

Note

This is the default choice.

This model component was created to force PISM with sampled (possibly periodic) climate data by reading ice upper surface boundary conditions from a file. These fields are provided directly to the ice dynamics code (see Climate inputs, and their interface with ice dynamics for details).

PISM will stop if variables ice_surface_temp (ice temperature at the ice surface but below firn) and climatic_mass_balance (top surface mass flux into the ice) are not present in the input file.

A file foo.nc used with -surface given -surface_given_file foo.nc may contain several records. If this file contains one record (i.e. fields corresponding to one time value only), provided forcing data is interpreted as time-independent. Variables time and time_bnds should specify model times corresponding to individual records.

For example, to use monthly periodic forcing with a period of 1 year starting at the beginning of \(1980\) (let’s use the \(360\)-day calendar for simplicity), create a file (say, “foo.nc”) with 12 records. The time variable may contain \(15, 45, 75, \dots, 345\) (mid-month for all \(12\) months) and have the units of “days since 1980-1-1”. (It is best to avoid units of “months” and “years” because their meanings depend on the calendar.) Next, add the time_bounds variable for the time dimension with the values \(0, 30, 30, 60, \dots\) specifying times corresponding to beginnings and ends of records for each month and set the time:bounds attribute accordingly. Now run

pismr -surface given -surface_given_file foo.nc -surface.given.periodic

See Using time-dependent forcing for more information.

Note

This surface model ignores the atmosphere model selection made using the option -atmosphere.
PISM can handle files with virtually any number of records: it will read and store in memory at most input.forcing.buffer_size records at any given time (default: 60, or 5 years’ worth of monthly fields).
when preparing a file for use with this model, it is best to use the t,y,x variable storage order: files using this order can be read in faster than ones using the t,x,y order, for reasons explained in the User’s Manual.

To change the storage order in a NetCDF file, use ncpdq:
```
ncpdq -a t,y,x input.nc output.nc
```
will copy data from input.nc into output.nc, changing the storage order to t,y,x at the same time.

Parameters

Prefix: surface.given.

file Name of the file containing climate forcing fields.
periodic (no) If true, interpret forcing data as periodic in time
smb_max (91000 kg m-2 year-1) Maximum climatic mass balance value (used to check input data). Corresponds to 100 m/year ice equivalent.

Elevation-dependent temperature and mass balance ¶

Options: -surface elevation
Variables: none
C++ class: pism::surface::Elevation

This surface model component parameterizes the ice surface temperature \(T_{h}\) = ice_surface_temp and the mass balance \(m\) = climatic_mass_balance as piecewise-linear functions of surface elevation \(h\).

The option -ice_surface_temp (list of 4 numbers) determines the surface temperature using the 4 parameters \(\T{min}\), \(\T{max}\), \(\h{min}\), \(\h{max}\). Let

\[\diff{T}{h} = (\T{max} - \T{min}) / (\h{max} - \h{min})\]

be the temperature gradient. Then

\[\begin{split}T(x,y) = \begin{cases} \T{min}, & h(x,y) \le \h{min}, \\ \T{min} + \diff{T}{h} \, (h(x,y) - \h{min}), & \h{min} < h(x,y) < \h{max}, \\ \T{max}, & \h{max} \le h(x,y). \end{cases}\end{split}\]

The option -climatic_mass_balance (list of 5 numbers) determines the surface mass balance using the 5 parameters \(\m{min}\), \(\m{max}\), \(\h{min}\), \(\h{ELA}\), \(\h{max}\). Let

\[\diff{\m{abl}}{h} = -\m{min} / (\ELA - \h{min})\]

and

\[\diff{\m{acl}}{h} = \m{max} / (\h{max} - \ELA)\]

be the mass balance gradient in the ablation and in the accumulation area, respectively. Then

\[\begin{split}m(x,y) = \begin{cases} \m{min}, & h(x,y) \le \h{min}, \\ \diff{\m{abl}}{h} \, (h(x,y) - \ELA), & \h{min} < h(x,y) \le \ELA, \\ \diff{\m{acl}}{h} \, (h(x,y) - \ELA), & \ELA < h(x,y) \le \h{max}, \\ \m{max}, & \h{max} < h(x,y). \end{cases}\end{split}\]

The option -climatic_mass_balance_limits (list of 2 numbers) limits the mass balance below \(\h{min}\) to \(\ms{min}\) and above \(\h{max}\) to \(\ms{max}\), thus

\[\begin{split}m(x,y) = \begin{cases} m^{*}_{\text{min}}, & h(x,y) \le \h{min}, \\ \diff{\m{abl}}{h} \, (h(x,y) - \ELA), & \h{min} < h(x,y) \le \ELA, \\ \diff{\m{acl}}{h} \, (h(x,y) - \ELA), & \ELA < h(x,y) \le \h{max}, \\ m^{*}_{\text{max}}, & \h{max} < h(x,y). \end{cases}\end{split}\]

Note: this surface model ignores the atmosphere model selection made using the -atmosphere option.

Temperature-index scheme ¶

Options: -surface pdd
Variables: air_temp_sd, snow_depth
C++ class: pism::surface::TemperatureIndex

The default PDD model used by PISM, turned on by option -surface pdd, is based on [153] and EISMINT-Greenland intercomparison (see [148]).

Our model computes the solid (snow) precipitation rate using the air temperature threshold with a linear transition. All precipitation during periods with air temperatures above surface.pdd.air_temp_all_precip_as_rain (default of \(2^\circ C\)) is interpreted as rain; all precipitation during periods with air temperatures below surface.pdd.air_temp_all_precip_as_snow (default of \(0^\circ C\)) is interpreted as snow.

For long-term simulations, a PDD model generally uses an idealized seasonal temperature cycle. “White noise” is added to this cycle to simulate additional daily variability associated to the vagaries of weather. This additional random variation is quite significant, as the seasonal cycle may never reach the melting point but that point may be reached with some probability, in the presence of the daily variability, and thus melt may occur. Concretely, a normally-distributed, mean zero random temperature increment is added to the seasonal cycle. There is no assumed spatial correlation of daily variability. The standard deviation of the daily variability is controlled by configuration parameters with the prefix surface.pdd.std_dev.:

file The name of the file to read air_temp_sd (standard deviation of air temperature) from.
lapse_lat_base (72 degree_north) standard deviationis is a function of latitude, with value surface.pdd.std_dev.value at this latitude; this value is only active if surface.pdd.std_dev.lapse_lat_rate is nonzero
lapse_lat_rate (0 Kelvin / degree_north) standard deviation is a function of latitude, with rate of change with respect to latitude given by this constant
param_a (-0.15) Parameter \(a\) in \(\Sigma = aT + b\), with \(T\) in degrees C. Used only if surface.pdd.std_dev.use_param is set to yes.
param_b (0.66 Kelvin) Parameter \(b\) in \(\Sigma = aT + b\), with \(T\) in degrees C. Used only if surface.pdd.std_dev.use_param is set to yes.
periodic (no) If true, interpret air_temp_sd read from surface.pdd.std_dev.file as periodic in time
use_param (no) Parameterize standard deviation as a linear function of air temperature over ice-covered grid cells. The region of application is controlled by geometry.ice_free_thickness_standard.
value (5 Kelvin) standard deviation of daily temp variation; = EISMINT-Greenland value [148]

A file foo.nc used with -surface pdd -pdd_sd_file foo.nc should contain standard deviation of near-surface air temperature in variable air_temp_sd, and the corresponding time coordinate in variable time. If -pdd_sd_file is not set, PISM uses a constant value for standard deviation, which is set by the configuration parameter surface.pdd.std_dev.value. The default value is \(5.0\) degrees [148]. However, this approach is not recommended as it induces significant errors in modeled surface mass balance in both ice-covered and ice-free regions [154], [155].

Over ice-covered grid cells, daily variability can also be parameterized as a linear function of near-surface air temperature \(\sigma = a \cdot T + b\) using the surface.pdd.std_dev.use_param configuration flag, and the corresponding parameters surface.pdd.std_dev.param_a and surface.pdd.std_dev.param_b. This parametrization replaces prescribed standard deviation values over glacierized grid cells as defined by the mask variable (see geometry.ice_free_thickness_standard). Default values for the slope \(a\) and intercept \(b\) were derived from the ERA-40 reanalysis over the Greenland ice sheet [156].

The number of positive degree days is computed as the magnitude of the temperature excursion above \(0\!\phantom{|}^\circ \text{C}\) multiplied by the duration (in days) when it is above zero.

In PISM there are two methods for computing the number of positive degree days. The first computes only the expected value, by the method described in [153]. This is the default when a PDD is chosen (i.e. option -surface pdd). The second is a Monte Carlo simulation of the white noise itself, chosen by adding the option -pdd_method random_process. This Monte Carlo simulation adds the same daily variation at every point, though the seasonal cycle is (generally) location dependent. If repeatable randomness is desired use -pdd_method repeatable_random_process instead.

../_images/pdd-model-flowchart.png — Fig. 47 PISM’s positive degree day model. \(F_s\) and \(F_i\) are PDD factors for snow and ice, respectively; \(\theta_{\text{refreeze}}\) is the refreeze fraction.¶

By default, the computation summarized in Fig. 47 is performed every week. (This frequency is controlled by the parameter surface.pdd.max_evals_per_year.) To compute mass balance during each week-long time-step, PISM keeps track of the current snow depth (using units of ice-equivalent thickness). This is necessary to determine if melt should be computed using the degree day factor for snow (surface.pdd.factor_snow) or the corresponding factor for ice (surface.pdd.factor_ice).

A fraction of the melt controlled by the configuration parameter surface.pdd.refreeze (\(\theta_{\text{refreeze}}\) in Fig. 47, default: \(0.6\)) refreezes. The user can select whether melted ice should be allowed to refreeze using the configuration flag surface.pdd.refreeze_ice_melt.

Since PISM does not have a principled firn model, the snow depth is set to zero at the beginning of the balance year. See surface.mass_balance_year_start_day. Default is \(274\), corresponding to October 1^st.

Our PDD implementation is meant to be used with an atmosphere model implementing a cosine yearly cycle such as searise_greenland (section SeaRISE-Greenland), but it is not restricted to parameterizations like these.

This code also implements latitude- and mean July temperature dependent ice and snow factors using formulas (6) and (7) in [145]; set -pdd_fausto to enable. The default standard deviation of the daily variability (option -pdd_std_dev) is 2.53 degrees when -pdd_fausto is set [145]. See also configuration parameters with the prefix surface.pdd.fausto.:

T_c (272.15 Kelvin) = -1 + 273.15; for formula (6) in [145]
T_w (283.15 Kelvin) = 10 + 273.15; for formula (6) in [145]
beta_ice_c (0.015 meter / (Kelvin day)) water-equivalent thickness; for formula (6) in [145]
beta_ice_w (0.007 meter / (Kelvin day)) water-equivalent thickness; for formula (6) in [145]
beta_snow_c (0.003 meter / (Kelvin day)) water-equivalent thickness; for formula (6) in [145]
beta_snow_w (0.003 meter / (Kelvin day)) water-equivalent thickness; for formula (6) in [145]
enabled (false) Set PDD parameters using formulas (6) and (7) in [145]
latitude_beta_w (72 degree_north) latitude below which to use warm case, in formula (6) in [145]

Note that when used with periodic climate data (air temperature and precipitation) that is read from a file (see section Boundary conditions read from a file), use of time_stepping.hit_multiples is recommended: set it to the length of the climate data period in years.

This model provides the following scalar:

surface_accumulation_rate
surface_melt_rate
surface_runoff_rate

and these 2D diagnostic quantities (averaged over reporting intervals; positive flux corresponds to ice gain):

surface_accumulation_flux
surface_melt_flux
surface_runoff_flux

This makes it easy to compare the surface mass balance computed by the model to its individual components:

SMB = surface_accumulation_flux - surface_runoff_flux

Parameters

Prefix: surface.pdd.

air_temp_all_precip_as_rain (275.15 Kelvin) threshold temperature above which all precipitation is rain; must exceed surface.pdd.air_temp_all_precip_as_snow to avoid division by zero, because difference is in a denominator
air_temp_all_precip_as_snow (273.15 Kelvin) threshold temperature below which all precipitation is snow
factor_ice (0.00879121 meter / (Kelvin day)) EISMINT-Greenland value [148]; = (8 mm liquid-water-equivalent) / (pos degree day)
factor_snow (0.0032967 meter / (Kelvin day)) EISMINT-Greenland value [148]; = (3 mm liquid-water-equivalent) / (pos degree day)
fausto.T_c (272.15 Kelvin) = -1 + 273.15; for formula (6) in [145]
fausto.T_w (283.15 Kelvin) = 10 + 273.15; for formula (6) in [145]
fausto.beta_ice_c (0.015 meter / (Kelvin day)) water-equivalent thickness; for formula (6) in [145]
fausto.beta_ice_w (0.007 meter / (Kelvin day)) water-equivalent thickness; for formula (6) in [145]
fausto.beta_snow_c (0.003 meter / (Kelvin day)) water-equivalent thickness; for formula (6) in [145]
fausto.beta_snow_w (0.003 meter / (Kelvin day)) water-equivalent thickness; for formula (6) in [145]
fausto.enabled (false) Set PDD parameters using formulas (6) and (7) in [145]
fausto.latitude_beta_w (72 degree_north) latitude below which to use warm case, in formula (6) in [145]
firn_compaction_to_accumulation_ratio (0.75) How much firn as a fraction of accumulation is turned into ice
firn_depth_file The name of the file to read the firn_depth from.
interpret_precip_as_snow (no) Interpret precipitation as snow fall.
max_evals_per_year (52) maximum number of times the PDD scheme will ask for air temperature and precipitation to build location-dependent time series for computing (expected) number of positive degree days and snow accumulation; the default means the PDD uses weekly samples of the annual cycle; see also surface.pdd.std_dev.value
method (expectation_integral) PDD implementation method
positive_threshold_temp (273.15 Kelvin) temperature used to determine meaning of “positive” degree day
refreeze (0.6) EISMINT-Greenland value [148]
refreeze_ice_melt (yes) If set to “yes”, refreeze surface.pdd.refreeze fraction of melted ice, otherwise all of the melted ice runs off.
std_dev.file The name of the file to read air_temp_sd (standard deviation of air temperature) from.
std_dev.lapse_lat_base (72 degree_north) standard deviationis is a function of latitude, with value surface.pdd.std_dev.value at this latitude; this value is only active if surface.pdd.std_dev.lapse_lat_rate is nonzero
std_dev.lapse_lat_rate (0 Kelvin / degree_north) standard deviation is a function of latitude, with rate of change with respect to latitude given by this constant
std_dev.param_a (-0.15) Parameter \(a\) in \(\Sigma = aT + b\), with \(T\) in degrees C. Used only if surface.pdd.std_dev.use_param is set to yes.
std_dev.param_b (0.66 Kelvin) Parameter \(b\) in \(\Sigma = aT + b\), with \(T\) in degrees C. Used only if surface.pdd.std_dev.use_param is set to yes.
std_dev.periodic (no) If true, interpret air_temp_sd read from surface.pdd.std_dev.file as periodic in time
std_dev.use_param (no) Parameterize standard deviation as a linear function of air temperature over ice-covered grid cells. The region of application is controlled by geometry.ice_free_thickness_standard.
std_dev.value (5 Kelvin) standard deviation of daily temp variation; = EISMINT-Greenland value [148]

Diurnal Energy Balance Model “dEBM-simple”¶

Options: -surface debm_simple
Variables: surface_albedo
C++ class: pism::surface::DEBMSimple

This PISM module implements the “simple” version of the diurnal energy balance model developed by [157]. It follows [158] and includes parameterizations of the surface albedo and the atmospheric transmissivity that make it possible to run the model in a standalone, prognostic mode.

It is designed to use time-dependent forcing by near-surface air temperature and total (i.e. liquid and solid) precipitation provided by one of PISM’s atmosphere models. The temperature forcing should resolve the annual cycle, i.e. it should use monthly or more frequent temperature records if forced using Boundary conditions read from a file. In cases when only annual temperature records are available we recommend using the Cosine yearly cycle approximation.

Note

We recommend setting time_stepping.hit_multiples to the length of the climate data period in years when forcing dEBM-simple with periodic climate data that are read from a file.

Similarly to other surface models, the outputs are

ice temperature at its top surface
climatic mass balance (SMB)

dEBM-simple re-interprets near-surface air temperature to produce the ice temperature at its top surface.

SMB is defined as

\[ \begin{align}\begin{aligned}\text{SMB} &= \text{Solid accumulation} - \text{Runoff},\\\text{SMB} &= \text{Solid accumulation} - \left(\text{Melt} - \text{Refreeze} \right).\end{aligned}\end{align} \]

Solid accumulation is approximated using provided total precipitation and a linear transition from interpreting it as “all snow” when the air temperature is below surface.debm_simple.air_temp_all_precip_as_snow (default of \(0^\circ C\)) to “all rain” when the air temperature is above surface.debm_simple.air_temp_all_precip_as_rain (default of \(2^\circ C\)). Alternatively, all precipitation is interpreted as snow if surface.debm_simple.interpret_precip_as_snow is set.

Note

Part of the precipitation that is interpreted as rain is assumed to run off instantaneously and does not contribute to reported modeled runoff.

A fraction \(\theta\) (surface.debm_simple.refreeze) of computed melt amount is assumed to re-freeze:

(44)¶\[\text{Refreeze} = \theta\, \text{Melt}.\]

By default only snow melt is allowed to refreeze; set surface.debm_simple.refreeze_ice_melt to refreeze both snow and ice melt. To distinguish between melted snow and ice dEBM-simple keeps track of the evolving snow depth during a balance year and resets it to zero once a year on the day set using surface.mass_balance_year_start_day.

Let \(T\) be the near-surface air temperature and define

(45)¶\[ \begin{align}\begin{aligned}M_I &= \frac{\Delta t_{\Phi}}{\Delta t \rho_{\text{w}} L_{\text{m}}} \tau_\text{A} ( 1 - \alpha ) \bar S_{\Phi},\\M_T &= \frac{\Delta t_{\Phi}}{\Delta t \rho_{\text{w}} L_{\text{m}}} c_1 T_\text{eff},\\M_O &= \frac{\Delta t_{\Phi}}{\Delta t \rho_{\text{w}} L_{\text{m}}} c_2,\end{aligned}\end{align} \]

then the average daily melt rate is approximated by

(46)¶\[\begin{split}M &= \begin{cases} M_I + M_T + M_O,& T \ge T_{\text{min}},\\ 0,& T < T_{\text{min}}. \end{cases},\end{split}\]

Here

\(M_I\) is the insolation-driven melt contribution representing the net uptake of incoming solar shortwave radiation of the surface during the diurnal melt period,
\(M_T\) is the temperature-driven melt contribution representing the air-temperature-dependent part of the incoming longwave radiation as well as turbulent sensible heat fluxes,
\(M_O\) is the negative melt offset representing the outgoing longwave radiation and the air temperature-independent part of the incoming longwave radiation [159].

See Table 36 for details and note that (46) is equation 1 in [157].

The following two sections describe implementations of insolation-driven and temperature-driven melt contributions.

Table 36 Notation used in (45) and (46)¶
Quantity	Description
\(\Phi\)	Threshold for the solar elevation angle (`surface``.debm_simple``.phi`). It is assumed that melt can occur only when the sun is above this angle. \(\Phi\) should be treated as a tuning parameter since its value is not well constrained.
\(\Delta t_{\Phi} / \Delta t\)	Fraction of the day during which the sun is above the elevation angle \(\Phi\) and melt can occur
\(\tau_A\)	Transmissivity of the atmosphere
\(\alpha\)	Surface albedo
\(\bar S_{\Phi}\)	Mean top of the atmosphere insolation during the part of the day when the sun is above the elevation angle \(\Phi\).
\(T_{\text{eff}}\)	“Effective air temperature” computed using provided air temperature forcing and additional stochastic variations used to model the effect of daily temperature variations (see (52), [157] and [153])
\(c_1\)	Tuning parameter that controls the temperature influence on melt (`surface``.debm_simple``.c1`)
\(c_2\)	Tuning parameter that controls the (negative) melt offset (`surface``.debm_simple``.c2`)
\(\rho_w\)	Fresh water density (`constants``.fresh_water``.density`)
\(L_m\)	Latent heat of fusion (`constants``.fresh_water``.latent_heat_of_fusion`)
\(T_{\text{min}}\)	Threshold temperature (`surface``.debm_simple``.melting_threshold_temp`). Melt is prohibited if the air temperature is below \(T_{\text{min}}\) to avoid melt rates from high insolation values and low albedo values when it is too cold to actually melt.

Insolation-driven melt contribution ¶

(47)¶\[M_I = \frac{\Delta t_{\Phi}}{\Delta t \rho_{\text{w}} L_{\text{m}}} \left( \tau_\text{A} \left( 1 - \alpha \right) \bar S_{\Phi} \right),\]

This term models the influence of the mean insolation during the melt period \(\bar S_{\Phi}\), the atmosphere transmissivity \(\tau_{A}\) and the surface albedo \(\alpha\).

Mean top of the atmosphere insolation ¶

The mean top of the atmosphere insolation during the part of the day when the sun is above \(\Phi\) degrees is approximated by

(48)¶\[\bar S_{\Phi} = \frac{S_0}{h_{\Phi}} \left( \frac{\bar d}{d} \right)^2 \left( h_{\Phi} \sin(\phi) \sin(\delta) + \cos(\phi) \cos(\delta) \sin(h_{\Phi}) \right)\]

where

\(S_0\) is the solar constant surface.debm_simple.solar_constant [152],
\(\bar d / d\) is the ratio of the length of the semimajor axis of the Earth’s orbit to the Earth-Sun distance,
\(h_{\Phi}\) is the hour angle when the sun has an elevation angle of at least \(\Phi\),
\(\phi\) is the latitude
\(\delta\) is the solar declination angle.

In short, \(\bar S_{\Phi}\) is a function of latitude, the factor \(\bar d / d\), and the solar declination angle \(\delta\). In the “present day” case both \(\bar d / d\) and \(\delta\) have the period of one year and are approximated using trigonometric expansions (see [160]).

Paleo simulations

Trigonometric expansions for \(\bar d / d\) and \(\delta\) mentioned above are not applicable when modeling times far from present; in this case we use more general (and more computationally expensive) formulas ([160], chapter 2). Set surface.debm_simple.paleo.enabled to switch to using the “paleo” mode.

In this case

the ratio \(\bar d / d\) is a function of the eccentricity of the Earth’s orbit and the perihelion longitude,
the solar declination \(\delta\) is a function of the eccentricity of the Earth’s orbit, the perihelion longitude, and the Earth’s obliquity.

The values of these are set using the following configuration parameters (prefix: surface.debm_simple.paleo.):

eccentricity (0.0167) Eccentricity of the Earth’s orbit
obliquity (23.44 degree) Mean obliquity (axial tilt) of the Earth
perihelion_longitude (102.947 degree) Mean longitude of the perihelion relative to the vernal equinox

Alternatively, PISM can read in scalar time series of variables eccentricity, obliquity, and perihelion_longitude from a file specified using surface.debm_simple.paleo.file.

Note

We provide a script (examples/debm_simple/orbital_parameters.py) that can be used to generate time series of these parameters using trigonometric expansions due to [161] with corrections made by the authors of the GISS GCM ModelE (expansion coefficients used in orbital_parameters.py come from the GISS ModelE version 2.1.2). See https://data.giss.nasa.gov/modelE/ar5plots/srorbpar.html for details.

These expansions are considered to be valid for about 1 million years.

Surface albedo ¶

To capture melt processes driven by changes in albedo without requiring a more sophisticated surface process model (including the firn layer, for example), dEBM-simple assumes that the surface albedo is a piecewise linear function of the modeled melt rate.

(49)¶\[\alpha = \max\left( \alpha_{\text{max}} + \alpha_{s} \cdot M, \alpha_{\text{min}}\right).\]

Here \(M\) is the estimated melt rate from the previous time step (meters liquid water equivalent per second) and \(\alpha_s\) is a negative tuning parameter (surface.debm_simple.albedo_slope).

In this approach, the albedo decreases linearly with increasing melt from the maximum value \(\alpha_{\text{max}}\) (the “fresh snow” albedo surface.debm_simple.albedo_max) for regions with no melting to the minimum value \(\alpha_{\text{min}}\) (the “bare ice” albedo surface.debm_simple.albedo_min).

Alternatively, albedo (variable surface_albedo; no units) can be read from a file specified using surface.debm_simple.albedo_input.file.

Note

It is recommended to use monthly records of albedo in surface.debm_simple.albedo_input.file.
The fresh snow albedo is also treated as a tuning parameter. Default values of \(\alpha_{\text{max}}\) and \(\alpha_s\) were obtained by fitting this approximation to the output of a regional climate model.

Atmosphere transmissivity ¶

dEBM-simple assumes that the transmissivity of the atmosphere \(\tau_A\) is a linear function of the local surface altitude. Similarly to the albedo parameterization, the default values of coefficients \(a\) and \(b\) below were obtained using linear regression of an RCM output. This parameterization also relies on the assumption that no other processes (e.g. changing mean cloud cover in a changing climate) affect \(\tau_A\).

(50)¶\[\tau_A = a + b\cdot z,\]

where \(a\) is set by surface.debm_simple.tau_a_intercept, \(b\) is set by surface.debm_simple.tau_a_slope, and \(z\) is the ice surface altitude in meters.

Temperature-driven melt contribution ¶

(51)¶\[M_T = \frac{\Delta t_{\Phi}}{\Delta t \rho_{\text{w}} L_{\text{m}}} c_1 T_\text{eff},\]

where

(52)¶\[T_{\text{eff}}(T, \sigma) = \frac{1}{\sigma \sqrt{2 \pi}} \int_{T_{\text{pos}}}^{\infty} \xi\, \exp\left( -\frac{(\xi - T)^2}{2 \sigma^2} \right)\, d\xi.\]

The “effective temperature” \(T_{\text{eff}}\) is the expected value of “positive” excursions, i.e. excursions above the positivity threshold \(T_{\text{pos}}\) (surface.debm_simple.positive_threshold_temp, usually \(0\!\phantom{|}^\circ \text{C}\)) of stochastic temperature variations added to the provided temperature forcing.

Similarly to the PDD Temperature-index scheme, these stochastic variations are assumed to follow the normal distribution with the mean of zero and the standard deviation \(\sigma\) and are used to model the effect of daily temperature variations not resolved by this model either because of the temporal resolution of the provided forcing or the chosen time step length.

Note

The standard deviation \(\sigma\) of added daily variations should be treated as a tuning parameter. The appropriate value may change depending on the application domain (for example: Greenland vs Antarctica), the temporal resolution of the air temperature forcing and lengths of time steps taken by dEBM-simple; see surface.debm_simple.max_evals_per_year.

Here \(\sigma\) can be

constant in time and space (the default; set using surface.debm_simple.std_dev),
read from a file containing the two dimensional variable air_temp_sd that can be constant in time or time-dependent (units: Kelvin; specify the file name using surface.debm_simple.std_dev.file), or
parameterized as a function of air temperature \(T\): \(\sigma = \max(a\, (T - T_{\text{melting}}) + b, 0)\) with \(T_{\text{melting}} = 273.15\) Kelvin.

These mechanisms are controlled by parameters with the prefix surface.debm_simple.std_dev.:

file The file to read air_temp_sd (standard deviation of air temperature) from
param.a (-0.15) Parameter \(a\) in \(\sigma = \max(a\, (T - T_{\text{melting}}) + b, 0)\). Used only if surface.debm_simple.std_dev.param.enabled is set to yes.
param.b (0.66 Kelvin) Parameter \(b\) in \(\sigma = \max(a\, (T - T_{\text{melting}}) + b, 0)\). Used only if surface.debm_simple.std_dev.param.enabled is set to yes.
param.enabled (no) Parameterize standard deviation as a linear function of air temperature over ice-covered grid cells. The region of application is controlled by geometry.ice_free_thickness_standard.
periodic (no) If true, interpret forcing data as periodic in time

Tuning parameters ¶

Default values of many parameters come from [157] and are appropriate for Greenland; their values will need to change to use this model in other contexts. See Table 1 in [159] for parameter values more appropriate in an Antarctic setting and for the description of a calibration procedure that can be used to obtain some of these values.

Table 37 Notable tuning parameters in the order of decreasing importance¶
Parameter	Equation	Configuration parameters
\(c_1\)	(51)	`surface``.debm_simple``.c1`
\(c_2\)	(45)	`surface``.debm_simple``.c2`
\(T_{\text{min}}\)	(46)	`surface``.debm_simple``.melting_threshold_temp`
\(\sigma\)	(52)	`surface``.debm_simple``.std_dev` and others at the end of Temperature-driven melt contribution
\(\theta\)	(44)	`surface``.debm_simple``.refreeze`, `surface``.debm_simple``.refreeze_ice_melt`
\(\alpha_{\text{max}}\)	(49)	`surface``.debm_simple``.albedo_max`
\(\alpha_s\)	(49)	`surface``.debm_simple``.albedo_slope`
\(a\)	(50)	`surface``.debm_simple``.tau_a_intercept`
\(b\)	(50)	`surface``.debm_simple``.tau_a_slope`
\(\Phi\)	(47), (48)	`surface``.debm_simple``.phi`

PIK ¶

Options: -surface pik
Variables: climatic_mass_balance \(kg / (m^{2} s)\), lat (latitude), (degrees north)
C++ class: pism::surface::PIK

This surface model component implements the setup used in [52]. The climatic_mass_balance is read from an input (-i) file; the ice surface temperature is computed as a function of latitude (variable lat) and surface elevation (dynamically updated by PISM). See equation (1) in [52].

Scalar temperature offsets ¶

Options: -surface ...,delta_T
Variables: delta_T
C++ class: pism::surface::Delta_T

The time-dependent scalar offsets delta_T are added to ice_surface_temp computed by a surface model.

Please make sure that delta_T has the units of “Kelvin”.

This modifier is identical to the corresponding atmosphere modifier, but applies offsets at a different stage in the computation of top-surface boundary conditions needed by the ice dynamics core.

Parameters

Prefix: surface.delta_T.

file Name of the file containing temperature offsets.
periodic (no) If true, interpret forcing data as periodic in time

Adjustments using modeled change in surface elevation ¶

Options: -surface ...,elevation_change
Variables: surface_altitude (CF standard name),
C++ class: pism::surface::LapseRates

The elevation_change modifier adjusts ice-surface temperature and surface mass balance using modeled changes in surface elevation relative to a reference elevation read from a file.

The surface temperature is modified using an elevation lapse rate \(\gamma_T =\) surface.elevation_change.temperature_lapse_rate. Here

\[\gamma_T = -\frac{dT}{dz}.\]

Two methods of adjusting the SMB are available:

Scaling using an exponential factor

\[\mathrm{SMB} = \mathrm{SMB_{input}} \cdot \exp(C \cdot \Delta T),\]

where \(C =\) surface.elevation_change.smb.exp_factor and \(\Delta T\) is the temperature difference produced by applying surface.elevation_change.temperature_lapse_rate.

This mechanisms increases the SMB by \(100(\exp(C) - 1)\) percent for each degree of temperature increase.

To use this method, set -smb_adjustment scale.
Elevation lapse rate for the SMB

\[\mathrm{SMB} = \mathrm{SMB_{input}} - \Delta h \cdot \gamma_M,\]

where \(\gamma_M =\) surface.elevation_change.smb.lapse_rate and \(\Delta h\) is the difference between modeled and reference surface elevations.

To use this method, set -smb_adjustment shift.

Parameters

Prefix: surface.elevation_change..

file Name of the file containing the reference surface elevation field (variable usurf).
periodic (no) If true, interpret forcing data as periodic in time
smb.exp_factor (0 Kelvin-1) Exponential for the surface mass balance.
smb.lapse_rate (0 (m / year) / km) Lapse rate for the surface mass balance.
smb.method (shift) Choose the SMB adjustment method. scale: use temperature-change-dependent scaling factor. shift: use the SMB lapse rate.
temperature_lapse_rate (0 K / km) Lapse rate for the temperature at the top of the ice.

Mass flux adjustment ¶

Options: -surface ...,forcing
Variables: thk (ice thickness), ftt_mask (mask of zeros and ones; 1 where surface mass flux is adjusted and 0 elsewhere)
C++ class: pism::surface::ForceThickness

The forcing modifier implements a surface mass balance adjustment mechanism which forces the thickness of grounded ice to a target thickness distribution at the end of the run. The idea behind this mechanism is that spinup of ice sheet models frequently requires the surface elevation to come close to measured values at the end of a run. A simpler alternative to accomplish this, namely option -no_mass, represents an unmodeled, frequently large, violation of the mass continuity equation.

In more detail, let \(H_{\text{tar}}\) be the target thickness. Let \(H\) be the time-dependent model thickness. The surface model component described here produces the term \(M\) in the mass continuity equation:

\[\frac{\partial H}{\partial t} = M - S - \nabla\cdot \mathbf{q}.\]

(Other details of this equation do not concern us here.) The forcing modifier causes \(M\) to be adjusted by a multiple of the difference between the target thickness and the current thickness,

\[\Delta M = \alpha (H_{\text{tar}} - H)\]

where \(\alpha>0\). We are adding mass (\(\Delta M>0\)) where \(H_{\text{tar}} > H\) and ablating where \(H_{\text{tar}} < H\).

Option -force_to_thickness_file identifies the file containing the target ice thickness field thk and the mask ftt_mask. A basic run modifying surface model given would look like

pismr -i foo.nc -surface given,forcing -force_to_thickness_file bar.nc

In this case foo.nc contains fields climatic_mass_balance and ice_surface_temp, as normal for -surface given, and bar.nc contains fields thk which will serve as the target thickness and ftt_mask which defines the map plane area where this adjustment is applied. Option -force_to_thickness_alpha adjusts the value of \(\alpha\), which has a default value specified in the Configuration parameters.

In addition to this one can specify a multiplicative factor \(C\) used in areas where the target thickness field has less than -force_to_thickness_ice_free_thickness_threshold meters of ice; \(\alpha_{\text{ice free}} = C \times \alpha\). Use the -force_to_thickness_ice_free_alpha_factor option to set \(C\).

Using climate data anomalies ¶

Options: -surface ...,anomaly
Variables: ice_surface_temp_anomaly, climatic_mass_balance_anomaly \(kg / (m^{2} s)\)
C++ class: pism::surface::Anomaly

This modifier implements a spatially-variable version of -surface ...,delta_T which also applies time-dependent climatic mass balance anomalies.

See also -atmosphere ...,anomaly (section Using climate data anomalies), which is similar but applies anomalies at the atmosphere level.

Parameters

Prefix: surface.anomaly.

file Name of the file containing climate forcing fields.
periodic (no) If true, interpret forcing data as periodic in time

The caching modifier ¶

Options: -surface ...,cache
C++ class: pism::surface::Cache
See also: The caching modifier

This modifier skips surface model updates, so that a surface model is called no more than every surface.cache.update_interval 365-day “years”. A time-step of \(1\) year is used every time a surface model is updated.

This is useful in cases when inter-annual climate variability is important, but one year differs little from the next. (Coarse-grid paleo-climate runs, for example.)

Parameters

Prefix: surface.cache.

update_interval (10 365days) Update interval (in 365-day years) for the -surface cache modifier.

Preventing grounding line retreat ¶

Options: -surface ...,no_gl_retreat
C++ class: pism::surface::NoGLRetreat

This modifier adjust the surface mass balance to prevent the retreat of the grounding line. See Till friction angle optimization for an application.

Note

This modifier adds mass in violation of mass conservation. Save the diagnostic no_gl_retreat_smb_adjustment to get an idea about the amount added. Note, though, that this is an imperfect measure: it includes mass added to maintain non-negativity of ice thickness.
We assume that the sea level and the bed elevation remain constant throughout the simulation.
This does not prevent grounding line retreat caused by the thinning of the ice due to the melt at the base. Set geometry.update.use_basal_melt_rate to “false” to ensure that basal melt has no effect on the position of the grounding line