Articles | Volume 16, issue 8
Model description paper
27 Apr 2023
Model description paper |  | 27 Apr 2023

Glacier Energy and Mass Balance (GEMB): a model of firn processes for cryosphere research

Alex S. Gardner, Nicole-Jeanne Schlegel, and Eric Larour

This paper provides the first description of the open-source Glacier Energy and Mass Balance model. GEMB models the ice sheet and glacier surface–atmospheric energy and mass exchange, as well as the firn state. It is a column model (no horizontal communication) of intermediate complexity that includes those processes deemed most relevant to glacier studies. GEMB prioritizes computational efficiency to accommodate the very long (thousands of years) spin-ups necessary for initializing deep firn columns and sensitivity experiments needed to characterize model uncertainty on continental scales. The model is one-way coupled with the atmosphere, which allows the model to be run offline with a diversity of climate forcing but neglects feedback to the atmosphere. GEMB provides numerous parameterization choices for various key processes (e.g., albedo, subsurface shortwave absorption, and compaction), making it well suited for uncertainty quantification and model exploration. The model is evaluated against the current state of the art and in situ observations and is shown to perform well.

1 Introduction

The near-surface (uppermost tens to hundreds of meters of depth) energy and mass budget of mountain glaciers, ice fields, ice caps, and ice sheets (i.e., glaciers) is controlled by complex interactions between clouds, the atmospheric boundary layer, the ice surface, and processes internal to the ice–air matrix (Munro and Davies, 1977; Colbeck, 1982; Wiscombe and Warren, 1980; van den Broeke et al., 1994; Greuell and Konzelmann, 1994). It all starts with the nucleation of supercooled water vapor around impurities in the atmosphere that form highly dendritic ice crystals that become heavy and fall from the atmosphere and deposit on the glacier surface, forming a highly reflective, insulative, and low-density surface layer. Over time, ice crystals tend towards a shape that minimizes surface area through vapor diffusion and mechanical breakage, rounding the crystals, reducing reflectance (Brun, 1989; Flanner and Zender, 2016; Gardner and Sharp, 2010), and increasing both density and thermal conductivity (e.g., Calonne et al., 2019). The rate at which the metamorphism takes place depends on both the mean temperature and vertical gradients in temperature (e.g., Herron and Langway, 1980; Arthern et al., 2010). This layer will become buried by successive snowfalls, subjecting it to increasing overburden stress that causes the crystals to slide and compact, further increasing both density and conductivity as the snow transitions to firn. Sliding of crystals tends to control the rate of compaction as the firn approaches a density of ∼550 kg m−3 after which the migration of grain boundaries, through sintering processes, controls the rate of compaction to the point at which the air within the ice matrix becomes sealed off from the surrounding pore space at a density of ∼830 kg m−3 (Herron and Langway, 1980; Alley, 1987). Beyond a density of 830 kg m−3, compaction is regulated by the compression of the air within the sealed pore space.

Under most conditions, net solar radiation is the largest input of energy for melting of ice (Male and Granger, 1981). The amount of solar energy absorbed is largely governed by the effective grain size of the snow crystals within the top few centimeters of the surface, in combination with the concentration and placement of light-absorbing impurities (Gardner and Sharp, 2010; Warren and Wiscombe, 1980). This dependency can create strong positive feedbacks in the energy balance wherein increased solar input leads to enhanced grain growth that in turn results in increased shortwave absorption, modified thermal gradients, enhanced compaction, and increased thermal conductivity. The introduction of melting decreases the number of ice–air boundaries, which reduces scattering and further enhances absorption of shortwave radiation (Gardner and Sharp, 2010). If melting is sufficient to overcome capillary forces it will descend vertically within the snow–firn column, redistributing large amounts of energy though latent heat release upon freezing or mass though encountering an impermeable surface and moving horizontally within the firn (Colbeck, 1974; Coléou and Lesaffre, 1998; Marchenko et al., 2017; Hirashima et al., 2017). This complex interplay of surface processes creates nonlinear responses to changes in surface forcing that require detailed modeling of the underlying physical processes.

Modern firn modeling draws heavily on the model physics implemented within seasonal snow models that have been developed for hydrology applications and avalanche forecasting. A few of the more widely used snow models include SNOWPACK (Bartelt and Lehning, 2002), CROCUS (Brun et al., 1989), and SnowModel (Liston and Elder, 2006). For a more comprehensive discussion of processes relevant for modeling of seasonal snow, readers should refer to the provided references and those that follow. Here we review those aspects relevant to modeling of perennial firn over ice sheets.

Early numerical modeling of firn was motivated by ice core research with an emphasis on understanding borehole temperature and pore close-off that is needed to determine the ice-age gas-age differential (e.g., Herron and Langway, 1980; Greuell and Oerlemans, 1989). Since then, firn models have become critical to the estimation of the surface mass budget of the ice sheets (e.g., Pfeffer et al., 1991; Janssens and Huybrechts, 2000) as they are needed to model the vertical movement of meltwater within the snow and firn and to determine if the meltwater refreezes in place or moves horizontally following the hydrologic gradient (Marsh and Woo, 1984; Pfeffer et al., 1990). Subsequent firn modeling efforts suggest that a warming climate will result in a steady decrease in the Greenland Ice Sheet's capacity to retain meltwater due to a reduction in firn pore space (van Angelen et al., 2013), a finding that has been supported by in situ observations (Vandecrux et al., 2019). Further, impermeable ice layers formed at depth within the firn can support perched aquifers or enhance meltwater runoff (Culberg et al., 2021; Miller et al., 2022; Macferrin et al., 2019). Recent observations of extensive “firn aquifers” in Greenland that persist throughout the winter, when snow accumulation and melt rates are high (Forster et al., 2014), have proven the utility of firn modeling for understanding newly discovered phenomena (Steger et al., 2017).

With the launch of the first satellite altimeters in the late 1970s it became possible, for the first time, to measure large-scale changes in ice sheet topography (Zwally et al., 1989). However, the interpretation of such changes was challenging due to unknown changes in near-surface snow and firn density (Van Der Veen, 1993). The challenge is that changes in the near-surface density (i.e., changes in the firn air content: FAC) can cause changes in surface elevation without any corresponding changes in mass. By the late 1990s, numerical firn models were being used to estimate the uncertainties in the conversion of elevation change to mass change (Arthern and Wingham, 1998) and by the mid-2000s were being used to estimate the ice-sheet-wide changes in FAC for ice-sheet-wide estimation of mass change from satellite altimetry data (Zwally et al., 2005). Despite significant advances in firn modeling over the past decade and a half (e.g., Bougamont and Bamber, 2005; Ligtenberg et al., 2011), models still have significant deficits (Lundin et al., 2017; Vandecrux et al., 2020) and remain the largest source of uncertainty when inferring ice sheet mass change from satellite altimetry (Smith et al., 2020).

Here we present a new coupled surface energy balance and firn model, the Glacier Energy and Mass Balance model (GEMB, the “B” is silent), which has been integrated as a module into the Ice-sheet and Sea-level System Model (ISSM). Here we describe the state of GEMB as of version 1.0. ISSM is an open-source software framework for modeling ice sheets, solid earth, and sea level response that is developed at NASA's Jet Propulsion Laboratory (JPL) in conjunction with the University of California, Irvine, and Dartmouth College (Larour et al., 2012c). The Dakota software (, last access: 15 May 2015) is embedded within the ISSM framework, facilitating uncertainty quantification and sensitivity studies (Larour et al., 2012a, b; Schlegel et al., 2013, 2015; Schlegel and Larour, 2019). Currently GEMB is a stand-alone module that is responsible for the calculation of ice sheet and glacier surface–atmospheric energy and mass exchange as well as firn state within ISSM. GEMB provides the ice sheet flow model with near-surface ice temperature and mass flux boundary conditions. It is a column model (no horizontal communication) of intermediate complexity that includes those processes deemed most relevant to glacier studies. GEMB prioritizes computational efficiency to accommodate the very long (thousands of years) spin-ups necessary for initializing deep firn columns and sensitivity experiments needed to characterize model uncertainty on continental scales. GEMB is not coupled with a model of the atmosphere and instead runs offline, forced with climate reanalysis or climate model data. This approach allows flexibility in selection of forcing datasets at the expense of simulating surface–atmosphere feedbacks.

GEMB is currently run independently of the ice flow model. The goal is to eventually couple the ice flow model with GEMB. In this situation the flow model would provide the 3-D surface displacement vectors to GEMB to inform the advection of the column nodes relative to the climate forcing and to allow for solving of longitudinal stretching of the firn. Currently the major hurdle for coupling the models is the implementation of downscaling routines for the climate forcing.

2 Model description

GEMB is a vertical 1-D column model, i.e., no horizontal communication between model nodes, that simulates atmosphere–surface mass and energy exchanges and the internal evolution of snow, firn, and ice. The model shares many characteristics with earlier published firn models that also simulate atmosphere–surface exchanges (e.g., Bassford, 2002; Bougamont and Bamber, 2005; Greuell and Konzelmann, 1994). The model is a finite-difference model with tens to hundreds of layers, the thickness of which is managed dynamically. It is forced at its surface with near-surface (2–10 m) estimates of precipitation, air temperature, wind speed, vapor pressure, surface pressure, downwelling longwave and shortwave radiation fluxes, and optional inputs of solar zenith angle, cloud optical thickness, and bare ice albedo. At its bottom boundary, the model applies a constant thermal flux. Internally, the model simulates thermal diffusion, shortwave subsurface penetration, meltwater retention, percolation and refreeze, effective snow grain size, dendricity, sphericity, and compaction. The model does not yet account for changes in firn due to horizontal advection or ice divergence (Horlings et al., 2021). GEMB also does not account for changes in thermal properties caused by debris cover that are important for the modeling of valley glaciers. In this section we detail specific implementation of various processes and their options within the model.

Figure 1Diagram showing the model layer initialization as described in Sect. 2.1: dz is the layer thickness, dztop is the maximum near-surface layer thickness, ztop is the depth of the near-surface, zmax is the maximum column depth, β is a unitless scaling by depth parameter, and n is the layer number.


2.1 Layer initialization

