Hostname: page-component-cd9895bd7-gxg78 Total loading time: 0 Render date: 2024-12-22T17:26:36.213Z Has data issue: false hasContentIssue false
  • English
  • Français

Surface mass balance and energy balance of the 79N Glacier (Nioghalvfjerdsfjorden, NE Greenland) modeled by linking COSIPY and Polar WRF

Published online by Cambridge University Press:  07 June 2021

M. T. Blau Open the ORCID record for M. T. Blau [Opens in a new window] ,
J. V. Turton Open the ORCID record for J. V. Turton [Opens in a new window] ,
T. Sauter Open the ORCID record for T. Sauter [Opens in a new window]  and
T. Mölg Open the ORCID record for T. Mölg [Opens in a new window]
M. T. Blau*
Affiliation:
Climate System Research Group, Institute for Geography, Friedrich-Alexander University, Erlangen, Germany Center for Climate Physics, Institute for Basic Science, Busan 46241, South Korea Department of Climate System, Pusan National University, Busan 46241, South Korea
J. V. Turton
Affiliation:
Climate System Research Group, Institute for Geography, Friedrich-Alexander University, Erlangen, Germany
T. Sauter
Affiliation:
Climate System Research Group, Institute for Geography, Friedrich-Alexander University, Erlangen, Germany
T. Mölg
Affiliation:
Climate System Research Group, Institute for Geography, Friedrich-Alexander University, Erlangen, Germany
*
Author for correspondence: M. T. Blau, E-mail: manuel.blau@pusan.ac.kr
Rights & Permissions [Opens in a new window]

Abstract

To get a better overview of atmosphere-driven mass changes at the 79N Glacier (Nioghalvfjerdsfjorden Glacier), the largest outlet glacier of the northeast Greenland ice stream, the surface mass balance (SMB) is modeled by linking the COupled Snowpack and Ice surface energy and mass-balance model in PYthon (COSIPY) with the output of a regional atmospheric model (Polar WRF) for the years 2014–2018. After a manual model optimization, the model produces reliable results when compared to observations in the region and to values from the literature. High spatial resolution (1 km) simulations reveal strong interannual variability of the SMB. Stronger surface melting increased the ablation and runoff in years with high mass loss (2016 and 2017) whereas in other years (2015 and 2018) melting and refreezing inside the snowpack dominated the mass balance (MB). A cooler regional climate with higher snowfall-driven accumulation, higher albedo and reduced surface melt in the ablation period of 2018 resulted in a positive SMB in 2018, however, the annual total MB remained negative. The results suggest a promising new dataset for gaining more insights into mass-balance processes and their contribution to the acceleration of glacier retreat in northeast Greenland.

Keywords

Glacier mass balanceglacier modelingice/atmosphere interactionsmass-balance reconstructionsurface mass budget
Type
Article
Information
Journal of Glaciology , Volume 67 , Issue 266 , December 2021 , pp. 1093 - 1107
Check for updates
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution, and reproduction in any medium, provided the original work is properly cited.
Copyright
Copyright © The Author(s), 2021. Published by Cambridge University Press

1. Introduction

The northeast Greenland ice stream (NEGIS) covers an area of ~200 000 km2 of the Greenland ice sheet (GrIS) and drains through three marine-terminating outlet glaciers (Fig. 1). Together, those drain ~12% of the GrIS (Fig. 1a; Mayer and others, Reference Mayer2018). The largest of those is called the Nioghalvfjerdsfjorden Glacier, also referred to in English (including throughout this paper) as 79N Glacier. Approximately 8.4% of the GrIS drains through the 79N Glacier to the ocean (Huybrechts and others, Reference Huybrechts, Mayer, Oerter and Jung-Rothenhäusler1999). With a drainage basin of ~120 000 km2, the 79N Glacier has the largest outlet of the GrIS and is the Arctic's largest remaining ice shelf (Seroussi and others, Reference Seroussi2011).

Fig. 1. Ice velocity of the GrIS from 14 December 2019 to 7 January 2020 (a) and the NEGIS from 6 October to 26 October 2016 (b) (Solgaard and Kusk, Reference Solgaard and Kusk2019). Site map of the 79N and the lower drainage basin of the NEGIS (c). The map in (c) shows the elevation taken from the PWRF output and presents the modeled domain. The isolines indicate 100 m (gray) and 1000 m (black) altitude steps. Stars show locations of the AWSs from PROMICE (Fausto and van As, Reference Fausto and van As2019) and points show the location of surface ablation stakes (Zeising and others, Reference Zeising, Neckel, Steinhage, Scheinert and Humbert2020).

The floating tongue of the 79N Glacier directly connects the GrIS to the Atlantic Ocean (Schaffer and others, Reference Schaffer2017). It is ~8 km wide at the grounding line and extends ~80 km long into the Nioghalvfjerds fjord. At the terminus, the glacier tongue widens to ~30 km (Thomsen and others, Reference Thomsen1997; Reeh and others, Reference Reeh, Mayer, Miller, Thomsen and Weidick1999; Reeh and others, Reference Reeh, Mayer, Olesen, Christensen and Thomsen2000). The glacier thickness increases from ~50 m at the glacier margins to 250–600 m at the grounding line (Huybrechts and others, Reference Huybrechts, Mayer, Oerter and Jung-Rothenhäusler1999; Reeh and others, Reference Reeh, Mayer, Olesen, Christensen and Thomsen2000; Schaffer and others, Reference Schaffer2016). The inland extent of the glacier is difficult to outline, as it joins with Zachariae Isstrøm (ZI) and Storstrømmen Glaciers to become the NEGIS. However, the glacier outlines can be approximately determined by the high ice velocity compared to the surrounding ice sheet (Fig. 1b; Krieger and others, Reference Krieger, Floricioiu and Neckel2020).

The 79N Glacier was assumed to be stable until the early 21st century (Joughin and others, Reference Joughin, Smith, Howat, Scambos and Moon2010). However, recent observations revealed increasing surface velocities and glacier thinning (Khan and others, Reference Khan2014). Dynamical mass loss through numerous small calving events led to an area reduction of ~65.4 km2 between 1999 and 2013 (Jensen and others, Reference Jensen, Box and Hvidberg2016). This is ~4% of the total area of the floating portion of the glacier. Andersen and others (Reference Andersen2019) reported an average annual area change of ~−2.7 km2 a−1 from 1999 to 2018. Additional retreat is now likely following the 100 km2 calving of Spaltegletsjer (a tributary of 79N) in 2019 and 2020 (personal communication from Jason Box, 4 September 2020). The glacier is exposed to both changes in atmosphere and ocean circulation. The trapping of relatively warm Atlantic water under the glacier tongue by a sill leads to a basal melt rate of 10.4 ± 3.1 m a–1 (Schaffer and others, Reference Schaffer2020). Enhanced basal melting in particular years is likely driven by atmospheric processes which drive additional warm water into the cavity of beneath the floating tongue (Münchow and others, Reference Münchow, Schaffer and Kanzow2020). Additional variability in the mass balance (MB) comes from the surface mass balance (SMB) changes (Mouginot and others, Reference Mouginot2019). The 79N Glacier region has experienced 3°C warming of the average air-temperature over the last 40 years (Turton and others, Reference Turton, Mölg and van As2019a), and persistent supraglacial melt lakes now feature on both the floating tongue and grounded part of the glacier (Leeson and others, Reference Leeson2015; Hochreuther and others, Reference Hochreuther, Neckel, Reimann, Humbert and Braun2021). The equilibrium line altitude (ELA) is likely to increase in altitude within the 21st century (Leeson and others, Reference Leeson2015; Noёl and others, Reference Noёl, van de Berg, Lhermitte and van den Broeke2019).

Due to the remoteness of the 79N Glacier there are only a couple of SMB measurements and climate observations from automatic weather stations (AWSs) available (Fig. 1). Therefore, to estimate SMB of the glacier, models are valid tools. Several studies have been conducted recently concerning SMB of Greenland. Those attempts use diverse models such as the regional atmospheric climate model RACMO2 (e.g. Ettema and others, Reference Ettema2009; Noёl and others, Reference Noёl2018; Mouginot and others, Reference Mouginot2019), the Modèle Atmosphérique Régional (MAR) (e.g. Fettweis, Reference Fettweis2007), a degree-day model (e.g. Mayer and others, Reference Mayer2018) or other climate-SMB models (e.g. The IMBIE Team, 2020) to detect mass changes, and flowline models to simulate ice velocity and glacier outlines (Krieger and others, Reference Krieger, Floricioiu and Neckel2020). Others achieved their results concerning the MB or ice thickness change with data from satellite missions or from other air-borne surveys (e.g. Seroussi and others, Reference Seroussi2011; Khan and others, Reference Khan2014). However, common focus was set on mass changes of the entire GrIS on a rather coarse grid and NEGIS was treated as a part of GrIS. This leads to an underrepresentation of the large spatial and temporal variability of SMB in the 79N region which cannot be captured in a coarse grid in its full entity (Mayer and others, Reference Mayer2018; Noёl and others, Reference Noёl2018; Turton and others, Reference Turton, Mölg and van As2019a).