GEMB is a finite-difference model that simulates the snow–firn–ice as a number of discrete plane-parallel horizontal layers, each with their own state parameters. The thickness (dz) of each horizontal layer is initialized according to default or user-supplied specification of the minimum and maximum thickness. Firn properties and energy fluxes have more vertical heterogeneity nearer the surface–atmosphere boundary than deeper within the firn column. Because of this the near-surface firn must be modeled using finer model layers than those deeper within the firn column, where energy fluxes are small and gradients in ice properties are more diffuse. To accommodate this, finer model layers are specified for the near-surface. Layers below the near-surface are assigned increasing layer thickness with depth. Layer thickness is increased with depth to minimize the overall number of layers needed to simulate the firn column, improving computational efficiency without sacrificing accuracy. Users can specify the maximum near-surface thickness (dztop, default = 0.05 m), the depth of the near-surface (ztop, default = 10 m), the maximum column depth (zmax, default = 250 m), and a unitless scaling by depth parameter (β, default = 1.025). The thickness of each layer and the depth of the snow-firn-ice column are then determined as shown in Fig. 1. The user can also specify a minimum near-surface thickness (dzmin) and a minimum column depth (zmin). GEMB will combine layers, with the layer directly below, if they are thinner than dzmin when located within ztop and dzminβ(nn_top) when located below ztop. GEMB will split layers if they are thicker than dzmax when located within ztop and dzmaxβ(nn_top) when located below ztop. If the total depth of the column is less than zmin then an additional ice layer is added to the bottom of the column.

2.2 Grain growth

The model's main time step (i.e., the time step for all processes excluding the thermal model) is set by the time step of the input data (typically on the order of hours). For all time steps, the first calculation of the GEMB model is done by the grain growth module. This module tracks and evolves the effective grain radius, dendricity, and sphericity of all model layers through time. These properties evolve according to published laboratory estimates for dendritic (Brun et al., 1992, Eqs. 3 and 4), non-dendritic (Marbouty, 1980, Fig. 9), and wet snow metamorphism (Brun, 1989, Fig. 6) that are dependent on mean layer temperature and density and on gradients in temperature.

2.3 Shortwave flux

After the grain properties of all model layers are determined, GEMB runs the albedo module. GEMB provides five methods for calculating broadband surface albedo. The default albedo scheme is after Gardner and Sharp (2010), where albedo is calculated as a function of grain specific surface area (S^), concentration of light-absorbing carbon (c: optional), solar zenith angle (u: optional), and cloud optical thickness (τ: optional):

(1) α = α S ^ + d α c + d α u + d α τ ,

where αS^ is the pure snow albedo, dαc is the change in albedo due to the presence of light-absorbing carbon, dαu is the change in albedo due changes in the solar zenith angle, and dατ is the change in albedo due to changes in the cloud optical thickness. Equations for each of the four components that sum to the net broadband albedo are provided in Eqs. (7)–(11) in Gardner and Sharp (2010).

Surface broadband albedo can also be calculated based on the effective grain radius (re), where broadband albedo is determined as a summation of the albedo within three spectral bands of solar irradiance (Brun et al., 1992; Lefebre et al., 2003). Effective grain radius can be related to specific surface area as

(2) r e = 3 ρ i S ^ ,

where ρi is the density of ice (910 kg m−3).

Additionally, the surface albedo can be parameterized as a function of snow density and cloud amount according to Greuell and Konzelmann (1994) or as a combination of exponential time decay and firn wetness following Bougamont and Bamber (2005). A detailed review of all mentioned albedo schemes can be found in Sect. 4.2 of Gardner and Sharp (2010).

The albedo of bare ice (defined by a density threshold, default of 820 kg m−3) can be set as a constant everywhere or spatially varying when specified for each model node (e.g., as derived from MODIS). Alternatively, the ice albedo and shallow snow-covered (<10 cm in depth) ice albedo can be parameterized as a function of model-estimated accumulation of surface meltwater following Eq. (5) of Alexander et al. (2014).

When shortwave radiation reaches the glacier surface it is scattered and absorbed. The fraction of energy that is scattered then reflected back to the atmosphere is dictated by the modeled broadband albedo. The remainder of the energy is absorbed within the snow, firn, and ice. By default, GEMB will allocate all absorbed energy to the top model layer (subsurface absorption turned off) but also allows the absorbed radiation to be distributed within the near-surface (i.e., across multiple near-surface layers: subsurface absorption turned on). When the subsurface absorption is turned on, the amount of shortwave radiation absorbed within each model layer is dependent on which albedo scheme is used. If the albedo scheme is based on effective grain radius, then the subsurface absorption is calculated in three spectral bands, dependent on effective grain radius as described by Brun et al. (1989, p. 333–334). In all other cases, the subsurface albedo is treated as a function of each layer's snow density (Bassford, 2002, Eq. 4.15). In this case a specified fraction of the shortwave radiation (approximately 36 %) is absorbed by the surface layer (Greuell and Konzelmann, 1994), with the rest of the energy being absorbed at depth with consideration of Beer's law (Bassford, 2002).

2.4 Thermodynamics

The thermodynamics module is responsible for determining the temperature of the snow, firn, and ice. Temperature (T) evolves according to thermal diffusion in response to radiative, sensible, and latent heat fluxes. The thermal conductivity of snow and ice is calculated according to Sturm et al. (1997) (default) or Calonne et al. (2011). GEMB calculates thermal diffusion using an explicit forward-in-time, central-in-space method (Patankar, 1980, chap. 3 and 4):

(3) T p = A u T u + A d T d + A p - A u - A d T p + S A p ,

where neighbor coefficients Au, Ap, and Ad are


subscripts “u” and “d” represent grid points up and down from the center grid point “p”, and identifies previous time step values. S is a source term. dz is layer thickness, ρ is layer density, dt is the time step, K is effective thermal conductivity, and CI is the heat capacity of ice.

The thermal diffusion calculation is executed at a much finer time step (dt) than the main time step: tens of seconds versus hours. The finer time step is required to satisfy the Courant–Friedrichs–Lewy condition, a necessary condition to ensure numeric stability of the thermal solution when using an explicit forward-in-time method (Courant et al., 1928). See Patankar (1980) for an overview of numerical heat transfer and stability criterion for explicit forward methods. The small time step for the thermal diffusion makes this module, by far, the most computationally expensive component of GEMB. The thermal module is looped within each GEMB time step, and therefore the thermal time step must divide evenly into the model's main time step. The effective thermal conductivity (keff) can be calculated as a function of density according to Sturm et al. (1997, Eq. 4) for densities <910 kg m−3:

(4) k eff = 0.138 - 1.01 × 10 - 3 ρ + 3.233 × 10 - 6 ρ 2 ,

where ρ is the layer density (kg m−3).

The effective thermal conductivity can also be calculated according to Calonne et al. (2011).

(5) k eff = 0.024 - 1.23 × 10 - 4 ρ + 2.5 × 10 - 6 ρ 2

For densities ≥910 kg m−3 effective thermal conductivity is calculated as a function of temperature:

(6) k eff = 9.828 e - 5.7 × 10 - 3 T ,

where T is the layer temperature in Kelvin.

The maximum acceptable thermal time step is calculated, dependent upon the thermal conductivity, and then divided by a scaling factor to achieve numerical stability. The single maximum acceptable thermal time step is calculated for each GEMB time step. Within every time step, the thermal module calculates the diffusion of temperature within the snow, firn, and ice.

At each thermal diffusion time step the module determines the radiative and turbulent fluxes (Paterson, 1994; Murphy and Koop, 2005), rates of evaporation and condensation, and the new temperature profile throughout the model depth. Layer temperatures are allowed to artificially exceed the melting point of ice (273.15 K), the thermal energy above which is later converted to melt by the melt module. That said, the surface temperature used for computing the outgoing longwave energy and turbulent heat flux is taken as the lesser of the temperature of the top model layer or 273.15 K. The snow–ice is assumed to be a grey body (default emissivity = 0.97) with the incoming longwave radiation absorbed within the top layer and outgoing longwave radiation calculated using the Stefan–Boltzmann constant (5.67×10-8 W m−2 K−4). Calculation of turbulent fluxes is defined in the next subsection. At the end of the thermodynamics module, any mass lost or added by evaporation or condensation is removed or added to the top layer.

GEMB uses a Lagrangian framework. As such the advective heat transport by firn and ice is handled implicitly. Redistribution of mass and energy by vertical movement of meltwater within the firn matrix is handled using a bucket scheme and is described in Sect. 2.7.

2.5 Turbulent flux

Turbulent fluxes as well as evaporation, sublimation, and condensation are computed within the thermodynamics module using the bulk method. Turbulent fluxes are dependent on the atmospheric stratification (i.e., stability), which can be quantified through, for example, the bulk Richardson number, Ri:

(7) Ri = 2 g z v - z 0 T a - T s 1 × 10 5 p air 0.286 z T T a + T s v z v 2 ,

where g is the acceleration due to gravity (9.81 m s−2), zv is the height in meters at which the wind speed (v) was measured, z0 is the surface roughness in meters, Ta is the air temperature at zv, and Ts is the surface temperature, with the latter two in Kelvin. pair is the air pressure in pascals. This is the bulk version of Högström (1996: Eq. 9).

The bulk transfer coefficient is calculated based on momentum roughness length and the height at which air temperature and wind speed inputs are provided. The momentum roughness length is set to 0.12 mm for dry snow, 1.3 mm for wet snow, and 3.2 mm for ice, after Bougamont and Bamber (2005), and the humidity and temperature roughness lengths are set to 10 % of these values (Foken, 2008). Within each thermodynamics time step, the model calculates the bulk Richardson number (Ri) and determines the Monin–Obukhov stability parameters for weighting the sensible and latent heat flux. In stable conditions (i.e., Ri>0), the stability weighting parameters are determined after Beljaars and Holtslag (1991) and Ding et al. (2020, Eqs. 12 and 13), and in unstable conditions (i.e., Ri≤0), the stability weighting parameters are determined after Högström (1996) as in Sjöblom (2014, Eqs. 26 and 27). To determine the latent heat flux, the model calculates the surface vapor pressure for liquid water or ice (Murphy and Koop, 2005), depending on whether the surface is wet (Murray, 1967) or dry (Bolton, 1980), respectively. The sensible heat and latent heat are then determined as a function of the bulk transfer coefficient and the Monin–Obukhov stability parameters, as well as model inputs of surface temperature (lesser of the temperature of the top model layer or 273.15 K), near-surface air temperature, near-surface wind speed, near-surface vapor pressure, and the air pressure at the surface. Mass transfer due to evaporation, sublimation, and condensation is based on the latent heat flux along with the latent heat of evaporation and condensation (gas  liquid) or the latent heat of sublimation (solid  gas).

2.6 Accumulation

The GEMB accumulation module is responsible for the addition of mass due to precipitation and related modification of uppermost layer properties (i.e., density, temperature, dendricity, sphericity, and grain radius). When precipitation occurs and the air temperature is below the melting point, the model accumulates precipitation as snow at the surface. The density of newly fallen (ρ0) snow can be set to a default value, i.e., 350 kg m−3 (Weinhart et al., 2020) for Antarctica or 315 kg m−3 (Fausto et al., 2018) for Greenland, or it can be parameterized as a function of the annual temperature, accumulation, surface pressure, and wind speed according to Kaspers et al. (2004):

(8) ρ 0 = 7.36 × 10 - 2 + 1.06 × 10 - 3 T s + 6.69 × 10 - 2 A + 4.77 × 10 - 3 v ,

where Ts is the annual surface temperature in Kelvin, A is the average annual accumulation (m w.e.), and v is the annual mean wind speed (m s−1) measured 10 m above the surface.

ρ0 can also be modeled as a function of annual air temperature according to Kuipers Munneke et al. (2015).

(9) ρ 0 = 481 + 4.834 ( T s - 273.15 )

By default, fresh snow dendricity is set to 1 and fresh snow sphericity is set to 0.5. For Greenland, we consider wind effects on the snow's initial dendricity and sphericity following Guyomarc'h and Merindol (1998) as in Vionnet et al. (2012). In the case in which the air temperature is at or exceeds the melting point, precipitation is treated as rain, and liquid water is added to the uppermost model layer. When this occurs, the surface layer properties, specifically temperature and density, are updated to account for the addition of mass as ice and the release of latent heat. Any thermal energy exceeding the melting point of ice will be converted back to liquid water by the melt module and allowed to percolate. The approach of adding rain as ice plus latent heat of refreeze, and later converting back to liquid water if the top firn layer does not possess enough cold content to freeze the rain in place, is simply for computational convenience.

If the newly accumulated mass exceeds the minimum allowable top-layer thickness (dzmin), the accumulation is added as a new surface layer. Otherwise, the snow is added to the existing surface layer, and the surface properties are adjusted to accommodate the accumulation by averaging properties, weighting by mass, with the exception of dendricity and sphericity, which take the value of the newly accumulated snow. New snow that is accumulated is assumed to have the same temperature as the near-surface air temperature. GEMB does not account for mass transport due to wind.

2.7 Melt

After new mass is accumulated within the top layer(s), the melt module is run. This module is responsible for calculating how much melt will occur throughout the model column, how much of the melt will percolate into layers below, how much of this melt will refreeze, how much of the melt will run off, and how the temperature of each layer will evolve accordingly. GEMB uses a bucket scheme (Steger et al., 2017) to model the vertical movement of meltwater within the firn. The first step of this scheme is to determine how much of the current pore water in the firn column can be refrozen without heating the firn layer above the freezing point of ice. If this water does exist, it is refrozen locally, and the model layers' physical and thermal properties are updated in response to this process. The next step is to determine where melt will occur. Beginning with the surface layer, the module determines if the local thermal energy is capable of melting any of the ice within the layer, and if it is, this portion of the layer is melted and any excess thermal energy (i.e., all ice in layer is melted) or meltwater exceeding the irreducible water content (i.e., water held in place by capillary effects) is redistributed to the layer below. The irreducible water content is assumed to be 7 % of the pore volume according to Colbeck (1973). As the meltwater is distributed to lower layers, it may reach a layer of impermeability (density of pore hole close-off, ∼830 kg m−3), and at this point, that meltwater exits the system and is considered runoff. If the layer into which the meltwater flows is permeable and the capillary effects can retain the incoming meltwater, the water is held within that layer and is combined with any pre-existing liquid water. If any portion of the incoming meltwater cannot be accommodated within the layer, it percolates (instantaneously) to the layer below. Water that remains is refrozen locally, within each layer, but only until the local temperature reaches the freezing point (273.15 K). Meltwater that reaches the deepest layer of the column exits the system and is considered by the model to be runoff. Runoff values from all layers of a column are summed to determine the total amount of melt that runs off (i.e., exits the model domain). GEMB does not currently model firn aquifers as any water exceeding the irreducible water content is removed from the system once it encounters an impermeable ice layer. Modifying GEMB to simulate firn aquifers would simply require the addition of a slope-dependent runoff criterion.

2.8 Layer management

After the melt, refreeze, and runoff calculations are complete, the melt module ensures that the thickness of any single layer does not exceed thresholds set for the minimum and maximum allowable layer thickness. This is done through the merging or splitting of model layers and associated changes in layer properties. Specifically, this module ensures that layering within the column adheres to the user-defined minimum thickness for a layer (dzmin), the user-defined maximum change in depth between adjacent layers (β, default 1.025 % or 2.5 % change between layers), and the user-defined maximum (zmax) and minimum (zmin) depths of the total column (with default values of 250 and 130 m, respectively). GEMB will combine layers, with the layer directly below, if they are thinner than dzmin when located within ztop and dzminβ(nn_top) when located below ztop. GEMB will split layers if they are thicker than dzmax when located within ztop and dzmaxβ(nn_top) when located below ztop.

Before completion, the melt module checks for the conservation of mass and energy and throws an error if these values are not conserved. If the maximum column depth is exceeded, mass is removed from the bottom layer. If the depth of the model falls below the minimum allowable depth, then ice is added to the bottom layer of the model with identical properties as possessed by the bottom layer. All additions and subtractions of mass are cataloged and accounted for in final mass change estimates (i.e., not included in estimates of surface mass balance).

2.9 Compaction

After the merging or splitting of model layers, GEMB runs the densification module. This module calculates how much the snow and firn layers compact (increase in density) over the main time step. Compaction (dρ) is determined following Herron and Langway (1980):

(10) d ρ = c ( ρ ice - ρ ) d t ,

where dρ is the change in layer density (kg m−3), c is the rate parameter, ρice is the density of ice (910 kg m−3), ρ is the layer density from the previous time step, and dt is the time interval. c=c0 for density ≤550 kg m−3 and c=c1 for densities between 550 and 830 kg m−3.

The densification module offers five different approaches for computing compaction rate factors c0 and c1.

  1. Herron and Langway (1980)

    (11) c 0 = 11 e - 10 160 T R b ˙ ρ water c 1 = 575 e - 21 400 T R b ˙ ρ water

    T is the layer temperature in Kelvin, R is the gas constant (8.314 J K−1 mol−1), b˙ is the average annual accumulation (kg m−2 yr−1), and ρwater is the density of water (1000 kg m−3).

  2. Arthern et al. (2010) (semi-empirical)

    (12) c 0 = 0.07 e - E c T R + E g T mean R b ˙ g c 1 = 0.03 e - E c T R + E g T mean R b ˙ g

    Ec is the activation energy for self-diffusion of water molecules through the ice lattice (60 kJ mol−1), Eg is the activation energy for grain growth (42.4 kJ mol−1), Tmean is the mean annual temperature in Kelvin, and g is the acceleration due to gravity (9.81 m s−1).

  3. Arthern the al. (2010) (physical model)

    (13) c 0 = 9.2 × 10 - 9 e - E c T R σ r e 2 c 1 = 3.7 × 10 - 9 e - E c T R σ r e 2

    re is the effective grain radius in meters (see Eq. 2) and σ is the overburden pressure.

  4. Li and Zwally (2004)

    (14) c 0 = c 1 = b ˙ ρ ice 139.21 - 0.542 T mean 8.36 ( 273.15 - T ) - 2.061
  5. Helsen et al. (2008)

    (15) c 0 = c 1 = b ˙ ρ ice 76.138 - 0.28965 T mean 8.36 ( 273.15 - T ) - 2.061

In addition to the densification models, GEMB implements a calibration of the Arthern et al. (2010) model that was developed by Ligtenberg et al. (2011) for the Antarctic and applied by Kuipers Munneke et al. (2015) to Greenland. Using the semi-empirical model for dry snow compaction, as described by Arthern et al. (2010) in Appendix B, the c0 and c1 rate parameters are scaled by model-to-observed calibration values (MO550 for density ≤550 kg m−3 and MO830 for densities between 550 and 830 kg m−3) that are trained by a comparison between firn core density profiles and modeled density profiles with respect to a climatological mean accumulation rate (b˙, kg m−2 yr−1):

(16) MO 550 / 830 = b 550 / 830 + m 550 / 830 ln ( b ˙ ) ,

where 550 and 830 indicate coefficients for modeling of densities ≤550 kg m−3 and between 550 and 830 kg m−3, respectively.

Offset (b550/830) and scale (m550/830) coefficients are estimated by spinning up uncalibrated GEMB firn profiles at node locations that are closest to the location of firn cores. Offset and scale coefficients are then determined by minimizing the least-squares fit between modeled and observed densities. These calibration coefficients are then applied in the ice-sheet-wide simulations.

Table 1GEMB model parameters used in this study.

Download Print Version | Download XLSX

Figure 2GEMB Greenland (a) and Antarctic (b) model node locations used for this study.

3 Study-specific model setup

3.1 Parameterization selection

GEMB v1.0 has many options for various model parameterizations. For this study we use a simplified version of the Gardner and Sharp (2010) albedo scheme in which albedo is modeled solely as a function of snow–ice grain specific surface area (S^, surface area per unit mass). Subsurface shortwave penetration is turned off. Bare ice and shallow snow-covered ice albedo (areas with a surface density >820 kg m−3) are spatially varying and derived from the MODIS MCD43C3 16 d black-sky albedo product (Schaaf and Wang, 2015). Here, we determine the bare ice albedo by taking all summer (JJA) MCD43C3 16 d black-sky albedo values and calculating the lowest 10th percentile of measured albedos for every pixel. Bare ice albedos are not allowed to be less than 0.4. Albedos are then smoothed to a 0.25 (∼32 km) resolution and bilinearly interpolated onto the GEMB model nodes. The smoothing is preformed to match the resolution of the ERA5 climate reanalysis data (Copernicus Climate Change Service, 2019), which is the forcing of choice when not comparing to the Institute for Marine and Atmospheric Research Utrecht Firn Densification Model (IMAU-FDM). The full list of selected model parameters used for this study is provided in Table 1.

3.2 Model grid

GEMB has no horizontal communication and thus easily supports nonuniform node spacing for computational efficiency; i.e., the model does not need to be run at the same spatial resolution as the climate forcing. This allows the model to use a coarser node spacing for areas with small spatial gradients in surface forcing (e.g., flat ice sheet interiors) and to use refined node spacing in areas of steep spatial gradients (e.g., areas of complex topography or near ice–ocean boundaries). For the Greenland simulations presented here, we use 10 990 nodes with node spacing ranging from 0.7 km along the coast to 21.1 km in the interior. For the Antarctic we simulate 50 390 nodes with node spacing ranging from 3.7 km along the periphery to 33.3 km in the interior. Node locations for both ice sheets are shown in Fig. 2.

3.3 Atmospheric forcing