Prior to the current study, high-resolution SMB modeling down to a spatial resolution of 2.5 km was only available until 1997 (Huybrechts and others, Reference Huybrechts, Mayer, Oerter and Jung-Rothenhäusler1999). Therefore, recent intensive melt events or the effects of warm air advection during winter (Turton and others, Reference Turton, Mölg and van As2019a), are not incorporated or poorly resolved by coarser-resolution SMB models which simulate the whole of Greenland. Other unsolved questions include the ratio of meltwater that runs off into the ocean, refreezes or is stored as liquid water in ice pores or meltwater ponds. Previously, Turton and others (Reference Turton, Mölg and Collier2020) applied the Polar version of the Weather Research and Forecasting Model (PWRF) with 1 km spatial resolution on the domain shown in Figure 1c, to investigate the dominant atmospheric processes and climatology of the region. The current study attempts to link the PWRF atmospheric output with a recently developed SMB model called the COupled Snowpack and Ice surface energy and mass-balance model in PYthon (COSIPY) (Sauter and others, Reference Sauter and Arndt2020). COSIPY was previously tested for Himalayan glaciers (Sauter and others, Reference Sauter and Arndt2020) with observations from AWS as atmospheric input. However, as yet, COSIPY has not been tested in a one-way, offline coupled setting, using high-resolution 4-D model output from an atmospheric model as input. Here, we use the PWRF data as input and explore which parameter settings for COSIPY are appropriate for our region of interest. The goals of this study are twofold. First, we evaluate the accuracy of COSIPY in its ability to simulate the SMB when linked with PWRF. Second, once we have confidence in the model for this region, we use it to produce a high-resolution (1 × 1 km) dataset of SMB for the 79N Glacier. This study will provide a basis for studying the acceleration of mass loss in the 79N region further and to investigate the relationship between specific atmospheric processes and SMB. The modeled data (available from www.doi.org/10.5281/zenodo.4434259) can provide the foundations of a deeper understanding of the processes effecting the MB of this remote region.

The paper continues with a description of the COSIPY model and the observations used to manually optimize and evaluate the model (Section 2). The results are split into different sections, beginning with model optimization and evaluation (Section 3.1) and followed by presentation and discussion of the surface mass and energy balance for the 79N Glacier from 2014 to 2018 (Section 3.2). The paper concludes in Section 4.

2. Data and methods

The concept of driving a process-based MB model with high-resolution atmospheric model output over multiple years was successful for mountain glaciers (e.g. Mölg and others, Reference Mölg, Maussion and Scherer2014), and hence it seems a promising approach for the current study. Furthermore, due to the sparse atmospheric observations in the region and the large area of 79N, using high-resolution input from an atmospheric model, as opposed to localized in situ observations, should yield more accurate results with regard to spatial variability. Several in situ observations of SMB, radiation and near-surface meteorological data are used for model optimization and evaluation of the results. Table 1 lists all data used in this study.

Table 1. List of input data from PWRF, AWS and stake observations including the time span of data usage and used variables

a 2 m air-temperature, 2 m air-pressure, incoming shortwave radiation, incoming longwave radiation, snowfall, total precipitation, wind-speed (horizontal wind vectors), initial information (skin temperature, water content, snow height, total height).

b Surface temperature, outgoing longwave radiation, shortwave albedo, snow height.

c Surface height change.

2.1. Study site

The regional climate is dominated by the characteristic annual cycle of the high-Arctic region, with the polar night (October to April), when no sunlight reaches the area, and the polar day (April to September). Two AWSs from the PROMICE network are located on the Kronprins Christian Land, just north of the 79N Glacier. A lower station is located in the ablation zone (KPC-L; where L stands for lower) and an upper station is near the equilibrium line (KPC-U; where U stands for upper; Fig. 1; Fausto and van As, Reference Fausto and van As2019). Cold conditions with annual mean temperatures (2009–18) of −13.4 ± 1.5°C (KPC-L) and −16.6 ± 0.5°C (KPC-U) characterize the local climate (Turton and others, Reference Turton, Mölg and van As2019a). The average summer temperature does not exceed the melting point at KPC-U, although daily maximum values can reach up to 4°C (Turton and others, Reference Turton, Mölg and van As2019a). Katabatic winds dominate the wind conditions. The strong temperature variability during winter is due to a combination of regional scale meteorological conditions, such as warm air advection from passing storms (Turton and others, Reference Turton, Mölg and van As2019a), and large-scale variability associated with teleconnections (Lim and others, Reference Lim2016; Hahn and others, Reference Hahn, Ummenhofer and Kwon2018).

2.2. Input data from Polar WRF