For this study we force GEMB with climate data from the regional atmospheric climate model RACMO2.3p2 for the Antarctic (van Wessem et al., 2018) and RACMO2.3 for Greenland (Noël et al., 2016). Atmospheric fields of 10 m wind speed, 2 m temperature, surface pressure, incoming longwave and shortwave radiation, vapor pressure, and precipitation are supplied at 3 h resolution. The Antarctic product is provided at 27 km horizontal resolution for the years 1979 through 2014 and Greenland data at 11 km resolution for the years 1979 through 2014.

Figure 3Locations of shallow (green) and deep (blue) firn cores. Firn cores are classified as shallow if they reach the 550 kg m−3 density horizon and deep if they extend to the 830 kg m−3 density horizon.

3.4 Firn model calibration

GEMB is run with the uncalibrated semi-empirical Arthern et al. (2010) firn model (Eq. 12) and is forced with a repeated cycle of historical climatology until a steady-state density profile is reached, e.g. 5000 years for Greenland and 6000 years for Antarctica. The depths of the modeled density horizons are extracted at the end of this spin-up procedure from the nodes that are closest in location to the firn cores used for calibration (Fig. 3). Here we use 64 shallow firn cores and 15 deep firn cores for calibration of Greenland firn and 117 shallow and 29 deep firn cores for calibration of Antarctic firn (Montgomery et al., 2018; Smith et al., 2020; Medley et al., 2022). Firn cores are classified as shallow if they reach beyond the 550 kg m−3 horizon up to the 830 kg m−3 horizon and deep if they reach beyond the 830 kg m−3 density horizon. Appendix A provides the full list of cores used in the calibration and validation along with the shallow or deep classification and citation.

For these firn core locations, MO550 and MO830 are plotted independently as functions of the natural log of b˙ (see Eq. 16). A line is then fit through the points for each, representing the linear expression of how the model-to-observed depths vary as a function of b˙ at each horizon. Next, the model steady-state spin-up is repeated, but this time multiplying the rate coefficients of c0 and c1 of the semi-empirical model of Arthern et al. (2010) (Eq. 12) by MO550 and MO830, respectively. Following Ligtenberg et al. (2011), MO550 and MO830 are allowed a minimum value of 0.25. Derived calibration offset (b550/830) and scale (m550/830) coefficients (Eq. 16) are provided in Table 1.

The new steady-state model density profiles are then verified against the calibration subset of firn cores as well as an additional withheld subset. Figure 4 shows the modeled vs. observed 550 and 830 kg m−3 depths for the uncalibrated and calibrated model runs.

Figure 4Comparison between modeled and observed 550 kg m−3 (a, b) and 830 kg m−3 (c, d) depths for the Greenland (a, c) and Antarctic (b, d) Ice Sheets for the locations of firn cores shown in Fig. 3. Blue markers show pre-calibrated model agreement with observations, and red markers show the post-calibrated model fit. The one-to-one line is shown by the black diagonal lines. Pre- and post-model calibration fits to observations (root mean square error: RMSE) are shown in the top left of each panel.


4 Comparison to other models

4.1 Surface mass balance

Here we compare the GEMB v1.0 surface mass balance (SMB) components to those computed within the RACMO surface model: RACMO2.3p2 for the Antarctic (van Wessem et al., 2018) and RACMO2.3 for Greenland (Noël et al., 2016). RACMO values are bilinearly interpolated to GEMB nodes. Since both GEMB and RACMO surfaces are forced with the same atmospheric data (RACMO), differences between modeled components can be directly attributed to differences in how surface processes are treated (e.g., albedo, roughness lengths, thermal diffusion). We also note that the RACMO surface model accounts for feedbacks between the atmosphere and the surface, while GEMB does not.

Figure 5Mean 1979–2014 Greenland Ice Sheet GEMB (a, d, g) and RACMO (b, e, h) model output in meters of ice equivalent and their difference (c, f, i).

Figure 6Mean 1979–2014 Greenland Ice Sheet GEMB (a) and IMAU-FDM (b) firn air content and their difference (c).

Figure 7Mean 1979–2014 Antarctic Ice Sheet GEMB (left) and RACMO/IMAU-FDM (middle) model output in meters of ice equivalent and their difference (right). FDM is short for IMAU-FDM, the results of which are shown in the bottom row.

Mean spatial patterns of surface mass balance (i.e., accumulation (precipitation + deposition)  ablation (runoff + evaporation + sublimation)) and firn air content are shown for Greenland in Figs. 5 and 6, respectively, and for the Antarctic Ice Sheet in Fig. 7. The SMB patterns for the Antarctic Ice Sheet are driven largely by snowfall and sublimation, with little contribution from meltwater runoff. In Greenland, snowfall and runoff dominate. This largely explains the much closer agreement between models in the Antarctic vs. Greenland.

In Greenland the largest differences between models are concentrated in areas of high melt and low elevation that comprise the ice sheet periphery. In general, GEMB produces more negative surface mass balance in areas of maximum melt and slightly less negative surface mass balance near the equilibrium altitude. The more negative surface mass balance at lower elevations is most likely caused by a lower “bare ice” albedo in GEMB than in RACMO. Higher surface mass balance near the equilibrium line altitude is likely due to lower rates of fresh snowmelt. For the Antarctic, differences can largely be summarized as GEMB having higher rates of melt and runoff along the ice sheet periphery and lower rates of sublimation in areas conducive to katabatic winds (Parish and Bromwich, 1987). GEMB's higher rates of melt are most likely driven by lower surface albedo. Lower rates of sublimation can be attributed to the fact that GEMB does not yet (as of v1.0) include a model for snowdrift sublimation, while RACMO does (Lenaerts et al., 2010).

Figure 8The 1979–2014 Greenland Ice Sheet GEMB (blue) and RACMO (red) monthly (a) and cumulative (b) model output. GEMB (blue) and IMAU-FDM (red) average fin air content anomaly, relative to 1979, is shown in panel (c). E/C is the mass change due to evaporation, condensation, sublimation, and deposition.


Figure 9The 1979–2014 Antarctic Ice Sheet GEMB (blue) and RACMO (red) monthly (a) and cumulative (b) model output. GEMB (blue) and IMAU-FDM (red) average fin air content anomaly, relative to 1979, is shown in panel (c). E/C is the mass change due to evaporation, condensation, sublimation, and deposition.


Looking at Greenland (Fig. 8) and Antarctic (Fig. 9) monthly (a) and cumulative (b) time series of surface mass balance components we find very close agreement between GEMB and RACMO. Agreement is so close that several RACMO variables are not visible in the monthly time series (Figs. 8a and 9a) as they are overlaid by nearly identical GEMB output. The most notable difference is higher rates of melt for RACMO simulations in Greenland that are largely compensated for by increased meltwater retention within the firn. In the Antarctic lower rates of sublimation and higher rates of meltwater runoff in GEMB result in a slightly more positive surface mass balance trend relative to those simulated by RACMO.

4.2 Firn air content

Next, we compare firn properties as modeled by GEMB v1.0 to those modeled by IMAU-FDM (Ligtenberg et al., 2011) for Greenland (Ligtenberg et al., 2018) and the Antarctic (Ligtenberg et al., 2011) forced with the same RACMO data as used by the GEMB simulations (see the “Study-specific model setup” section and Table 1). Like GEMB, IMAU-FDM is an uncoupled firn model (i.e., not coupled with the atmospheric model). IMAU-FDM is a widely used firn model product. We note that there has been a recent and rapid increase in the number of ice-sheet-wide firn modeling studies (Brils et al., 2022; Dunmire et al., 2021; Keenan et al., 2021; Medley et al., 2022; Veldhuijsen et al., 2022). An in-depth comparison between all modeled results would be highly valuable but is not done here. We first look at the spatial patterns of total firn air content (FAC), which is the depth of air in meters per unit area. For Greenland, GEMB tends to generate lower FAC relative to IMAU-FDM in the west and southeast percolation zones along the ice sheet margins and a slightly higher FAC in the northeast (Fig. 6). This is counterintuitive as GEMB is shown to produce less melt than RACMO in areas where RACMO has a higher FAC; we would expect the opposite to be true. This difference is likely the result of a slightly more aggressive compaction calibration scaling coefficient used by GEMB (Table 1: b550/830=1.27/2.00, m550/830=-0.12/-0.25) relative to those used by IMAU-FDM (b550/830=1.04/1.73, m550/830=-0.09/-0.20). Differences in calibration coefficients can be attributed to differences in the firn cores used for the calibration, as IMAU-FDM had access to fewer cores when calibrated, and differences in the firn models themselves.

In the Antarctic there is generally better agreement between models and scaling coefficients of b550/830=1.64/2.00 and m550/830=-0.17/-0.24 for GEMB and b550/830=1.44/2.37 and m550/830=-0.15/-0.29 for IMAU-FDM (Kuipers Munneke et al., 2015). Differences between GEMB and IMAU-FDM can be characterized as GEMB having lower FAC at low elevations and slightly higher FAC at higher elevations (Fig. 7). This pattern can be attributed to GEMB having warmer surface temperatures (higher melt) at lower elevations and a smaller b830 scaling coefficient that will impact total FAC most in areas with deep firn (i.e., the Antarctic Plateau).

It should be noted that for the runs presented here GEMB assumes an incompressible ice density of 910 kg m−3, while IMAU-FDM uses 917 kg m−3. The overall impact of this assumption, when comparing total FAC between models, will at worst be a 0.76 % underestimate by GEMB for those columns that do not densify to ice by the time they reach the bottom of the model domain (i.e., Antarctic interior) relative to IMAU-FDM. For all other locations the difference will be significantly less.

Looking at the temporal evolution in FAC for both Greenland (Fig. 8) and Antarctica (Fig. 9) it can be seen that seasonal fluctuations and interannual variations are in very close agreement for both ice sheets, but deviations in long-term trends are apparent. For the model runs examined GEMB is spun up using the 1979–1988 climatology repeated for 5000 years for Greenland and the 1979–2009 climatology repeated for 6000 years for Antarctica. IMAU-FDM is spun up using the 1960–1979 climatology repeated for as many years as it takes for the mean annual accumulation rate multiplied by the spin-up length to equal the thickness of the firn layer (from the surface to the depth at which the ice density is reached: Kuipers Munneke et al., 2015).

For Greenland (Fig. 8c) GEMB has virtually no trend in FAC between 1979 and 2005, becoming more negative thereafter coincident with increases in summer melt. IMAU-FDM FAC trends positive between 1979 and 2005 after which the trend becomes negative, closely matching the rate of FAC loss simulated by GEMB. Differences in FAC trend between 1979 and 2005 can be attributed to differences in model spin-up that are known to be a major source of uncertainty in FAC trends (Kuipers Munneke et al., 2015). Since GEMB is spun up using the 1979–1988 climatology, there should be little trend in FAC over this period. GEMB was initialized to the 1979–1988 climatology as the 1960–1979 period was not included with the provided RACMO forcing. In addition, climate reanalyses are known to perform considerably better after the introduction of satellite data in 1978. This is especially true over the poles where in situ observations are sparse (Tennant, 2004). The difference in prescribed spin-up climatology between the two models is the most likely cause of the observed differences in FAC trend between 1979 and 2005. When the 1979–2005 FAC trend is removed from both products (not shown), the FAC time series are nearly identical until 2004 after which GEMB estimates ∼0.5 m of FAC loss and IMAU-FDM estimates ∼1.0 m of FAC loss between 2004 and 2015. Some of the difference in trend can be attributed to IMAU-FDM having higher rates of melt along the equilibrium altitude during this period of time (Fig. 5).

Figure 10Rate of firn air content change (dFAC) for the Antarctic Ice Sheet as modeled by GEMB (a, d) and IMAU-FDM (b, e), as well as their difference (c, f). Rates for the period 1979–2008 are shown in panels (a)(c), and rates for the 2008–2015 period are shown in panels (d)(f).

Changes in FAC for the Antarctic (Fig. 9c) are nearly identical between models with a notable divergence beginning in 2008, after which GEMB FAC trends slightly positive and IMAU-FDM FAC trends slightly negative. Spatial patterns of FAC trends for the period pre- and post-2008 (Fig. 10) show much larger rates of FAC change for the 2008–2015 period with IMAU-FDM tending to have the same sign but larger trends than GEMB. For the post-2008 period, GEMB FAC gains outcompete losses, resulting in a slight positive trend in FAC over this period that is not seen in IMAU-FDM.

Despite a desire to identify which model is closer to the truth, we do not yet have an objective way to determine that. Likely the most definitive analysis will be to compare FAC-corrected altimetry (e.g., ICESat/2) estimates of ice sheet mass change to mass changes derived from satellite gravimetry (i.e., GRACE/FO) data, but that is beyond the scope of this study.

5 Model sensitivity to vertical resolution

Firn and surface mass balance models are highly sensitive to the chosen setup and empirical parameterizations (Bougamont et al., 2007; Kuipers Munneke et al., 2015). Major sources of uncertainty include parameterizations of albedo, snow grain growth, surface roughness, densification and its calibration (Stevens et al., 2020), and thermal conductivity. These uncertainties become exacerbated with the introduction of liquid water into the firn column (Vandecrux et al., 2020). Model setup can also introduce large sources of uncertainty and error. Three particularly important decisions are the choice of spin-up climatology (Kuipers Munneke et al., 2015), the length of time the spin-up is run (needs to reach equilibrium to prevent model drift), and the vertical resolution of the model. Most of these sensitivities have been explored extensively, and we direct the reader to previous publications that explore these topics in detail (e.g., Bougamont et al., 2007; Kuipers Munneke et al., 2015; Vandecrux et al., 2020; Stevens et al., 2020; Lundin et al., 2017).

One model sensitivity that we have not seen covered elsewhere is the sensitivity of the model to the vertical resolution of the firn column. The vertical resolution can have a large impact on melt rates and surface temperatures. This is because all of the energy and mass transfer between the atmosphere and the surface are often allocated to the model's uppermost layer in the firn column. A thicker layer will tend to dampen fluctuations in thermal energy relative to a thinner layer. Since snow grain growth, snow–ice melt, and compaction have nonlinear relations to temperature, model results will diverge for differing vertical resolutions. This is most true for the near-surface layers as thermal gradients are attenuated with depth.

Figure 11GEMB change in the modeled rate of melt (a–c) as well as evaporation and condensation (d–f) for near-surface model layer depths (dztop) of 5, 10, and 20 cm relative to a near-surface model layer depth of 2 cm.

To demonstrate model sensitivity to vertical resolution we run GEMB v1.0 for Greenland using the same model setup (Table 1) but for four different dztop [2, 5, 10, 20 cm] and dzmin [1, 2.5, 5, 10 cm] pairs (i.e., four model runs each having dzmin=1/2dztop). Results of the simulations are shown in Fig. 11. Increasing the size of the top layers results in a progressive reduction in melt as diurnal peak temperatures are muted and vertical gradients in temperature are reduced, in turn reducing grain metamorphism and surface darkening. Increasing the size of the top layers also results in stronger thermal gradients between the atmosphere and the surface that drive larger (both positive and negative) surface latent and sensible heat and mass fluxes. This has more impact in summer when diurnal fluctuations are larger than during polar night when diurnal changes in temperature are small or nonexistent. The net effect of increasing the surface layer thickness is slightly lower evaporation and sublimation losses.

6 Comparison of modeled and observed near-surface temperatures

To demonstrate the model's skill at simulating near-surface temperature, we compare GEMB v1.0 model output to an observational dataset of near-surface firn temperature that was collected at Summit Station, Greenland (Miller et al., 2017). Here, subsurface temperatures are measured using thermistors buried within the snowpack every 20 cm from July 2013 to June 2014. We compare these observations against results produced by GEMB v1.0 forced with RACMO output (see Table 1), rerun to output concurrent daily model solutions for the location of the observations (72.580 N, 38.459 W). We correct the data for change in thermistor depth over time. By default, GEMB's thermal conductivity is calculated as a function of snow density, after Sturm et al. (1997). To investigate model sensitivity to thermal conductivity, we also run the same simulation using the density-dependent thermal conductivity relation suggested by Calonne et al. (2011, Eq. 12). As demonstrated by Calonne et al. (2011), their solution gives a higher thermal conductivity for an equivalent density, which sits just inside the upper bounds of the 95 % confidence interval reported by Sturm et al. (1997). Here we compare results between the two GEMB simulations and the observations. For simplicity, all results are plotted with respect to the modeled or observed surface. Note that since the first thermal probe was placed at a depth of 20 cm below the surface, observational values between the snow surface and 20 cm depth are extrapolations.

Figure 12Observed and modeled daily near-surface (upper 2 m) temperatures for Summit Station, Greenland (72.580 N and 38.459 W), from July 2013 to June 2014. Temperature profiles from (a) observations (Miller et al., 2017), (c) modeled using GEMB and Sturm et al. (1997) thermal conductivity, and (e) modeled using GEMB and Calonne et al. (2011) thermal conductivity. (b) Observed downward shortwave (dsw) and longwave (dlw) radiation provided for context. Differences between temperature profiles for (d) GEMB–Sturm and observations and for (f) GEMB–Calonne and observations. Note: all GEMB simulations are forced with RACMO model output (see Table 1), and therefore modeled temperature differences relative to observations are due to errors in forcing, model parameterizations, and observations.

Figure 12a, c, and e show the thermal profiles to a depth of 2 m for observed, GEMB–Sturm, and GEMB–Calonne temperatures, respectively. Figure 12b shows the observed radiative fluxes, and Fig. 12d and f show the difference between modeled and observed snow temperature profiles. Overall, there is very good agreement between the modeled and observed thermal profiles with differences resulting from errors in atmospheric forcing, errors in observations, and biases attributed to errors in model parameterizations. Both model simulations produce mean temperatures between 0.2 and 2 m depth that are ∼0.8 K warmer than the observations, with the Sturm parameterization of thermal conductivity producing slightly better agreement (0.77 K) than the Calonne parameterization (0.79 K). Looking at bias and root mean square error as a function of depth, it can be seen (Fig. 13) that the Sturm parameterization produces slightly better agreement with the observations in the top 0.5 m of the snowpack, while the Calonne parameterization produces a significantly better fit for depths below 0.8 m. From this single location comparison, the Calonne parameterization outperforms the Sturm parameterization. Drawing any more definitive conclusions is challenging given other sources of disagreement such as atmospheric forcing and density structure.

Figure 13Mean bias (left) and RMSE (right) in modeled near-surface temperatures as a function of depth relative to observations (Miller et al., 2017). Daily GEMB model results using Sturm et al. (1997) and Calonne et al. (2011) thermal conductivity parameterizations for Summit Station, Greenland (72.580 N and 38.459 W), from July 2013 to June 2014.

While differences in mean temperature are only ∼0.8 K, seasonal differences can be as large as 5 K (Fig. 13). A 5 K increase in snow–firn temperature increases the compaction rate by a factor of 1.6–1.8. Since compaction rates are strongly dependent on temperature, seasonal biases in temperate have a large impact on rates of firn compaction, even if annual biases are small. Therefore, future work should prioritize improving near-surface atmospheric temperatures over glacier surfaces, from both models and observations, and thermal diffusion within snow and ice.

Calonne et al. (2019) show that the quadric relations between effective thermal conductivity and density, which are often used for snow (e.g., Eqs. 4 and 5), do not perform well for firn and porous ice. Instead, they propose an alternative empirical relation to more accurately model effective conductivity for the full range of densities from 0 to 917 kg m−3 (see Eq. 5 in Calonne et al., 2019). We plan to include empirical relations proposed by Calonne et al. (2019) in a future version of GEMB.

7 Summary and conclusions

This paper provides the first description of the open-source Glacier Energy and Mass Balance (GEMB) model that has been integrated as a module into the Ice-sheet and Sea-level System Model (ISSM). The model is one-way coupled with the atmosphere, which allows it to be run offline with a diversity of climate forcing but neglects feedback to the atmosphere. GEMB is written in C++, which produces efficient compiled machine code for fast execution. GEMB provides numerous parameterization choices for various key processes (e.g., albedo, subsurface shortwave absorption, and compaction), making it well suited for uncertainty quantification and model exploration.

To evaluate output from GEMB we compare it to the model output from the Institute for Marine and Atmospheric Research Utrecht Firn Densification Model (IMAU-FDM: Ligtenberg et al., 2011) for both Greenland (Ligtenberg et al., 2018) and the Antarctic (Ligtenberg et al., 2011). Models are forced with the same RACMO climate data (see the “Study-specific model setup” section, Table 1) and are independently calibrated to ice core density profiles. By forcing GEMB with the same climate data as IMAU-FDM we are able to attribute differences in model output to differences in model parameterizations and setup.

Overall, we find good agreement between models for both ice sheets with a few notable differences. For Greenland GEMB produces considerably less melt than IMAU-FDM but nearly as much runoff. This is because, relative to IMAU-FDM, GEMB tends to generate more low-elevation melt that runs off and less high-elevation melt that refreezes within the firn. These differences are most likely due to inter-model differences in the bare ice albedo and fresh snowmelt. GEMB and IMAU-FDM produce considerably different trends in FAC that are the result of differing climatologies used for model spin-up. GEMB used the 1979–1988 climatology repeated for 5000 years, while IMAU-FDM used the 1960–1979 climatology. This results in GEMB having little trend in FAC prior to the onset of increased melting in 2005. In contrast, IMAU-FDM has a positive FAC trend from 1979–2005 and a negative trend thereafter. These differences can have a large impact on volume-to-mass conversions used to generate glacier mass change estimates from satellite altimetry. For the Antarctic, GEMB and IMAU-IMAU-FDM are nearly identical with GEMB having slightly lower rates of sublimation and slightly higher rates of meltwater runoff that lead to a slightly more positive surface mass balance trend than found in IMAU-FDM. Changes in FAC for the Antarctic are nearly identical between models with a notable divergence beginning in 2008, after which GEMB FAC trends slightly positive and IMAU-FDM FAC trends slightly negative.