Data for the initialization of COSIPY comes from a high-resolution (1 km) Polar WRF (PWRF) output. PWRF is a polar-optimized version of the WRF model (Hines and others, Reference Hines and Bromwich2008). The atmospheric modeling was conducted and evaluated by Turton and others (Reference Turton, Mölg and Collier2020) using version v3.9.1.1 of PWRF. The output is available online (https://osf.io/53e6z/; last accessed on 9 May 2020). The data are available at hourly resolution for the domain shown in Figure 1c. Model outputs for this study from July 2014 to December 2018 are used. July 2014 is treated as a spin-up period, to ensure a minimum snow height (or no fresh snow as observed by a monthly average albedo of 0.47 at KPC-L, representing bare ice; Ryan and others, Reference Ryan2019) at the glacier surface starting August 2014 (Section 2.4). The experiment setup and parameterization for the PWRF modeling are described in Turton and others (Reference Turton, Mölg and Collier2020).

2.3. Observations

2.3.1. PROMICE weather stations

COSIPY is optimized with observations from two AWSs just north of 79N Glacier at two altitudes at the Kronprins Christian Land (Fig. 1c). KPC-L, in the ablation zone has an elevation of 255.6 m a.s.l., KPC-U is located near the equilibrium line at an elevation of 907.2 m a.s.l. (Table 1). The data are corrected (e.g. tilt-correction) and quality controlled according to the procedure by van As (Reference van As2011) and are available online at the PROMICE homepage (http://www.promice.org/PromiceDataPortal/; last accessed on 10 May 2020; Fausto and van As, Reference Fausto and van As2019). However, measurement errors cannot be completely excluded (Turton and others, Reference Turton, Mölg and van As2019a). Data from these locations were previously used to evaluate the PWRF atmospheric and radiation variables (Turton and others, Reference Turton, Mölg and Collier2020). To handle the noise in the snow height data, daily changes are calculated as difference in ‘midnight’ values defined as mean of all hourly values between 22:00 and 02:00 local time (LT), following Mölg and others (Reference Mölg2020). For COSIPY optimization the PROMICE data from January to December 2015 are used, as there are no missing values for all required input variables; the evaluation of COSIPY-WRF spans from August 2014 to December 2018 (Table 1).

2.3.2. Stake observations

In situ data from a field campaign conducted by the Alfred Wegener Institute (AWI) as a contribution to the Greenland Ice Sheet Ocean Interaction project (https://groce.de/en/; last accessed on 10 May 2020) are also used for independent evaluation of COSIPY-WRF. The observations are previously unpublished surface ablation estimates using ablation stakes on the 79N Glacier (Fig. 1c). Data are available for download from https://doi.pangaea.de/10.1594/PANGAEA.922131 (Zeising and others, Reference Zeising, Neckel, Steinhage, Scheinert and Humbert2020). The ablation stakes were installed in July 2017 and observations of snow height were taken in July 2018, which give the height change at the stake sites for 1 year. At each site, four stakes were installed, to prevent data loss in case of a failure of one stake, and to provide uncertainties in surface height change estimates.

The elevation and location of the stakes was recorded with a GPS device which has an elevation accuracy of ~10 m. The elevation of the stakes on the floating tongue is also subjected to an additional error of 3 m due to the tidal movement. We used the coordinates of the GPS to select the closest gridpoint in the COSIPY-WRF output for comparison. The uncertainties in stake location do not translate to an uncertainty of the surface height change, but to the elevation of the stake observations (personal communication from Zeising, 5 February 2020). The measured GPS locations are listed in Table 4 in Appendix A.

2.4. COSIPY

2.4.1. Structure and optimization

COSIPY is available online in its current release as an open-source software, we use v1.1 for this study (https://github.com/cryotools/cosipy; last accessed on 9 May 2020; Sauter and Arndt, Reference Sauter and Arndt2020). It is a physical surface energy and MB model that assumes full mass conservation in the snowpack. Therefore, it calculates the different mass fluxes by distinguishing between internal mass balance (IMB) and SMB. The sum of IMB and SMB is the total MB (TMB). In COSIPY, IMB summarizes the subsurface melt and the refreezing of meltwater, and SMB summarizes all MB processes at the surface including surface melt, snowfall, sublimation and deposition. Excess meltwater can be held in the snow layers as liquid water content. If one layer is saturated, meltwater drains in the next layer below (tipping bucket approach) or runs off when it reaches the surface below the lowest model layer (Sauter and others, Reference Sauter and Arndt2020). By doing so, the model divides the snow and ice pack into several user-defined numbers of layers in the vertical (Sauter and others, Reference Sauter and Arndt2020).

The ablation component relies on the surface energy balance (SEB). An energy surplus at the surface can cause melting (if the surface temperatures reach 0°C). A negative SEB, on the other hand, causes cooling of the surface layers. Sublimation or deposition also occurs depending on the surface-air water vapor pressure gradient between the surface and the surface layer. In COSIPY, SMB is computed as the difference of the total accumulation at the surface (snowfall accumulation and deposition) and the total ablation at the surface (surface melt and sublimation). Penetrating energy into the snowpack triggers ablation processes in the snowpack. Here, ablation occurs due to internal melting if temperatures reach 0°C. Melt water from surface and subsurface melt can penetrate deeper snow layers and refreeze there too or is held in pores as liquid water content. Excess meltwater contributes to the runoff and leaves the domain (Sauter and others, Reference Sauter and Arndt2020).

Carefully chosen parameter settings are needed to produce an ‘optimal’ model performance of COSIPY in this Arctic setup. For this reason, and as COSIPY is a newly released SMB model, a manual trial-and-error approach is preferred over an automated Monte-Carlo approach which is used for long-term documented MB models (e.g. Mölg and others, Reference Mölg, Maussion, Yang and Scherer2012; Cullen and others, Reference Cullen, Mölg, Conway and Steffen2014). A manual testing of parameter tuning gives a better idea of the response of COSIPY results to changes in the configuration. The testing followed the value ranges in Tables 5 and 6 in Appendix B and the final parameter selection is presented in Table 2. Input data from the two PROMICE AWSs are used from January to December 2015 for optimization (Section 2.3.1). In a second step, the model is applied to the 79N domain using the PWRF data as input.

Table 2. Parameterization and parameters selection for the COSIPY-WRF run

References are provided where value/value ranges have been used from previous studies. Each parameter was tested according to the options and ranges indicated by Tables 5 and 6 in Appendix B.

The initial englacial vertical profiles of the state variables temperature, liquid water content and density needed for the optimization are from the observations. The top-layer and bottom-layer temperatures are the mean temperatures of December 2014 and January 2015 measured at the surface and 10 m depth by the PROMICE AWSs (248.40 ± 6.4 K top, and 260.46 ± 0.6 K bottom for KPC-L; 244.87 ± 7.0 K top, and 258.66 ± 0.1 K bottom for KPC-U) (Table 2). Even though the 79N Glacier is reported to be much thicker than 10 m (e.g. Schaffer and others, Reference Schaffer2016), the assumption regarding the constant bottom-layer temperature is reasonable. The conductive heat flux below that depth was regarded to be negligible over the timescale examined following Cullen and others (Reference Cullen, Mölg, Conway and Steffen2014).

In COSIPY, the albedo is calculated with the parameterization of Oerlemans and Knap (Reference Oerlemans and Knap1998). The strong variability of ice albedo, which can range from 0.1 to 0.6 depending on the surface conditions, makes it difficult to find values for the entire 79N Glacier (Bøggild and others, Reference Bøggild, Brandt, Brown and Watten2010). Manual testing and inspection of the AWS data revealed that the albedo decreases from a fresh snow albedo of 0.87 to a firn albedo of 0.78 with an aging factor of 2 d after the snowfall event. When there is no snow layer, the albedo drops to 0.525 (van de Wal and Oerlemans, Reference van de Wal and Oerlemans1994). The depth scaling factor for albedo was found to be 12.5 cm (Table 2). The non-reflected shortwave radiation heats the surface, penetrates the snowpack and supplies the interior with energy according to Bintanja and van den Broeke (Reference Bintanja and van den Broeke1995). The outgoing longwave radiation is calculated after updating the surface temperature assuming a final surface emissivity of 0.99, following Noёl and others (Reference Noёl2018) (Table 2).

To calculate the turbulent heat fluxes, COSIPY assumes a linear roughness length increase with time following Mölg and others (Reference Mölg, Cullen, Hardy, Winkler and Kaser2009). The maximum increase from a fresh snow value of 0.5 mm to a firn value of 1.7 mm (Smeets and van den Broeke, Reference Smeets and van den Broeke2008; Fausto and others, Reference Fausto2016) occurs if there is no new snow addition over a specified time window, which has the same length as the albedo aging factor (Mölg and others, Reference Mölg, Cullen, Hardy, Winkler and Kaser2009). A value of 3.5 mm is assumed for bare ice (Cullen and others, Reference Cullen, Mölg, Conway and Steffen2014). The stability influence on the sensible and latent heat fluxes is based on the Richardson number (Table 2; Braithwaite, Reference Braithwaite1995).

Surplus energy results in ablation. If the surface temperature is below 273.16 K, only sublimation can occur. If the surface temperature exceeds the freezing point, it is set to a melting temperature of 273.16 K and the excess energy is translated into melt energy (Sauter and others, Reference Sauter and Arndt2020). Aging snow layers decrease their height due to compaction and densification. The density calculation with the ‘Boone’ method is described in Essery and others (Reference Essery, Morin, Lejeune and Ménard2013). Thereby, the density threshold of ice is 800 kg m−3 (Herron and Langway, Reference Herron and Langway1980). The threshold of firn is derived after manually testing parameter combinations following the tested range in Cullen and others (Reference Cullen, Mölg, Conway and Steffen2014) from 400 to 500 kg m−3. A minimum firn density of 400 kg m−3 was selected as it produces results closest to the observations (Table 2).

2.4.2. COSIPY-WRF configuration

The final COSIPY-WRF run was initiated on 15 July 2014. For the COSIPY-WRF simulations, the constant bottom-layer temperature was 258.7 K (the mean value for the observations at both AWSs) and was run with only the final set of parameters after the manual optimization. The first 15 d are regarded as a spin-up and are not included in the analysis. The spin-up includes the last 2 weeks of July where no snow accumulation was modeled or observed. The analyzed COSIPY-WRF output starts from 1 August 2014 (Table 1). At this time there is no fresh snow and minimal old snow observed and modeled at the glacier surface facilitating the initialization process of the model. The comparison against SEB terms and SMB from PROMICE AWSs is conducted for hourly values, except for albedo, for which a daily mean value is calculated. Albedo values are calculated as the sum of the observed outgoing shortwave radiation divided by the sum of the observed incoming shortwave radiation in 1 d. This approach is suggested by Oerlemans and Knap (Reference Oerlemans and Knap1998) to reduce errors from hours of low solar elevation.

3. Results and discussion

3.1. Model evaluation

3.1.1 Optimization period (January to December 2015)

In order to optimize the model set up, all time steps of year 2015 (1 January to 31 December) are evaluated according to the correlation coefficient (R) and the root mean square error (RMSE) to the PROMICE AWSs (Table 3). COSIPY-WRF simulates the observed conditions at both AWSs well in the optimization period. The model results agree with the measurements with relatively high significant correlation coefficients and a low RMSE (Table 3). For surface temperature, COSIPY-WRF archives a performance comparable to Huintjes and others (Reference Huintjes2015) (R = 0.95, RMSE = 3.44 K for KPC-L, and R = 0.94, RMSE = 2.03 K for KPC-U, compared with R = 0.98, RMSE = 2.12 K in Huintjes and others, Reference Huintjes2015). The evaluation statistics for albedo (R = 0.78, RMSE = 0.10 for KPC-L, and R = 0.62, RMSE = 0.12 for KPC-U) compare well with the statistics archived by Mölg and others (Reference Mölg, Cullen, Hardy, Kaser and Klok2008) (R = 0.69, RMSE = 0.10). This suggests a reasonable level of optimization.

Table 3. Evaluation statistics (correlation coefficient – R and root mean square error – RMSE) for COSIPY-WRF evaluation against PROMICE AWSs in the optimization and the evaluation period

The statistics for the evaluation period were calculated one time for the entire time series and a second time for the time series excluding the optimization period.

3.1.2 Model evaluation from August 2014 to December 2018

A longer-term model evaluation of COSIPY-WRF is conducted for the period from 1 August 2014 to 31 December 2018. The statistical comparisons are also included in Table 3; a visual comparison is presented in Figures 2 and 3 for KPC-L and KPC-U, respectively. The comparison against the observation from PROMICE AWSs reveals a good agreement of COSIPY-WRF with the observed surface temperature (R = 0.92, RMSE = 4.46 K for KPC-L, and R = 0.90, RMSE = 5.13 K for KPC-U) and outgoing longwave radiation (R = 0.92, RMSE = 16.73 W m−2 for KPC-L, and R = 0.91, RMSE = 18.78 W m−2 for KPC-U) (Figs 2a, b, 3a, b). For both variables, the simulated time series follow the annual cycle of the measurements and meet the values relatively accurately. The statistical evaluation for the surface temperature is similar to the evaluation in Sauter and others (Reference Sauter and Arndt2020). Statistics for the outgoing longwave radiation are similar to those considered as good in Noёl and others (Reference Noёl2018; R = 0.98 and RMSE = 14.20 W m−2 at their station S5). When excluding the optimization period from the evaluation period there is no systematic decrease in the model performance (Table 3).

Fig. 2. Model evaluation against the KPC-L of hourly surface temperature (a), hourly outgoing longwave radiation (b), daily shortwave albedo (c) and accumulated hourly surface height change relative to the first time step (d). The time period is 1 August 2014 00:00 to 31 December 2018 23:00 LT. The blue-dashed lines indicate the start and end time of the optimization period. Due to a measurement error of the SR50 sensor for the snow height change (stripping), the time before 1 July 2015 was excluded from evaluation for the snow height evaluation. Snow height can drop below zero due to a small initial snow layer at the observation site.

Fig. 3. Model evaluation against PROMICE the KPC-L AWS of hourly surface temperature (a), hourly outgoing longwave radiation (b), daily shortwave albedo (c) and accumulated hourly surface height change relative to the first time step (d). The time period is 1 August 2014 00:00 to 31 December 2018 23:00 LT. The blue-dashed lines indicate the start and end time of the optimization period. Snow height can drop below zero due to a small initial snow layer at the observation site.

The albedo is modeled with a relatively high correlation coefficient of R = 0.85 for KPC-L and R = 0.67 for KPC-U (Figs 2c, 3c; Table 3). At KPC-L, the decrease in albedo during summer, when ice is exposed, is delayed in the simulations (Fig. 2c). This is a result of finding the best albedo parameters valid for the entire 79N Glacier. To meet the KPC-L albedo values, a lower ice albedo and a higher scaling factor would produce a better model fit. However, a higher ice albedo and a lower albedo aging scaling factor produce a better agreement with the KPC-U albedo values. A similar difficulty in albedo parameterization choice for a large heterogeneous ice surface was identified by Noël and others (Reference Noёl2018). In previous studies, a better agreement in albedo values between observations and simulations were achieved by advancing the albedo parameterization of Oerlemans and Knap (Reference Oerlemans and Knap1998). For instance, Helsen and others (Reference Helsen2017) included both wet and dry snow albedos and various albedo aging factors. However, COSIPY does not currently have the option to choose wet vs dry albedos; therefore, we selected the most appropriate albedo values based on manual optimization of the model (Section 2.4.1 and Table 2).

In terms of snow height change, an offset is evident at both stations (Figs 2d, 3d). Especially at KPC-L, the modeled snow height is too high in winter 2015. This is partly attributed to errors in the measurements during the accumulation period. The albedo and snow height at KPC-U remain high, indicating snow cover, whereas at KPC-L, snow height is much lower despite the high albedo value. Note that the values in Figures 2d, 3d can be negative as it is a relative change since the first date of the simulation (1 August 2014) when a small layer of snow was present on the surface. The sonic rangers used to measure the distance between the surface and the sensor have a membrane which deteriorates over time and produces erroneous values (personal communication from van As D, 22 August 2019), which is a known potential problem with these instruments (Mölg and others, Reference Mölg2020). Therefore, a number of time steps from August 2014 until July 2015 were left out while calculating the evaluation statistics resulting in an RMSE of 0.39 m (Fig. 2d). COSIPY-WRF still overestimates the snow height observations at KPC-U but results in a smaller RMSE of 0.32 m compared to the RMSE at KPC-L (Fig. 3d). After summer 2015, the agreement of modeled and observed snow height improves for both PROMICE AWSs. A positive model bias in snow height can also result from a number of things, including an underestimation of the density by the parameterization. Due to the large area of 79N Glacier, finding a valid density parameter is difficult, as no density observations are available for this location. We tested both the ‘Bonne’ and ‘Vionnet’ parameterizations available in COSIPY and found higher agreement with observations for the Boone method. For the final run, the ‘Boone’ method was used as this was also reported to produce the best fit in previous studies (e.g. Essery and others, Reference Essery, Morin, Lejeune and Ménard2013). We chose this parameterization, despite having no observations to assess its validity for the 79N Glacier region. Issues in snow height simulation also rise from the effects of snowdrift that are not simulated by PWRF. Previous studies have shown a significant effect of snowdrift on the snow height in eastern Greenland (Hasholt and others, Reference Hasholt, Liston and Knudsen2003).

Figure 4 compares the modeled height change with the observations from the ablation stakes. Negative values point to ablation and decreasing surface height and positive values indicate accumulation and increasing surface height between July 2017 and July 2018. COSIPY-WRF overestimates ablation at some of the stake locations. At the glacier margin (five stakes; red points in Fig. 4), the modeled ablation is stronger compared to the measurements resulting in an RMSE of 0.97 m. At the stakes near the grounding line (blue points in Fig. 4), this offset decreases. Here, the COSIPY-WRF results agree better with the observations (RMSE = 0.32 m). Some ablation stakes at higher altitudes point to either ablation or accumulation (five stakes: black points in Fig. 4). COSIPY-WRF simulates ablation well at the stakes affected by a decreasing surface. However, the model results in a positive MB at only one of the three sites with observed accumulation. This points to a modeled 2017–18 ELA that is too high. Altogether, this leads to an RMSE of 0.67 m in the high-altitude group. The errors can be attributed to various reasons, such as measurement error related to tidal movement of the floating parts of the glacier tongue and GPS errors in elevation (see Section 2.3.2) (Reeh and others, Reference Reeh, Mayer, Olesen, Christensen and Thomsen2000; Joughin and others, Reference Joughin, Smith, Howat, Scambos and Moon2010; personal communication Zeising, 5 February 2020). The inaccuracy of the GPS sensors at the ablation stakes makes it difficult to find the exact gridpoint for the comparison with COSIPY-WRF (personal communication Zeising, 5 February 2020).

Fig. 4. Comparison of the modeled glacier height change with the GPS measurements at the ablation stakes from AWI. At each site, four stakes were deployed to prevent data loss in case of the failure of one stake and to provide error estimates. The error bars indicate the SD of the measured values at each site. Coloring shows reference to the location of the measurement sites like in Figure 1 (red: ice margin, blue: grounding line, black: high altitude region).

In spite of these significant measurement uncertainties, the simulated height change agrees relatively well with the observations. The model fit is highest at locations near the grounding line. The RMSE increases toward the lower areas of the floating tongue and toward the ELA. The floating tongue is highly crevassed and features many small melt ponds and channelized melt areas (Hochreuther and others, Reference Hochreuther, Neckel, Reimann, Humbert and Braun2021) making it difficult to find appropriate places for stake observations. The ELA in the 2017–18 ablation period was lower than that in previous years, as the accumulation during winter 2017 to 2018 was more intense than previous years (Turton and others, Reference Turton, Mölg and Collier2020; Section 3.3), which could be partially responsible for a poorer agreement between observations and simulations at the lower elevations. All in all, the evaluation against PROMICE AWSs and ablation stake measurements shows that COSIPY-WRF in its current configuration can offer a realistic estimate of SMB and SEB of the 79N Glacier and is therefore analyzed further in Section 3.2.

3.2. SMB and IMB

The modeled SMB is illustrated in Figure 5. The 79N Glacier features distinct patterns regarding the spatial distribution of annual SMB. The floating tongue is prone to negative SMB in all simulated years. There are variations in the estimates regarding the ELA and the extent of ablation in higher regions. Between September 2014 and August 2017, negative SMB regions exceed 1000 m a.s.l. In the 2014–15 period, SMB is negative up to an altitude of 1300 m a.s.l., however annual SMB total is small in magnitude accumulating to ~−13.56 × 106 kg a−1 (Fig. 5a). The mass loss intensifies in the 2015–16 and 2016–17 periods, with annual SMB values of −25.51 × 106 and −29.68 × 106 kg a−1, respectively (Figs 5b, c). Lower annual ablation in 2017–18 at the glacier tongue is depicted in Figure 5d. In the 2017–18 period, annual accumulation at the surface is evident already at 500 m a.s.l. In this year, SMB accumulates to ~+5.75 × 106 kg a−1.

Fig. 5. Spatial distribution of the annual SMB from September to August of the following year starting 2014 (a), 2015 (b), 2016 (c) and 2017 (d). The contours indicate 100 m (gray) and 1000 m (black) altitude steps taken from WRF.

Often, a negative TMB results in enhanced meltwater production. The surface water collects into meltwater ponds or supraglacial lakes, as seen on the surface of 79N (Leeson and others, Reference Leeson2015; Ignéczi and others, Reference Ignéczi2016; Hochreuther and others, Reference Hochreuther, Neckel, Reimann, Humbert and Braun2021). From there, a portion of meltwater evaporates, refreezes or runs off. Meltwater that penetrates the snowpack can either refreeze englacially, be collected in subglacial lakes or runs off at the base. The latter can lubricate the bedrock and can accelerate the ice flow inducing instability to other parts of the glacier (Carr and others, Reference Carr, Stokes and Vieli2013; Rathmann and others, Reference Rathmann2017; Bevis and others, Reference Bevis2019). Figure 6 illustrate the portion of surface meltwater that refreezes, runs off, evaporates or is stored in the snowpack (Figs 6a, b, c, d). Furthermore, the contribution of SMB and IMB to TMB is illustrated. Figures 6e, f, g h depict the simulated cumulative annual IMB – refreeze minus subsurface melt – and SMB – snowfall and deposition minus surface melt, sublimation and evaporation – contributing to the TMB of the 79N Glacier. Surface melt is a dominant factor determining ablation at the surface (Section 3.3).

Fig. 6. The ratio of surface meltwater that contributes to runoff, refreeze, evaporation or remains in the snow layer as liquid water content from September to August of the following year starting 2014 (a), 2015 (b), 2016 (c) and 2017 (d). The bottom panels show the contribution of SMB and IMB to the TMB of the same periods starting 2014 (e), 2015 (f), 2016 (g) and 2017 (h).

The major portion of the surface melt contributes to the runoff in all years, whereas evaporation is suppressed because of the low temperatures in the region (Fig. 9b in Appendix C). As a result, only 1% of the surface meltwater evaporates. Furthermore, a considerable part of meltwater is stored inside the snow layers. Schröder and others (Reference Schröder, Neckel, Zindler and Humbert2020) found evidence for liquid water retention under ice lenses and snow throughout the winter in our study area, which confirms that not all meltwater refreezes. Refreezing is more dominant in 2014–15 and 2017–18. This may be related to cooler climate conditions and low surface temperatures, which promote fast refreezing of meltwater before it runs off (Figs 9a, b in Appendix C).

Although during the 2014–15 and the 2017–18 periods (Figs 6e, h), IMB is the dominant factor for the TMB, the SMB does dominate over IMB in the periods from 2015–16 to 2016–17 (Figs 6f, g). The modeled TMB remains negative throughout the simulated time. In the years in which IMB is more negative than SMB, less than two-thirds of the meltwater contribute to the surface runoff (Figs 6a, d). In contrast, >85% of the total meltwater runs off in years when SMB contributes more to the TMB (Figs 6b, c). The only year with a positive component of the MB is 2017–18, but estimated positive SMB is compensated by negative IMB, resulting in an overall negative TMB.

SMB has a higher annual variability compared to IMB as it is more dependent on regional atmospheric conditions and its variability. Changes in meltwater runoff induce variabilities to the stability, the velocity and the total mass loss of the glacier (Fettweis, Reference Fettweis2007; Mayer and others, Reference Mayer2018; Bevis and others, Reference Bevis2019). The IMBIE Team (2020) supported that the meltwater production and runoff for the GrIS increased because of regional warming. Therefore, glacier-wide surface ablation can explain strong run off between 2015 and 2017. Furthermore, an overall negative MB in all periods agrees with Mouginot and others (Reference Mouginot2019). They identified an overall negative MB for the GrIS, which was largely due to internal and dynamic processes that did not allow accumulation due to surface processes even in years with positive SMB. This agrees with our findings for the 2017–18 period, when modeled positive SMB turns into a negative MB due to the IMB. Dynamical processes, such as iceberg calving or surging, are not included in the COSIPY model. However, during our study period, no major calving events were detected and 79N is not a surge-type glacier, unlike its southerly neighbor, Storstrømmen (Mouginot and others, Reference Mouginot2019).

3.3. Interannual variability of the SMB

The difference in SMB between the 2016–17 (Fig. 5c) and 2017–18 (Fig. 5d) years (which are the highest and lowest surface melt years in our study, respectively), prompted us to investigate the differences between the two years (hereafter ‘2017’ and ‘2018’). TMB and SMB exhibit distinct features. Ablation is strongly reduced in 2018 concerning both the total and the surface ablation (Figs 7b, d), whereas accumulation does not show large differences (Fig. 7a, c). Snowfall and surface melt during the ablation period are prominent components of SMB (Figs 7e, f, j), whereas deposition and sublimation are less dominant (Figs 7g, h, k, l). Due to the overall low surface and 2 m air-temperature, surface melting is suppressed and ablation due to sublimation gains importance in the accumulation period and in high altitude bands (Figs 7i, k; Figs 9a, b in Appendix C).

Fig. 7. Average vertical profile (altitude bands of 50 m) of the TMB (a and b) and SMB (c and d), and the components of the modeled SMB. This includes snowfall from PWRF (e and f), deposition (g and h), surface melt (i and j) and sublimation (k and l). The red line shows results for 2016–17 and the blue line shows results for 2017–18. The data are averaged over the accumulation periods (October 2016 to April 2017, and October 2017 to April 2018) in light blue boxes and over the ablation periods (May to September 2017 and May to September 2018) in brown boxes.

Distinct and significant differences are evident regarding the strength of surface melt and snowfall (paired t-test). Although in 2017 intensive melting appears in the ablation period, the melting is reduced in 2018 (Fig. 7j). Simultaneously, snowfall increased at elevations below 1200 m in the ablation period of 2018 compared to 2017 (Fig. 7f). However, there was no significant difference in snowfall amount during the accumulation periods of the two years (Fig. 7e). Only in lower altitudes of the ablation zone does the surface melt exceed mass gain by snowfall, while snowfall accumulation is persistent, leading to a distinct and significant reduction in negative SMB values in 2018 compared to 2017 (paired t-test).

SMB features are also expressed by the SEB (Fig. 8). Energy surplus during the ablation period is caused by rising incoming shortwave radiation and decreasing albedo resulting in a high positive net shortwave radiation (Figs 8b, d, and Fig. 10b in Appendix C), whereas in the accumulation period the shortwave radiation budget is near zero and albedo values are high (Figs 8a, c, and Fig. 10a in Appendix C; Fausto and others, Reference Fausto2016; Hahn and others, Reference Hahn, Storelvmo, Hofer, Parfitt and Ummenhofer2020). In the 2017 accumulation period, the longwave radiation is more negative at lower altitudes than in 2018. In the 2017 ablation period, a more negative longwave radiation is simulated at higher elevations. Overall, this leads to a weaker longwave radiation deficit in 2017 compared to 2018 (Figs 8e, f, and Fig. 10d in the Appendix C; Turton and others, Reference Turton, Mölg and van As2019a). Sensible heat adds energy to the surface throughout the year (Figs 8g, h). Surface temperatures at the melting point initiate the conversion of excess energy to melt energy during the ablation period (Figs 8l, Fig. 9a in Appendix C). The low longwave radiation flux coupled with a negative latent heat flux and a lack of energy input by shortwave radiation and by sensible heat suppresses melt energy, which approaches zero in all altitude bands in the accumulation period 2018 (Fig. 8k).

Fig. 8. Average vertical profile (altitude bands of 50 m) for the components of the SEB. This includes net shortwave radiation, SWnet (a and b), albedo (c and d), net longwave radiation, LWnet (e and f), sensible heat flux (g and h), latent heat flux (i and j) and melt energy (k and l) from COSIPY-WRF. The red line shows results for 2016–17 and the blue line shows results for 2017–18. The data are averaged over the accumulation periods (October 2016 to April 2017, and October 2017 to April 2018) in light blue boxes and over the ablation periods (May to September 2017 and May to September 2018) in brown boxes.

In the ablation period of 2017, the snow layer melted entirely, including at higher altitudes, resulting in the exposure of large extents of ice surface. As a consequence, the albedo was reduced and more shortwave radiation was absorbed and therefore available for ablation processes in summer 2017; a causality chain known to induce strong variability in MB also in other climatic zones (Mölg and others, Reference Mölg, Cullen, Hardy, Winkler and Kaser2009). This albedo feedback also explains positive SMB in the 2017–18 period. Higher snowfall in lower altitudes increased the snow height, which caused higher albedo and a decreased shortwave radiation budget at the surface. A lower snowline resulted in higher albedo, which reduced the melt energy and dampened surface melting in summer 2018. This snowline-albedo feedback at the margins of the GrIS was described for the period 2001–2017 by Ryan and others (Reference Ryan2019).

The COSIPY-WRF simulation points to a significantly warmer ablation period, especially in July, in terms of surface and air temperature in 2017 compared to 2018 in all altitude bands (paired t-test; Figs 9a, b in Appendix C). Additionally, wind speed and relative humidity from PWRF are significantly increased during both seasons in 2017 (paired t-test; Figs 9c, d in Appendix C). The differences in wind speed at the altitude bands near the grounding line are the only non-significant points (Fig. 9c in Appendix C). The strong and persistent wind speed in July 2017 promotes the intensification of turbulent energy fluxes in the ablation period (Figs 8h, j). Those increases in importance to SMB as other ablation processes (e.g. surface melt; Fig. 8i) are suppressed due to low temperatures. The respective mass turnover by sublimation and deposition has less of an impact on the TMB, however (Figs 7g, h, k, l). The snowfall (Fig. 7b) has a much larger impact on SMB, as seen over large parts of the glaciers in 2018 compared to 2017 (Figs 5c, d). In conclusion, less mass loss at the surface may be attributed to a significant increase in snowfall, and therefore albedo values, which was largely caused by a cooler regional climate in summer 2018 (The IMBIE Team, 2020).

4. Conclusion

The forcing of an SMB model with a high-resolution regional atmospheric model, PWRF, reproduced available observations reasonably well and enables more insights into mass changes in the remote area of the 79N Glacier. The manual model optimization conducted to meet the climatic observations from two PROMICE AWSs, however, also revealed problems in the parameterization for a large area. Often, input to SMB models, especially on smaller glaciers, comes from one or a few in situ observation sites (Mölg and others, Reference Mölg, Cullen, Hardy, Kaser and Klok2008; Sauter and others, Reference Sauter and Arndt2020). In this case, we attempted to find the best parameter selection for two stations, one in the ablation zone and one near the ELA, with quite different climatic characteristics. In addition, the simulation of snowfall and the related accumulation was a critical factor (Section 3.1). However, the modeled data from COSIPY-WRF were found to be suitable for a first detailed analysis of MB processes as response to atmospheric features.

SMB observations in this remote region are sparse and often plagued with difficulties. Independent stake data were only available for 1 year and were affected by tidal variations and GPS inaccuracies. However, our estimates of ~2 m mass loss on the floating tongue were also found in similar magnitude for NEGIS and 79N by Khan and others (Reference Khan2014). At the grounding line, the model results are more comparable with observations. Partially, this could be due to reduced tidal influence, but may also be due to parameter selections based on the KPC-L station, which is close to this location (Fig. 1).

The present results indicate stronger variability of ablation processes compared to accumulation processes on interannual timescale in the study domain. SMB features are dominated by ablation over the 79N Glacier in 2014–17 and turn positive in 2018. A major reason for this change in sign of the SMB is a strong variability in surface melt rates. In intense ablation years (September 2015 to August 2017) the mass loss due to surface ablation dominates the impact of IMB, and a major part of the meltwater contributes to the surface runoff. This finding is reversed in years with less overall ablation (September 2014 to August 2015) or positive SMB (September 2017 to August 2018). In those years, dominant factors inducing mass loss are components of IMB while SMB is small or positive. Also, most of the meltwater does not contribute to the surface runoff in those years. The TMB in the analyzed period remains negative (Fig. 6).

The COSIPY-WRF results reveal a substantial interannual variability in surface ablation, such as the reduced values in 2018, which can be explained by the effects of a cooler regional climate. The main differences between 2018 and the previous years with more intense overall ablation are increased snowfall, decreased surface melt and consequently higher albedo and lower net shortwave radiation budget in summer 2018. This difference effects large parts of the higher regions of the 79N Glacier, while being insignificant on the glacier tongue. Therefore, the climate in the lower ablation zone is potentially less influential over interannual SMB variability.

The low overall temperatures in this region suppress the ablation processes such as surface melting. Especially in the accumulation period, when the mean temperatures do not exceed the melting point, sublimation and deposition rise in importance, but remain in a low quantity. Melting at the NEGIS and 79N including the floating tongue occurs almost exclusively in June to August when net shortwave radiation is positive. In conclusion, the presented dataset reproduces SMB and SEB at the 79N Glacier reasonably well and highlights the large spatial and interannual variability of SMB and the individual components. The high resolution of SMB output is beneficial for future investigation into interactions between specific atmospheric features or extreme events, and the mass and energy balance of the 79N Glacier.

Acknowledgements

We thank Ole Zeising from the Alfred-Wegener Institute (AWI) for providing the data from the ablation stakes at the 79N Glacier and for additional information and communication throughout the project. Thanks are also to Prof. Jason Box for his information on the calving of 79N. This study was supported by the German Federal Ministry for Education and Research (BMBF) with grant number 03F0778F. We acknowledge the High-Performance Computing Centre (HPC) at the University of Erlangen-Nürnberg's Regional Computation Centre (RRZE), for their support and resources while running the COSIPY simulations. Daily average values of the SMB components from July 2014 to December 2018 are available online and open access at: www.doi.org/10.5281/zenodo.4434259. We express our thanks to the two anonymous reviewers and the editor for providing reviews of the paper, which benefitted from their critiques.

Author contributions

M.T.B. conducted optimization, simulation and evaluation of the COSIPY-WRF model and wrote the manuscript. J.V.T. designed the study. T.S. provided assistance in the setup and running of the model. T.M. assisted with model optimization and evaluation. J.V.T, T.S. and T.M. contributed to writing the manuscript and discussion of the results.

Conflict of interest

The authors confirm that they have no conflict of interest.

Appendix A: Ablation stake locations

See Table 4.

Table 4. GPS locations of the ablations stakes and the measurement results

Appendix B: Manual optimization scheme

See Tables 5 and 6.

Table 5. List of parameterization options available in v1.3 COSIPY (Sauter and others, Reference Sauter, Arndt and Schneider2020)

Table 6. List of key parameters and the range of values tested in model optimization

Appendix C: Local climate and radiative conditions

See Figures 9 and 10.

Fig. 9. Average vertical profiles (altitude bands of 50 m) for the surface temperature from COSIPY-WRF (a) and the mean climate conditions – 2 m air-temperature (b), 2 m wind speed (c) and relative humidity – from the PWRF output at the 79N region in July 2017 and July 2018. The surface temperature.

Fig. 10. Average vertical profile (altitude bands of 50 m) for the components of incoming shortwave radiation, SWin (a and b) and incoming longwave radiation, LWin (c and f) from COSIPY-WRF with energy flux density on the x-axis. The data are averaged over the accumulation periods (October 2016 to April 2017, and October 2017 to April 2018) and over the ablation periods (May to September 2017 and May to September 2018).

Footnotes

The altitude is measured with the GPS relative to an ellipsoid. Deviations can occur due to GPC inaccuracies (~10 m) and tidal movement of the glacier tongue (~3 m). ‘Difference’ is the difference of the snow height in July 2018 relative to July 2017 observations. At each site up to four measurements were takes. The exact number of measurements at each site are listed. The SD is the standard deviations of all measurements taken from one site (Zeising, and others 2020).

References

Andersen, JK and 17 others (2019) Update of annual calving front lines or 47 marine terminating outlet glaciers in Greenland (1999–2018). Geology of Greenland Survey Bulletin 43, e2019430202 doi: 10.34194/GEUSB- 201943-02-02.Google Scholar
Bevis, M and 13 others (2019) Accelerating changes in ice mass within Greenland, and the ice sheet's sensitivity to atmospheric forcing. Proceedings of the National Academy of Sciences of the United States of America 116(6), 1934–1939 doi: 10.1073/pnas.1806562116.CrossRefGoogle ScholarPubMed
Bintanja, R and van den Broeke, MR (1995) The surface energy balance of Antarctic snow and blue ice. Journal of Applied Meteorology 34(4), 902–924. doi: 10.1175/1520-0450(1995)034<0902:TSEBOA>2.0.CO;22.0.CO;2>CrossRefGoogle Scholar
Bøggild, CE, Brandt, RE, Brown, KJ and Watten, SG (2010) The ablation zone in northeast Greenland: ice types, albedos and impurities. Journal of Glaciology 56(195), 101–113. doi: 10.3189/002214310791190776CrossRefGoogle Scholar
Braithwaite, RJ (1995) Aerodynamic stability and turbulent sensible-eat flux over a melting ice surface, the Greenland ice sheet. Journal of Glaciology 41(139), 562–571. doi: 10.3189/S002214300003488CrossRefGoogle Scholar
Carr, RJ, Stokes, CR and Vieli, A (2013) Recent progress in understanding marine-terminating Arctic outlet glacier response to climatic and oceanic forcing: twenty years of rapid change. Progress in Physical Geography 37(4), 436–467. doi: 10.1177/0309133313483163CrossRefGoogle Scholar
Cullen, NJ, Mölg, T, Conway, J and Steffen, K (2014) Assessing the role of sublimation in the dry snow zone of the Greenland ice sheet in a warming world. Journal of Geophysical Research Atmosphere 119(11), 6563–6577. doi: 10.1002/2014JD021557CrossRefGoogle Scholar
Essery, R, Morin, S, Lejeune, Y and Ménard, CB (2013) A comparison of 1701 snow models using observations from an alpine site. Advances in Water Resources 55, 131–148. doi: 10.1016/j.advwatres.2012.07.013CrossRefGoogle Scholar
Ettema, J and 6 others (2009) Higher surface mass balance of the Greenland ice sheet revealed by high-resolution climate modeling. Geophysical Research Letters 36(12), D06116. doi: 10.1029/2009GL038110.CrossRefGoogle Scholar
Fausto, RS and 5 others (2016) The implication of nonradiative energy fluxes dominating Greenland Ice Sheet exceptional ablation area surface melt in 2012. Geophysical Research Letters 43(6), 2649–2658. doi: 10.1002/2016GL067720CrossRefGoogle Scholar
Fausto, RS and van As, D (2019) Programme for monitoring of the Greenland ice sheet (PROMICE): Automatic weather station data. Version: v03. Dataset published via Geological Survey of Denmark and Greenland doi: 10.22008/promice/data/aws.CrossRefGoogle Scholar
Fettweis, X (2007) Reconstruction of the 1979–2006 Greenland ice sheet surface mass balance using the regional climate model MAR. The Cryosphere 1(1), 21–40. doi: 1.5194/tc-1-21-2007CrossRefGoogle Scholar
Hahn, LC, Storelvmo, T, Hofer, S, Parfitt, R and Ummenhofer, CC (2020) Importance of orography for Greenland cloud and melt response to atmospheric blocking. Journal of Climate 33(10), 4187–4206. doi: 10.1175/JCLI- D-19-0527.1CrossRefGoogle Scholar
Hahn, L, Ummenhofer, CC and Kwon, YO (2018) North Atlantic natural variability modulates emergence of widespread Greenland melt in a warming climate. Geophysical Research Letters 45, 9171–9178. doi: 10.1029/2018GL079682CrossRefGoogle Scholar
Hasholt, B, Liston, GE and Knudsen, NT (2003) Snow-distribution modelling in the Ammassalik region, South East Greenland. Hydrological Research 34(1–2), 1–16. doi: 10.2166/nh.2003.0025CrossRefGoogle Scholar
Helsen, MM and 7 others (2017) On the importance of the albedo parameterization for the mass balance of the Greenland ice sheet in EC-earth. The Cryosphere 11(4), 1949–1965 doi: 10.5194/tc-11-1949-2017.CrossRefGoogle Scholar
Herron, MM and Langway, CC Jr (1980) Firn densification: an empirical model. Journal of Glaciology 25(93), 373–358 doi: 10.3189/S0022163000015239CrossRefGoogle Scholar
Hines, KM and Bromwich, DH (2008) Development and testing of polar weather research and forecasting (WRF) model. Part 1: Greenland Ice Sheet meteorology. Monthly Weather Review 136(6), 1971–1989. doi: 10.1175/2007MWR2112.1CrossRefGoogle Scholar
Hochreuther, P, Neckel, N, Reimann, N, Humbert, A and Braun, M (2021) Fully automated detection of supraglacial lake area for Northeast Greenland using sentinel-2 time-series. International Journal of Remote Sensing 13(2), 205. doi: 10.3390/rs13020205CrossRefGoogle Scholar
Huintjes, E and 9 others (2015) Evaluation of a coupled snow and energy balance model for Zhadang glacier, Tibetan Plateau, using glaciological measurements and time-lapse photography. Arctic Antarctic and Alpine Research 47(3), 573–590 doi: 10.1657/AAAR0014-073.CrossRefGoogle Scholar
Huybrechts, P, Mayer, C, Oerter, H and Jung-Rothenhäusler, F (1999) Climate change and sea level: Ice-dynamics and mass-balance studies on the Greenland ice sheet- EPIC3 Rep. on the Contribution of the Alfred Wegener Institute to EU Contract ENV4-CT95-0124, European Commission, DG XII, 18 pp. Available at https://epic.awi.de/id/eprint/3663/1/Huy1999b.pdf.Google Scholar
Ignéczi, A and 7 others (2016) Northeast sector of the Greenland Ice Sheet to undergo the greatest inland expansion of supraglacial lakes during the 21st century. Geophysical Research Letters 43(18), 9729–9738 doi: 10.1002/2016GL070338CrossRefGoogle Scholar
The IMBIE Team (2020) Mass balance of the Greenland Ice sheet from 1992 to 2018. Nature 579, 233–239. doi: 10.1038/s41586- 019-0855-2CrossRefGoogle Scholar
Jensen, TS, Box, JE and Hvidberg, CS (2016) A sensitivity study of annual area change for Greenland ice sheet marine terminating outlet glaciers: 1999–2013. Journal of Glaciology 62(231), 72–81. doi: 10.1017/jog.2016.12CrossRefGoogle Scholar
Joughin, I, Smith, BE, Howat, IM, Scambos, T and Moon, T (2010) Greenland Flow variability from ice-sheet-wide velocity mapping. Journal of Glaciology 56(197), 415–430. doi: 10.3189/00224310792447734CrossRefGoogle Scholar
Khan, SA and 12 others (2014) Sustained mass loss of the northeast Greenland ice sheet triggered by regional warming. Nature Climate Change 4, 292–299 doi: 10.1038/nclimate2161CrossRefGoogle Scholar
Krieger, L, Floricioiu, D and Neckel, N (2020) Drainage basin delineation for outlet glaciers of northeast Greenland based on sentinel-1 ice velocities and TanDEM-X elevations. Remote Sensing of the Environment 237, 111483. doi: 10.1016/J.rse.2019.111483CrossRefGoogle Scholar
Leeson, AA and 6 others (2015) Supraglacial lakes on the Greenland ice sheet advance inland under warming climate. Nature Climate Change 5, 51–55 doi: 10.1029/NCLIMATE2463.CrossRefGoogle Scholar
Lim, YK and 7 others (2016) Atmospheric summer teleconnections and Greenland Ice Sheet surface mass variations: insights from MERRA-2. Environmental Research Letters 11, 024002 doi: 10.1088/1748- 9326/11/2/024002.CrossRefGoogle Scholar
Mayer, C and 8 others (2018) Large ice loss variability at Nioghalvfjerdsfjorden Glacier, Northeast-Greenland. Nature Communication 9(2768), 1–11. doi: 10.138/s41467- 018-05180-x.CrossRefGoogle ScholarPubMed
Mölg, T and 8 others (2020) Mesoscale atmospheric circulation controls of local meteorological elevation gradients on Kersten Glacier near Kilimanjaro summit. Earth System Dynamics 11(3), 653–672 doi: 10.5194/esd- 11-653-2020.CrossRefGoogle Scholar
Mölg, T, Cullen, NJ, Hardy, DR, Kaser, G and Klok, L (2008) Mass balance of a slope glacier on Kilimanjaro and its sensitivity to climate. International Journal of Climatology 28(7), 881–892. doi:10.1002/joc.1589CrossRefGoogle Scholar
Mölg, T, Cullen, NJ, Hardy, DR, Winkler, M and Kaser, G (2009) Quantifying climate change in the Tropical Midtroposphere over East Africa from glacier shrinkage on Kilimanjaro. Journal of Climate 22(15), 4162–4181. doi: 10.1175/2009JCLI2954.1CrossRefGoogle Scholar
Mölg, T, Maussion, F and Scherer, D (2014) Mid-latitude westeries as a driver of glacier variability in monsoonal high Asia. Nature Climate Change 4, 68–73. doi: 10.1038/NCLIMATE2055CrossRefGoogle Scholar
Mölg, T, Maussion, F, Yang, W and Scherer, D (2012) The footprint of Asian monsoon dynamics in the mass and energy balance of a Tibetan glacier. The Cryosphere 6(6), 1445–1461. doi: 10.5194/tc- 6-1445-2012CrossRefGoogle Scholar
Mouginot, J and 8 others (2019) Forty-six years of Greenland Ice sheet mass balance from 1972 to 2018. Proceedings of the National Academy of Sciences of the United States of America 116(19), 9239–9244 doi: 10.1073/pnas.1904242116.CrossRefGoogle ScholarPubMed
Münchow, A, Schaffer, J and Kanzow, T (2020) Ocean circulation connecting fram strait to glaciers off Northeast Greenland: mean flows, topographic Rossby waves, and their forcing. Journal of Physical Oceanography 50(2), 509–530. doi: 1.1175/JPO-D-19-0085.1CrossRefGoogle Scholar
Noёl, B and 11 others (2018) Modelling the climate and surface mass balance of polar ice sheets using RACMO2 – Part 1: Greenland (1958–2016). The Cryosphere 12(3), 81–831 doi: 10.5194/tc- 12-811-2018.Google Scholar
Noёl, B, van de Berg, WJ, Lhermitte, S and van den Broeke, MR (2019) Rapid ablation zone expansion amplifies north Greenland mass loss. Science Advances 5(9), eaaw0123. doi: 10.1146/sciadv.aaw0123CrossRefGoogle Scholar
Oerlemans, J and Knap, WH (1998) A 1 year record of global radiation and albedo in the ablation zone of Morteratschgletscher, Switzerland. Journal of Glaciology 44(147), 231–238. doi: 10.3189/S0022143000002574CrossRefGoogle Scholar
Rathmann, NM and 7 others (2017) Highly temporally resolved response to seasonal surface melt of the Zachariae and 79N outlet glaciers in northeast Greenland. Geophysical Research Letters 44(19), 9805–9814 doi: 10.1002/2017GL074368.CrossRefGoogle Scholar
Reeh, N, Mayer, C, Miller, H, Thomsen, HH and Weidick, A (1999) Present and past climate control on fjord glaciations in Greenland: implications for IRD-deposition in the sea. Geophysical Research Letters 26(8), 1039–1042. doi: 10.1029/1999GL900065CrossRefGoogle Scholar
Reeh, N, Mayer, C, Olesen, OB, Christensen, EL and Thomsen, HH (2000) Tidal movement of Nioghalvfjerdsfjorden glacier, northeast Greenland: observations and modelling. Annals of Glaciology 31, 111–117.CrossRefGoogle Scholar
Ryan, J C and 6 others (2019). Greenland Ice sheet surface melt amplified by snowline migration and bare ice exposure. Science Advances 5(3), eaav3738 doi: 10.1126/sciadv.eaav3738.CrossRefGoogle ScholarPubMed
Sauter, T and Arndt, A (2020) Coupled snowpack and ice surface energy and mass balance model in Python (COSIPY) doi: 10.5281/zenodo.36132921 In URL: https://cryotools.github.io/cosipy/ (last accessed on April 14, 2020).CrossRefGoogle Scholar
Sauter, T, Arndt, A and Schneider, C (2020) COSIPY v1.3 – an open-source coupled snowpack and ice surface energy and mass balance model. Geoscientific Model Development 13(11), 5645–5662. doi: 10.5194/gmd- 13-5645-2020CrossRefGoogle Scholar
Schaffer, J and 6 others (2016) A global, high-resolution data set of ice sheet topography, cavity geometry, and ocean bathymetry. Earth System Science Data 8(2), 543–557 doi: 10.5194/essd- 8-543-2016.CrossRefGoogle Scholar
Schaffer, J and 6 others (2017) Warm water pathways toward Nioghalvfjerdsfjorden Glacier, Northeast Greenland. Journal of Geophysical Research Oceans 122(5), 4004–4020 doi: 10.1002/2016JC012462.CrossRefGoogle Scholar
Schaffer, J and 5 others (2020) Bathymetry constrains ocean heat supply to Greenland's largest glacier tongue. Nature Geoscience 13, 227–231. doi: 10.1038/s41561- 019-0529-xCrossRefGoogle Scholar
Schröder, L, Neckel, N, Zindler, R and Humbert, A (2020) Perennial supraglacial lakes in Northeast Greenland observed by Polarimetric SAR. Remote Sensing 12(17), 2798. doi: 10.3390/rs12172798CrossRefGoogle Scholar
Seroussi, H and 6 others (2011) Ice flux divergence anomalies on 79North glacier, Greenland. Geophysical Research Letters 38(9), L09501 doi: 10.1029/2011GL047338.CrossRefGoogle Scholar
Smeets, CJPP and van den Broeke, MR (2008) Temporal and spatial variation of the aerodynamic roughness length in the ablation zone of the Greenland Ice Sheet. Boundary Layer Meteorology 128, 315–338. doi: 1.1007/s10546-008-9291-0CrossRefGoogle Scholar
Solgaard, A and Kusk, A (2019) Programme for monitoring of the Greenland ice sheet (PROMICE): Greenland ice velocity. Geological survey of Denmark and Greenland (GEUS) doi: 10.22008/promice/data/sentinel1icevelocity/greenlandicesheet/v1.0.0.CrossRefGoogle Scholar
Thomsen, HH and 6 others (1997) The Nioghalvfjerdsfjorden glacier project, North-East Greenland: a study of ice sheet response to climatic change. Geology of Greenland Survey Bulletin 176, 95–103.CrossRefGoogle Scholar
Turton, JV, Mölg, T and Collier, E (2019b) NEGIS_WRF model output, Open Science Framework Repository (doi: 10.17605/OSF.IO/53E6Z ) In URL: https://osf.io/53e6z/ (last accessed on 09 May 2019).Google Scholar
Turton, JV, Mölg, T and Collier, E (2020) High-resolution (1 km) polar WRF output for 79°N glacier and the Northeast of Greenland from 2014–2018. Earth System Science Data 12(2), 1191–1202. doi: 10.5194/essd- 12-1191-2020CrossRefGoogle Scholar
Turton, JV, Mölg, T and van As, D (2019a) Atmospheric processes and climatological characteristics of the 79N glacier (Northeast Greenland). Monthly Weather Review 147(4), 1375–1394. doi: 10.1175/MWR- D-18-0366.1CrossRefGoogle Scholar
van As, D (2011) Warming, glacier melt and surface energy budget from weather station observations in the Melville Bay region of northwest Greenland. Journal of Glaciology 57(202), 208–220. doi: 10.3089/002214311796405898CrossRefGoogle Scholar
van de Wal, RSW and Oerlemans, J (1994) An energy balance model for the Greenland ice sheet. Global and Planetary Change 9(1–2), 115–131. doi: 10.10160921-8181(94)90011-6CrossRefGoogle Scholar
Vionnet, V and 7 others (2012) The detailed snowpack scheme Crocus and its implementation in SURFEX v7.2. Geoscientific Model Development 5(3), 773–791 doi: 10.1029/2010JD015013.CrossRefGoogle Scholar
Zeising, O, Neckel, N, Steinhage, D, Scheinert, M and Humbert, A (2020) Stake surface ablation/accumulation measurements at 79°N Glacier between 2017 and 2018 doi:10.1594/PANGAEA.922131 In URL: https://doi.pangaea.de/10.1594/PANGAEA.922131 (last accessed on 20 September 2020).Google Scholar
Figure 0

Fig. 1. Ice velocity of the GrIS from 14 December 2019 to 7 January 2020 (a) and the NEGIS from 6 October to 26 October 2016 (b) (Solgaard and Kusk, 2019). Site map of the 79N and the lower drainage basin of the NEGIS (c). The map in (c) shows the elevation taken from the PWRF output and presents the modeled domain. The isolines indicate 100 m (gray) and 1000 m (black) altitude steps. Stars show locations of the AWSs from PROMICE (Fausto and van As, 2019) and points show the location of surface ablation stakes (Zeising and others, 2020).

Figure 1

Table 1. List of input data from PWRF, AWS and stake observations including the time span of data usage and used variables

Figure 2

Table 2. Parameterization and parameters selection for the COSIPY-WRF run

Figure 3

Table 3. Evaluation statistics (correlation coefficient – R and root mean square error – RMSE) for COSIPY-WRF evaluation against PROMICE AWSs in the optimization and the evaluation period

Figure 4

Fig. 2. Model evaluation against the KPC-L of hourly surface temperature (a), hourly outgoing longwave radiation (b), daily shortwave albedo (c) and accumulated hourly surface height change relative to the first time step (d). The time period is 1 August 2014 00:00 to 31 December 2018 23:00 LT. The blue-dashed lines indicate the start and end time of the optimization period. Due to a measurement error of the SR50 sensor for the snow height change (stripping), the time before 1 July 2015 was excluded from evaluation for the snow height evaluation. Snow height can drop below zero due to a small initial snow layer at the observation site.

Figure 5

Fig. 3. Model evaluation against PROMICE the KPC-L AWS of hourly surface temperature (a), hourly outgoing longwave radiation (b), daily shortwave albedo (c) and accumulated hourly surface height change relative to the first time step (d). The time period is 1 August 2014 00:00 to 31 December 2018 23:00 LT. The blue-dashed lines indicate the start and end time of the optimization period. Snow height can drop below zero due to a small initial snow layer at the observation site.

Figure 6

Fig. 4. Comparison of the modeled glacier height change with the GPS measurements at the ablation stakes from AWI. At each site, four stakes were deployed to prevent data loss in case of the failure of one stake and to provide error estimates. The error bars indicate the SD of the measured values at each site. Coloring shows reference to the location of the measurement sites like in Figure 1 (red: ice margin, blue: grounding line, black: high altitude region).

Figure 7

Fig. 5. Spatial distribution of the annual SMB from September to August of the following year starting 2014 (a), 2015 (b), 2016 (c) and 2017 (d). The contours indicate 100 m (gray) and 1000 m (black) altitude steps taken from WRF.

Figure 8

Fig. 6. The ratio of surface meltwater that contributes to runoff, refreeze, evaporation or remains in the snow layer as liquid water content from September to August of the following year starting 2014 (a), 2015 (b), 2016 (c) and 2017 (d). The bottom panels show the contribution of SMB and IMB to the TMB of the same periods starting 2014 (e), 2015 (f), 2016 (g) and 2017 (h).

Figure 9

Fig. 7. Average vertical profile (altitude bands of 50 m) of the TMB (a and b) and SMB (c and d), and the components of the modeled SMB. This includes snowfall from PWRF (e and f), deposition (g and h), surface melt (i and j) and sublimation (k and l). The red line shows results for 2016–17 and the blue line shows results for 2017–18. The data are averaged over the accumulation periods (October 2016 to April 2017, and October 2017 to April 2018) in light blue boxes and over the ablation periods (May to September 2017 and May to September 2018) in brown boxes.

Figure 10

Fig. 8. Average vertical profile (altitude bands of 50 m) for the components of the SEB. This includes net shortwave radiation, SWnet (a and b), albedo (c and d), net longwave radiation, LWnet (e and f), sensible heat flux (g and h), latent heat flux (i and j) and melt energy (k and l) from COSIPY-WRF. The red line shows results for 2016–17 and the blue line shows results for 2017–18. The data are averaged over the accumulation periods (October 2016 to April 2017, and October 2017 to April 2018) in light blue boxes and over the ablation periods (May to September 2017 and May to September 2018) in brown boxes.

Figure 11

Table 4. GPS locations of the ablations stakes and the measurement results

Figure 12

Table 5. List of parameterization options available in v1.3 COSIPY (Sauter and others, 2020)

Figure 13

Table 6. List of key parameters and the range of values tested in model optimization

Figure 14

Fig. 9. Average vertical profiles (altitude bands of 50 m) for the surface temperature from COSIPY-WRF (a) and the mean climate conditions – 2 m air-temperature (b), 2 m wind speed (c) and relative humidity – from the PWRF output at the 79N region in July 2017 and July 2018. The surface temperature.

Figure 15

Fig. 10. Average vertical profile (altitude bands of 50 m) for the components of incoming shortwave radiation, SWin (a and b) and incoming longwave radiation, LWin (c and f) from COSIPY-WRF with energy flux density on the x-axis. The data are averaged over the accumulation periods (October 2016 to April 2017, and October 2017 to April 2018) and over the ablation periods (May to September 2017 and May to September 2018).