To demonstrate the impact of model setup we explore the impact of the model's vertical resolution on the results. We show that coarsening the model's vertical resolution decreases melt and increases surface latent and sensible heat fluxes due to stronger thermal gradients between the surface and the atmosphere. Lastly, we compare modeled near-surface thermal profiles, calculated using two different density–thermal conductivity relations, to in situ observations collected at Summit Station, Greenland. Our comparison shows good agreement between modeled and observed temperatures, with slightly improved agreement using the Calonne versus the Sturm thermal conductivity parameterization.

Our analysis shows that the GEMB model performs as expected and produces results that are comparable to an existing state-of-the-art firn model (IMAU-FDM). There are notable differences in output between the models, but it is difficult to judge if one model outperforms the other as there are insufficient observational data to distinguish between models. Future studies that compare altimetry-derived mass changes with those derived from satellite gravimetry should help to distinguish between models.

Future work will focus on developing atmospheric and firn downscaling routines that will allow GEMB to produce higher-resolution output when forced with medium-resolution (30–100 km) climate reanalysis or coarse-resolution (200–500 km) climate model output.

Appendix A

Table A1The following tables list the cores, taken from Montgomery et al. (2018) and Medley et al. (2022), and their citations. Those used for deep core calibration analysis (density of 830 kg m−3) and shallow core calibration analysis (density of 550 kg m−3) are marked in the corresponding “Deep” and “Shallow” columns. Cores without a deep or shallow designation were used for model evaluation.

Download XLSX

Code availability

GEMB is a module implemented inside of ISSM. The archived version of the source code used in this paper is made available as part of a Zenodo repository at (ISSM Team, 2022); code can be downloaded, compiled, and executed following the instructions available on the ISSM website: (last access: 26 August 2022). The public SVN repository for the ISSM code can also be found directly at (last access: 6 August 2022) and can be downloaded using username “anon” and password “anon”. The version of the code for this study, corresponding to ISSM release 4.21, is SVN version tag number 27238. An unmaintained MATLAB version of the model is also available (; Gardner, 2022). It is recommended to use the well-maintained ISSM version of GEMB.

Data availability

All GEMB and RACMO/IMAU-FDM model output used to generate Figs. 2 through 13 and Appendix A can be found at (Schlegel et al., 2022).

Author contributions

ASG wrote the model, co-wrote the paper, and made the figures. NJS ran all of the experiments and co-wrote the paper. EL translated the original model from MATLAB into C++ and integrated it into ISSM.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


We thank Brooke Medley for sharing the firn core data. We thank Brice Noël, Jan Melchior van Wessem, Stefan Ligtenberg, and Michiel van den Broeke for sharing the RACMO climate forcing and IMAU-FDM results. We thank Nathaniel Miller for directing us to the Greenland temperature profile observations, and we thank Martin Sharp for helping motivate early work on the model. We thank Charlotte Lang for helping with the migration and testing of GEMB in ISSM. The authors were supported by the NASA Cryospheric Sciences, ICESat-2, MEaSUREs, and the NASA Sea Level Change Team programs. Funding was also provided by the JPL Research and Technology Development program. We gratefully acknowledge computational resources and support from the NASA Advanced Supercomputing Division. This research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration and funded through the internal Research and Technology Development program.

Financial support

This research has been supported by the National Aeronautics and Space Administration, NASA Headquarters.

Review statement

This paper was edited by Andrew Wickert and reviewed by C. Max Stevens and one anonymous referee.


Alexander, P. M., Tedesco, M., Fettweis, X., van de Wal, R. S. W., Smeets, C. J. P. P., and van den Broeke, M. R.: Assessing spatio-temporal variability and trends in modelled and measured Greenland Ice Sheet albedo (2000–2013), The Cryosphere, 8, 2293–2312,, 2014. 

Alley, R. B.: Firn densification by grain-boundary sliding: a first model, Le Journal de Physique Colloques, 48, C1–249, 1987. 

Arthern, R. J. and Wingham, D. J.: The Natural Fluctuations of Firn Densification and Their Effect on the Geodetic Determination of Ice Sheet Mass Balance, Clim. Change, 40, 605–624,, 1998. 

Arthern, R. J., Vaughan, D. G., Rankin, A. M., Mulvaney, R., and Thomas, E. R.: In situ measurements of Antarctic snow compaction compared with predictions of models, J. Geophys. Res., 115, F03011,, 2010. 

Baker, I.: NEEM Firn Core 2009S2 Density and Permeability, Arctic Data Center [data set],, 2016. 

Bartelt, P. and Lehning, M.: A physical SNOWPACK model for the Swiss avalanche warning Part I: Numerical model, Cold Reg. Sci. Technol., 35, 123–145, 2002. 

Bassford, R. P.: Geophysical and numerical modelling investigations of the ice caps in Severnaya Zemlya, PhD thesis, University of Bristol, England, 2002. 

Beljaars, A. C. M. and Holtslag, A. A. M.: Flux parameterization over land surfaces for atmospheric models, J. Appl. Meteorol., 30, 327–341, 1991. 

Benson, C.: Greenland Snow Pit and Core Stratigraphy (Analog and Digital Formats), National Snow and Ice Data Center, Boulder, Colorado USA [data set],, 2013. 

Benson, C.: Greenland Snow Pit and Core Stratigraphy, Carl S. Benson Collection, Coll. 2010011, Roger G. Barry Archives and Resource Center, National Snow Data Center, 2017. 

Bolton, D.: The computation of equivalent potential temperature, Mon. Weather Rev., 108, 1046–1953, 1980. 

Bolzan, J. F. and Strobel, M.: Oxygen isotope data from snowpit at GISP2 Site 15, PANGAEA [data set],, 1999a. 

Bolzan, J. F. and Strobel, M.: Oxygen isotope data from snowpit at GISP2 Site 73, PANGAEA [data set],, 1999b. 

Bolzan, J. F. and Strobel, M.: Oxygen isotope data from snowpit at GISP2 Site 37, PANGAEA [data set],, 1999c. 

Bolzan, J. F. and Strobel, M.: Oxygen isotope data from snowpit at GISP2 Site 31, PANGAEA [data set],, 1999d. 

Bolzan, J. F. and Strobel, M.: Oxygen isotope data from snowpit at GISP2 Site 13, PANGAEA [data set],, 1999e. 

Bolzan, J. F. and Strobel, M.: Oxygen isotope data from snowpit at GISP2 Site 51, PANGAEA [data set],, 1999f. 

Bolzan, J. F. and Strobel, M.: Oxygen isotope data from snowpit at GISP2 Site 571, PANGAEA [data set],, 2001. 

Bougamont, M. and Bamber, J. L.: A surface mass balance model for the Greenland Ice Sheet, J. Geophys. Res.-Earth Surf., 110, F04018,, 2005. 

Bougamont, M., Bamber, J. L., Ridley, J. K., Gladstone, R. M., Greuell, W., Hanna, E., Payne, A. J., and Rutt, I.: Impact of model physics on estimating the surface mass balance of the Greenland Ice Sheet, Geophys. Res. Lett., 34, L17501,, 2007. 

Brils, M., Kuipers Munneke, P., van de Berg, W. J., and van den Broeke, M.: Improved representation of the contemporary Greenland ice sheet firn layer by IMAU-FDM v1.2G, Geosci. Model Dev., 15, 7121–7138,, 2022. 

Brun, E.: Investigation on wet-snow metamorphism in respect of liquid-water content, Ann. Glaciol., 13, 22–26, 1989. 

Brun, E., Martin, E., Simon, V., Gendre, C., and Coleou, C.: An energy and mass model of snow cover suitable for operational avalanche forecasting, J. Glaciol., 35, 333–342, 1989. 

Brun, E., David, P., Sudul, M., and Brunot, G.: A numerical model to simulate snow-cover stratigraphy for operational avalanche forecasting, J. Glaciol., 38, 13–22, 1992. 

Calonne, N., Flin, F., Morin, S., Lesaffre, B., du Roscoat, S. R., and Geindreau, C.: Numerical and experimental investigations of the effective thermal conductivity of snow, Geophys. Res. Lett., 38, L23501,, 2011. 

Calonne, N., Milliancourt, L., Burr, A., Philip, A., Martin, C. L., Flin, F., and Geindreau, C.: Thermal Conductivity of Snow, Firn, and Porous Ice From 3-D Image-Based Computations, Geophys. Res. Lett., 46, 13079–13089, 2019. 

Colbeck, S. C.: Theory of metamorphism of wet snow, Cold Regions Research and Engineering Laboratory (U.S.) Engineer Research and Development Center (U.S.),, 1973. 

Colbeck, S. C.: The capillary effects on water percolation in homogeneous snow, J. Glaciol., 13, 85–97,, 1974. 

Colbeck, S. C.: An overview of seasonal snow metamorphism, Rev. Geophys. Space Phys., 20, 45–61, 1982. 

Coléou, C. and Lesaffre, B.: Irreducible water saturation in snow: experimental results in a cold laboratory, Ann. Glaciol., 26, 64–68,, 1998. 

Conway, H.: Roosevelt Island Ice Core Density and Beta Count Data, U.S. Antarctic Program (USAP) Data Center [data set],, 2003. 

Copernicus Climate Change Service (C3S): ERA5: Fifth generation of ECMWF atmospheric reanalyses of the global climate. Copernicus Climate Change Service Climate Data Store (CDS),, last access: 25 January 2019. 

Courant, R., Friedrichs, K., and Lewy, H.: Über die partiellen Differenzengleichungen der mathematischen Physik, Mathematische Annalen, 100, 32–74, 1928. 

Culberg, R., Schroeder, D. M., and Chu, W.: Extreme melt season ice layers reduce firn permeability across Greenland, Nat. Commun., 12, 1–9, 2021. 

Ding, M., Yang, D., van den Broeke, M., Allison, I., Xiao, C., Qin, D., and Huai, B.: The surface energy balance at Panda 1 Station, Princess Elizabeth Land: A typical katabatic wind region in East Antarctica, J. Geophys. Res.-Atmos., 125, e2019JD030378,, 2020. 

Dunmire, D., Banwell, A. F., Wever, N., Lenaerts, J. T. M., and Datta, R. T.: Contrasting regional variability of buried meltwater extent over 2 years across the Greenland Ice Sheet, The Cryosphere, 15, 2983–3005,, 2021. 

Fausto, R. S., Box, J. E., Vandecrux, B., van As, D., Steffen, K., MacFerrin, M., Machguth, H., Colgan, W., Koenig, L. S., McGrath, D., Charalampidis, C., and Braithwaite, R. J.: A snow density dataset for improving surface boundary conditions in Greenland ice sheet firn modelling, Front. Earth Sci., 6,, 2018. 

Flanner, M. G. and Zender, C. S.: Linking snowpack microphysics and albedo evolution, J. Geophys. Res.-Atmos., 111, D12208,, 2006. 

Foken, T.: Micrometeorology, Springer Berlin, Heidelberg,, 2008. 

Forster, R. R., Box, J. E., van den Broeke, M. R., Miège, C., Burgess, E. W., van Angelen, J. H., Lenaerts, J. T. M., Koenig, L. S., Paden, J., Lewis, C., Gogineni, S. P., Leuschen, C., and McConnell, J. R.: Extensive liquid meltwater storage in firn within the Greenland ice sheet, Nat. Geosci., 7, 95–98,, 2014. 

Gardner, A. S.: Glacier Energy and Mass Balance – MATLAB (v0.1), Zenodo [code],, 2022. 

Gardner, A. S. and Sharp, M. J.: A review of snow and ice albedo and the development of a new physically based broadband albedo parameterization, J. Geophys. Res., 115, F01009,, 2010. 

Greuell, W. and Konzelmann, T.: Numerical modelling of the energy balance and the englacial temperature of the Greenland Ice Sheet, Calculations for the ETH-Camp location (West Greenland, 1155 m a.s.l.), Global Planet. Change, 9, 91–114, 1994. 

Greuell, W. and Oerlemans, J.: The Evolution of the Englacial Temperature Distribution in the Superimposed Ice Zone of a Polar Ice Cap During a Summer Season BT – Glacier Fluctuations and Climatic Change, edited by: Oerlemans, J., 289–303, Springer Netherlands, 1989. 

Guyomarc'h, G. and Merindol, L.: Validation of an application for forecasting blowing snow, Ann. Glaciol., 26, 138–143, 1998. 

Helsen, M. M., van den Broeke, M. R., van de Wal, R. S. W., van de Berg, W. J., van Meijgaard, E., Davis, C. H., Li, Y., and Goodwin, I.: Elevation changes in Antarctica mainly determined by accumulation variability, Science, 320, 1626–1629,, 2008. 

Herron, M. and Langway, C.: Firn Densification: An Empirical Model, J. Glaciol., 25, 373–385,, 1980. 

Hirashima, H., Avanzi, F., and Yamaguchi, S.: Liquid water infiltration into a layered snowpack: evaluation of a 3-D water transport model with laboratory experiments, Hydrol. Earth Syst. Sci., 21, 5503–5515,, 2017. 

Högström, U.: Review of some basic characteristics of the atmospheric surface layer, Bound.-Lay. Meteorol., 78, 215–246, 1996. 

Horlings, A. N., Christianson, K., Holschuh, N., Stevens, C. M., and Waddington, E. D.: Effect of horizontal divergence on estimates of firn-air content, J. Glaciol., 67, 287–296,, 2021. 

ISSM Team: Ice-sheet and Sea-level System Model source code, v4.21 r27238 (4.21), Zenodo [code],, 2022. 

Janssens, I. and Huybrechts, P.: The treatment of meltwater retention in mass-balance parameterizations of the Greenland ice sheet, Ann. Glaciol., 31, 133–140,, 2000. 

Kaspers, K. A., van de Wal, R. S. W., van den Broeke, M. R., Schwander, J., van Lipzig, N. P. M., and Brenninkmeijer, C. A. M.: Model calculations of the age of firn air across the Antarctic continent, Atmos. Chem. Phys., 4, 1365–1380,, 2004. 

Keenan, E., Wever, N., Dattler, M., Lenaerts, J. T. M., Medley, B., Kuipers Munneke, P., and Reijmer, C.: Physics-based SNOWPACK model improves representation of near-surface Antarctic snow and firn density, The Cryosphere, 15, 1065–1085,, 2021. 

Koenig, L. S., Miège, C., Forster, R. R., and Brucker, L.: Initial in situ measurements of perennial meltwater storage in the Greenland firn aquifer: Measurements of Greenland Aquifer, Geophys. Res. Lett., 41, 81–85,, 2014. 

Kreutz, K.: Microparticle, Conductivity, and Density Measurements from the WAIS Divide Deep Ice Core, Antarctica, U.S. Antarctic Program (USAP) Data Center [data set],, 2011. 

Kuipers Munneke, P., Ligtenberg, S. R. M., Noël, B. P. Y., Howat, I. M., Box, J. E., Mosley-Thompson, E., McConnell, J. R., Steffen, K., Harper, J. T., Das, S. B., and van den Broeke, M. R.: Elevation change of the Greenland Ice Sheet due to surface mass balance and firn processes, 1960–2014, The Cryosphere, 9, 2009–2025,, 2015. 

Larour, E., Morlighem, M., Seroussi, H., Schiermeier, J., and Rignot, E.: Ice flow sensitivity to geothermal heat flux of Pine Island Glacier, Antarctica, J. Geophys. Res., 117, F04023,, 2012a. 

Larour, E., Schiermeier, J., Rignot, E., Seroussi, H., and Morlighem, M.: Sensitivity Analysis of Pine Island Glacier ice flow using ISSM and DAKOTA, J. Geophys. Res., 117, F02009,, 2012b. 

Larour, E., Seroussi, H., Morlighem, M., and Rignot, E.: Continental scale, high order, high spatial resolution, ice sheet modeling using the Ice Sheet System Model (ISSM), J. Geophys. Res.-Earth Surf., 117, F01022,, 2012c. 

Lefebre, F., Gallée, H., van Ypersele, J.-P., and Greuell, W.: Modeling of snow and ice melt at ETH Camp (West Greenland): A study of surface albedo, J. Geophys. Res., 108, 4231,, 2003. 

Lenaerts, J. T. M., van den Broeke, M. R., Déry, S. J., König-Langlo, G., Ettema, J., and Munneke, P. K.: Modelling snowdrift sublimation on an Antarctic ice shelf, The Cryosphere, 4, 179–190,, 2010. 

Li, J. and Zwally, H.: Modeling the density variation in the shallow firn layer, Ann. Glaciol., 38, 309–313,, 2004. 

Ligtenberg, S. R. M., Helsen, M. M., and van den Broeke, M. R.: An improved semi-empirical model for the densification of Antarctic firn, The Cryosphere, 5, 809–819,, 2011. 

Ligtenberg, S. R. M., Kuipers Munneke, P., Noël, B. P. Y., and van den Broeke, M. R.: Brief communication: Improved simulation of the present-day Greenland firn layer (1960–2016), The Cryosphere, 12, 1643–1649,, 2018. 

Liston, G. E. and Elder, K.: A Distributed Snow-Evolution Modeling System (SnowModel), J. Hydrometeorol., 7, 1259–1276,, 2006. 

Lundin, J. M. D., Stevens, C. M., Arthern, R., Buizert, C., Orsi, A., Ligtenberg, S. R. M., Simonsen, S. B., Cummings, E., Essery, R., Leahy, W., Harris, P., Helsen, M. M., and Waddington, E. D.: Firn Model Intercomparison Experiment (FirnMICE), J. Glaciol., 63, 401–422,, 2017. 

MacFerrin, M., Machguth, H., van As, D., Charalampidis, C., Stevens, C. M., Heilig, A., Vandecrux, B., Langen, P. L., Mottram, R., Fettweis, X., van den Broeke, M. R., Pfeffer, W. T., Moussavi, M. S., and Abdalati, W.: Rapid expansion of Greenland's low-permeability ice slabs, Nature, 573, 403–407, 2019. 

Male, D. H. and Granger, R. J.: Snow surface energy exchange, Water Resour. Res., 17, 609–627, 1981. 

Marbouty, D.: An experimental study of temperature-gradient metamorphism, J. Glaciol., 26, 303–312, 1980. 

Marchenko, S., van Pelt, W. J. J., Claremar, B., Pohjola, V., Pettersson, R., Machguth, H., and Reijmer, C.: Parameterizing Deep Water Percolation Improves Subsurface Temperature Simulations by a Multilayer Firn Model, Front. Earth Sci., 5, 16,, 2017. 

Marsh, P. and Woo, M. K.: Wetting front advance and freezing of meltwater within a snow cover: 2. A simulation model, Water Resour. Res., 20, 1865–1874,, 1984. 

Mayewski, P. and Whitlow, S.: Snow Pit and Ice Core Data from Southern Greenland, 1984, Version 1.0, Arctic Data Center [data set],, 2009. 

Medley, B., Joughin, I., Das, S. B., Steig, E. J., Conway, H., Gogineni, S., Criscitiello, A. S., McConnell, J. R., Smith, B. E., van den Broeke, M. R., Lenaerts, J. T. M., Bromwich, D. H., and Nicolas, J. P.: Airborne-radar and ice-core observations of annual snow accumulation over Thwaites Glacier, West Antarctica confirm the spatiotemporal variability of global and regional atmospheric models, Geophys. Res. Lett., 40, 3649–3654,, 2013. 

Medley, B., Neumann, T. A., Zwally, H. J., Smith, B. E., and Stevens, C. M.: Simulations of firn processes over the Greenland and Antarctic ice sheets: 1980–2021, The Cryosphere, 16, 3971–4011,, 2022. 

Miège, C., Forster, R. R., Box, J. E., Burgess, E. W., McConnell, J. R., Pasteris, D. R., and Spikes, V. B.: Southeast Greenland high accumulation rates derived from firn cores and ground-penetrating radar, Ann. Glaciol., 54, 322–332,, 2013. 

Miller, H. and Schwager, M.: Density of ice core ngt37C95.2 from the North Greenland Traverse, PANGAEA [data set],, 2000a. 

Miller, H. and Schwager, M.: Density of ice core ngt42C95.2 from the North Greenland Traverse, PANGAEA [data set],, 2000b. 

Miller, J. Z., Culberg, R., Long, D. G., Shuman, C. A., Schroeder, D. M., and Brodzik, M. J.: An empirical algorithm to map perennial firn aquifers and ice slabs within the Greenland Ice Sheet using satellite L-band microwave radiometry, The Cryosphere, 16, 103–125,, 2022. 

Miller, N. B., Shupe, M. D., Cox, C. J., Noone, D., Persson, P. O. G., and Steffen, K.: Surface energy budget responses to radiative forcing at Summit, Greenland, The Cryosphere, 11, 497–516,, 2017. 

Montgomery, L., Koenig, L., and Alexander, P.: The SUMup dataset: compiled measurements of surface mass balance components over ice sheets and sea ice with analysis over Greenland, Earth Syst. Sci. Data, 10, 1959–1985,, 2018. 

Mosley-Thompson, E., McConnell, J. R., Bales, R. C., Li, Z., Lin, P.-N., Steffen, K., Thompson, L. G., Edwards, R., and Bathke, D.: Local to regional-scale variability of annual net accumulation on the Greenland ice sheet from PARCA cores, J. Geophys. Res., 106, 33839–33851,, 2001. 

Munro, D. S. and Davies, J. A.: An experimental study of the glacier boundary layer over melting ice, J. Glaciol., 18, 425–436,, 1977. 

Murphy, D. M. and Koop, T.: Review of the Vapour Pressures of Ice and Supercooled Water for Atmospheric Applications, Q. J. Roy. Meteor. Soc., 131, 1539–1565,, 2005. 

Murray, F. W.: On the Computation of Saturation Vapor Pressure. J. Appl. Meteorol., 6, 203–204, 1967. 

Noël, B., van de Berg, W. J., Machguth, H., Lhermitte, S., Howat, I., Fettweis, X., and van den Broeke, M. R.: A daily, 1 km resolution data set of downscaled Greenland ice sheet surface mass balance (1958–2015), The Cryosphere, 10, 2361–2377,, 2016. 

Parish, T. R. and Bromwich, D. H.: The surface windfield over the Antarctic ice sheets, Nature, 328, 51–54,, 1987. 

Patankar, S.: Numerical Heat Transfer and Fluid Flow, 1st edn., CRC Press,, 1980. 

Paterson, W.: The Physics of Glaciers, 3rd edn., Pergamon Press, Oxford, London, New York, 1994. 

Pfeffer, W. T., Illangasekare, T. H., and Meier, M. F.: Analysis and Modeling of Melt-Water Refreezing in Dry Snow, J. Glaciol., 36, 238–246,, 1990. 

Pfeffer, W. T., Meier, M. F., and Illangasekare, T. H.: Retention of Greenland Runoff by Refreezing – Implications for Projected Future Sea-Level Change, J. Geophys. Res.-Oceans, 96, 22117–22124, 1991. 

Schaaf, C. and Wang, Z.: MCD43C3 MODIS/Terra+Aqua BRDF/Albedo Albedo Daily L3 Global 0.05Deg CMG V006, NASA EOSDIS Land Processes DAAC [data set],, 2015. 

Schlegel, N.-J., Larour, E., Seroussi, H., Morlighem, M., and Box, J. E.: Decadal-scale sensitivity of northeast Greenland ice flow to errors in surface mass balance using ISSM, J. Geophys. Res.-Earth Surf., 118, 667–680,, 2013. 

Schlegel, N.-J., Larour, E., Seroussi, H., Morlighem, M. and Box, J. E.: Ice discharge uncertainties in Northeast Greenland from boundary conditions and climate forcing of an ice flow model, J. Geophys. Res., 120, 29–54,, 2015. 

Schlegel, N.-J. and Larour, E. Y.: Quantification of surface forcing requirements for a Greenland Ice Sheet model using uncertainty analyses, Geophys. Res. Lett., 46, 9700–9709,, 2019. 

Schlegel, N.-J., Gardner, A. S., and Larour, E.: Output from the Glacier Energy and Mass Balance (GEMB v1.0) forced with 3-hourly RACMO fields, Greenland and Antarctica 1979–2014 (1.1), Zenodo [data set],, 2022. 

Sjöblom, A.: Turbulent fluxes of momentum and heat over land in the High-Arctic summer: the influence of observation techniques, Polar Res., 33, 21567,, 2014. 

Smith, B., Fricker, H. A., Gardner, A. S., Medley, B., Nilsson, J., Paolo, F. S., Holschuh, N., Adusumilli, S., Brunt, K., Csatho, B., Harbeck, K., Markus, T., Neumann, T., Siegfried, M. R., and Zwally, H. J.: Pervasive ice sheet mass loss reflects competing ocean and atmosphere processes, Science, 12, 1239–1242, 2020. 

Steger, C. R., Reijmer, C. H., van den Broeke, M. R., Wever, N., Forster, R. R., Koenig, L. S., Kuipers Munneke, P., Lehning, M., Lhermitte, S., Ligtenberg, S. R. M., Miège, C., and Noël, B. P. Y.: Firn Meltwater Retention on the Greenland Ice Sheet: A Model Comparison, Front. Earth Sci., 5, 3,, 2017. 

Stevens, C. M., Verjans, V., Lundin, J. M. D., Kahle, E. C., Horlings, A. N., Horlings, B. I., and Waddington, E. D.: The Community Firn Model (CFM) v1.0, Geosci. Model Dev., 13, 4355–4377,, 2020. 

Sturm, M., Holmgren, J., König, M., and Morris, K.: The thermal conductivity of seasonal snow, J. Glaciol., 43, 26–41, 1997. 

Tennant, W: Considerations when using pre-1979 NCEP/NCAR reanalyses in the southern hemisphere, Geophys. Res. Lett., 31, L11112,, 2004. 

US International Trans-Antarctic Scientific Expedition (US ITASE): Glaciochemical Data, Version 2, edited by: Mayewski, P. A. and Dixon, D. A., US International Trans-Antarctic Scientific Expedition (US ITASE) Glaciochemical Data, Version 2, US_ITASE_Core Info-SWE-Density_2013.xlsx, National Snow and Ice Data Center, Boulder, Colorado USA, 2013. 

van Angelen, J. H., Lenaerts, J. T. M., van den Broeke, M. R., Fettweis, X., and van Meijgaard, E.: Rapid loss of firn pore space accelerates 21st century Greenland mass loss, Geophys. Res. Lett., 40, 2109–2113,, 2013. 

Vandecrux, B., MacFerrin, M., Machguth, H., Colgan, W. T., van As, D., Heilig, A., Stevens, C. M., Charalampidis, C., Fausto, R. S., Morris, E. M., Mosley-Thompson, E., Koenig, L., Montgomery, L. N., Miège, C., Simonsen, S. B., Ingeman-Nielsen, T., and Box, J. E.: Firn data compilation reveals widespread decrease of firn air content in western Greenland, The Cryosphere, 13, 845–859,, 2019. 

Vandecrux, B., Mottram, R., Langen, P. L., Fausto, R. S., Olesen, M., Stevens, C. M., Verjans, V., Leeson, A., Ligtenberg, S., Kuipers Munneke, P., Marchenko, S., van Pelt, W., Meyer, C. R., Simonsen, S. B., Heilig, A., Samimi, S., Marshall, S., Machguth, H., MacFerrin, M., Niwano, M., Miller, O., Voss, C. I., and Box, J. E.: The firn meltwater Retention Model Intercomparison Project (RetMIP): evaluation of nine firn models at four weather station sites on the Greenland ice sheet, The Cryosphere, 14, 3785–3810,, 2020. 

van den Broeke, M. R., Duynkerke, P. G., and Henneken, E. A. C.: Heat, momentum and moisture budgets of the katabatic layer over the melting zone of the West Greenland Ice Sheet in summer, Bound.-Lay. Meteorol., 71, 393–413, 1994. 

Van Der Veen, C. J.: Interpretation of short-term ice-sheet elevation changes inferred from satellite altimetry, Clim. Change, 23, 383–405,, 1993. 

van Wessem, J. M., van de Berg, W. J., Noël, B. P. Y., van Meijgaard, E., Amory, C., Birnbaum, G., Jakobs, C. L., Krüger, K., Lenaerts, J. T. M., Lhermitte, S., Ligtenberg, S. R. M., Medley, B., Reijmer, C. H., van Tricht, K., Trusel, L. D., van Ulft, L. H., Wouters, B., Wuite, J., and van den Broeke, M. R.: Modelling the climate and surface mass balance of polar ice sheets using RACMO2 – Part 2: Antarctica (1979–2016), The Cryosphere, 12, 1479–1498,, 2018. 

Veldhuijsen, S. B. M., van de Berg, W. J., Brils, M., Kuipers Munneke, P., and van den Broeke, M. R.: Characteristics of the contemporary Antarctic firn layer simulated with IMAU-FDM v1.2A (1979–2020), The Cryosphere Discuss. [preprint],, in review, 2022. 

Vionnet, V., Brun, E., Morin, S., Boone, A., Faroux, S., Le Moigne, P., Martin, E., and Willemet, J.-M.: The detailed snowpack scheme Crocus and its implementation in SURFEX v7.2, Geosci. Model Dev., 5, 773–791,, 2012. 

Warren, S. G. and Wiscombe, W. J.: A model for the spectral albedo of snow. II: Snow containing atmospheric aerosols, J. Atmos. Sci., 37, 2734–2745, 1980. 

Weinhart, A. H., Freitag, J., Hörhold, M., Kipfstuhl, S., and Eisen, O.: Representative surface snow density on the East Antarctic Plateau, The Cryosphere, 14, 3663–3685,, 2020.  

Wilhelms, F.: Density of ice core ngt03C93.2 from the North Greenland Traverse, PANGAEA [data set],, 2000a. 

Wilhelms, F.: Density of ice core ngt06C93.2 from the North Greenland Traverse, PANGAEA [data set],, 2000b. 

Wilhelms, F.: Density of ice core ngt14C93.2 from the North Greenland Traverse, PANGAEA [data set],, 2000c. 

Wilhelms, F.: Density of ice core ngt27C94.2 from the North Greenland Traverse, PANGAEA [data set],, 2000d. 

Wiscombe, W. J. and Warren, S. G.: A model for the spectral albedo of snow. I: Pure snow, J. Atmos. Sci., 37, 2712–2733, 1980. 

Zwally, H. J., Giovinetto, M. B., Li, J., Cornejo, H. G., Beckley, M. A., Brenner, A. C., Saba, J. L., and Yi, D.: Mass changes of the Greenland and Antarctic ice sheets and shelves and contributions to sea-level rise: 1992–2002, J. Glaciol., 51, 509–527,, 2005. 

Zwally, H. J. A. Y., Bindschadler, R. A., Brenner, A. C., Major, J. A., and Marsh, J. G.: Growth of Greenland Ice Sheet: Measurement, Science, 246, 1587–1589,, 1989. 

Short summary
This is the first description of the open-source Glacier Energy and Mass Balance (GEMB) model. GEMB models the ice sheet and glacier surface–atmospheric energy and mass exchange, as well as the firn state. The model is evaluated against the current state of the art and in situ observations and is shown to perform well.