Modelling Glaciers in the HARMONIE-AROME NWP model

HARMONIE-AROME is a convection-permitting non-hydrostatic model that includes the multipurpose SURFEX surface model. It is developed for high resolution (1–3 km) weather forecasting and applied in a number of regions in Europe and the Mediterranean. A version of HARMONIE-AROME is also under development for regional climate modelling. Here we run HARMONIE-AROME for a domain over Greenland that includes a significant portion of the Greenland ice sheet. The model output reproduces temperature, wind speed and direction and relative humidity over the ice sheet well when compared with the observations from PROMICE automatic weather stations (AWS) operated within the model domain on the ice sheet (mean temperature bias 1.31± 3.6 K) but we identified a much lower bias (−0.16± 2.3 K) at PROMICE sites on days where melt does not occur at the ice sheet surface and is thus an artefact of the simplified surface scheme over glaciers in the existing HARMONIE-AROME operational set-up. The bias in summer time temperature also affects wind speed and direction as the dominant katabatic winds are caused by the cold ice surface and slope gradient. By setting an upper threshold to the surface temperature of the ice surface within SURFEX we show that the weather forecast error over the Greenland ice sheet can be reduced in summer when glacier ice is exposed. This improvement will facilitate accurate ice melt and run-off computations, important both for ice surface mass budget estimation and for commercial applications such as hydro-power forecasting. Furthermore, the HCLIM regional climate model derived from HARMONIE-AROME will need to accurately account for glacier surface processes in these regions in order to be used to accurately compute the surface mass budget of ice sheets and glaciers, a key goal of regional climate modelling studies in Greenland.


Introduction
The HARMONIE-AROME model is the current operational weather forecasting model used by the Danish Meteorological Institute and other weather forecasting institutions in Europe and North Africa (Termonia et al., 2017) for numerical weather prediction (Bengtsson et al., 2017).Operational HARMONIE-AROME numerical weather prediction (NWP) domains include both Greenland and Iceland but little or no model evaluation of HARMONIE-AROME has previously been performed over ice covered surfaces.
NWP models have frequently been adapted for climate modelling purposes (e.g.Christensen et al., 2001;Samuelsson et al., 2011;Van Meijgaard et al., 2008) including for regional climate modelling in Greenland where there is a particular focus on using models to reconstruct the surface mass budget of the Greenland ice sheet (e.g.Ettema et al., 2009;Fettweis et al., 2008;Langen et al., 2015Langen et al., , 2017;;Lucas-Picher et al., 2012;Rae et al., 2012).The HCLIM model, a climate variant of HARMONIE-AROME, though with a different dynamical core, is already in use for climate simulations within Europe (Lindstedt et al., 2015).A future target is to apply HCLIM more broadly to domains such as Greenland where glaciation is an important environmental characteristic, in order to accurately characterize the surface mass budget of the ice sheet and peripheral glaciers.
The surface mass budget, sometimes also called the climatic mass budget (Cogley et al., 2010) is the sum of the accumulation of precipitation and the amount of melt and run-off at the surface, taking into account the refreezing of liquid water in the surface layers of the snowpack.It is an important variable that drives both the dynamics of ice flow and determines in large part the contribution of mass to the oceans from glaciers (Church et al., 2013).Glacier surface properties are less important for NWP purposes due to the relatively small human population affected by glaciers and they have thus often been neglected in NWP model development.However, HARMONIE-AROME uses the SURFEX land surface model that contains parameterizations valid over snow covered surfaces.In this paper we examine the model performance over glacier surfaces in an operational NWP version of HARMONIE-AROME with the aim of assessing its utility as a climate model for applications in Greenland in the future, including in the HCLIM model format.

HARMONIE-AROME Model Configuration
The ALADIN-HIRLAM NWP system is a modelling framework with multiple options in terms of physics and dynamical schemes, it is thus very adaptable to study different processes in different locations as well as being designed to provide numerical weather predictions.In this project we use a recent version of the operational DMI-HARMONIE-AROME, which includes AROME physics, in a configuration for south Greenland, known as GLB (Yang et al., 2017).In addition, the corresponding single-column framework of HARMONIE-AROME, MUSC (Malardel et al., 2010;Rontu et al., 2016) has been used in an auxiliary study.
GLB is the updated version of DMI-HARMONIE (Yang et al., 2012) with source code based on the reference harmonie-38h1.2 as released by the HIRLAM research programme.GLB is configured on a model domain covering southwest Greenland (Fig. 1), with a horizontal grid of 2 km, 65 vertical levels and with a model top at 10 hPa.Main features of the forecast model are described in Seity et al.(2011) with some recent updates included as described in Bengtsson et al. (2017).The SURFEX externalised surface scheme (Masson et al., 2013) is used for surface data assimilation and the modelling of surface processes.In this study, we examine forecasts and observational data for 2014 and use summer 2015 for hindcast runs given the availability of relevant operational data archives and observational data.For the selected period, GLB cycling was run at 6 h intervals, with 6 h forecast integration for each cycle.Surface observations of temperature and humidity at coastal Greenland weather stations were assimilated but no observational data was included over the ice sheet.

Modifications to glacier surface scheme
We introduced a modification to the HARMONIE-AROME surface scheme that better simulates a glacier surface.The surface routine in SURFEX is known as ISBA (Interactions Soil -Biosphere -Atmosphere) and inside this routine the skin temperature is calculated from the energy budget at the surface.Where there is snow on the ground, this is allowed to melt, if the skin temperature is above freezing point until it has all disappeared.Within ISBA, we restricted the radiative surface temperature, PTS_RAD, such that the skin temperature is not permitted to rise above 273.15K in grid cells where glacier ice (termed "permanent snow" in the ECO-CLIMAP physiographic database used in SURFEX v.6) is present, even when snow is no longer there.As the skin temperature of ice cannot be higher than the melting point, this is a reasonable approximation when snow has disappeared.The fix does not affect snow covered ground or unglaciated areas.As a default setting, the depth of permanent snow is set at 10 m at the start.
A similar fixed skin temperature has been exploited in regional climate model surface schemes previously (for example, Ettema et al., 2009;Langen et al., 2015) and used to calculate the amount of energy available to melt both snow and glacier ice from the surface energy budget.Melt of snow and ice on glaciers that runs off the glacier surface is termed ablation and, in combination with accumulated snowfall, determines the mass budget of the glacier surface.HARMONIE-AROME does not yet contain a surface mass budget model but by introducing this correction, output from the model can be used to drive an offline model (for example, as in Langen et al., 2015Langen et al., , 2017)).However, it is important to remember that this fix does not adjust any of the other important parameters that determine melt rates over a glacier surface such as albedo and roughness lengths for turbulent energy fluxes (Fausto et al., 2016).This will be the focus of future work.

PROMICE Automatic Weather Station (AWS) Data
In order to assess the performance of HARMONIE-AROME over the ice sheet we examined model output and compared it with data from the PROMICE weather station observations (Van As et al., 2011) on the ice sheet at a number of locations around the ice sheet where the domain overlapped with the AWS sites that were active in 2014.
We here concentrate on locations (see Fig. 1) in southern Greenland, close to the village of Qassimiut, known as QAS_L (280 m a.s.l.) and QAS_U (900 m).Observations are maintained year round and in the observational records (going back to 2008) QAS_L has the highest observed melt rate in Greenland (Van As et al., 2016).We also include observations in southwest Greenland close to Nuuk at the sites NUK_L, NUK_N and NUK_U (530, 920 and 1120 m a.s.l.respectively).The AWS at these sites measure air temperature, wind speed and direction, relative humidity, downwards and upwards longwave and shortwave radiation, snow depth and ice height using a sonic ranger, ice temperatures and melt rates derived from pressure and temperature sensors drilled into the glacier ice.This means that the AWS data can be used to calculate the energy budget at the surface as well as other variables such as cloud cover and albedo that are not directly measured (Van As et al., 2011).The extreme operating environment and very high rates of melt, particularly at the lower sites (QAS_L and NUK_L) can lead to obvious errors in the data.In particular, the AWS may tilt during the ablation season due to surface melt but can be corrected using the inbuilt tilt-meter.We used a quality controlled and corrected version of the observational data provided by the PROMICE research group at GEUS, the Geological Survey of Denmark and Greenland (www.promice.org).

MODIS Albedo Data
We examine the surface properties of the glacier surface in the model and compared it to albedo derived from the MODIS satellite product, MOD10A1, developed and described by Stroeve et al. (2006).Since the MODIS product is based on optical data and thus subject to interference from clouds, these are filtered out using an automated algorithm and a time dependent interpolation scheme is used (Box et al., 2012).For this reason we show data from the first week of August to compare with the end of July.Data is available at a maximum resolution of 500 m and were provided by Jason Box (personal communication, 2017).and September in Figs. 2, 3 and 4, respectively.In order to compare the winds, the measured winds are scaled to 10 m from the assumption of a logarithmic wind profile as given in Eq. ( 1)

We
Here h boom is the height of the PROMICE station measurement boom, U (h boom ) is the wind measured at this height, U (10 m) is the wind scaled to 10 m height, and the roughness length z 0 is taken from HARMONIE-AROME to be 0.001 m for snow and ice.The model reproduces the observed magnitude and variability of the 2 m temperature, 2 m relative humidity, 10 m wind speed and 10 m wind direction for both the 12 and 24 h forecasts at QAS_L in January (Fig. 2) (for statistics see Fig. 6).The performance of the model is particularly impressive when taking into consideration the deep inversions that often form over the ice sheet in winter and which many regional climate models fail to capture (Rae et al., 2012).
However, in early June, the good fit of modelled temperature and wind speed and direction with observations, breaks down at QAS_L (Fig. 3).The modelled 2 m temperature is not only consistently much higher than the observed temperature, it also shows much larger variability, most likely due to the diurnal cycle being represented in the surface temperature.Similarly, the modelled wind speed becomes much lower compared to the observations and the wind direction also fits less well with the observations.Around this time the surface of the glacier becomes snow free in the model.Due to the model surface set-up this exposes a surface that is very dark and absorbing of radiation rather than one which uses the incoming energy to generate melt.Consequently the surface temperature is unrealistically high in the model and thus 2 m temperatures also start to have substantial swings over the diurnal cycle.This also affects the ability of the model to simulate katabatic winds, NWP simulations over glaciers in Iceland, have shown similar issues in HARMONIE-AROME (Nawri, 2015).
A similar situation occurs at the other PROMICE sites when the snow cover at the surface disappears in the model, though the timing of snow loss is different at the different locations.In September, when the melt season ends and snow starts to covers the ice surface again at QAS_L, the model again reproduces the meteorological observations reasonably as shown in Fig. 4.
In order to test that the model discrepancies occur when the surface is melting, we have plotted time series of the temperature and wind direction differences for each of the PROMICE stations within the GLB domain that we have model output data for in Fig. 5.
In Fig. 5 we show data where the observed surface temperature was less than 273 K separately (green dots), as these are periods of no surface melt.Positive temperature biases of more than 15 K are seen at times when surface melt occurs (the red dots in Fig. 5), while the temperature differences are more evenly distributed about 0 K when melting does not occur.Table 1.A summary of the verification statistics for both the 12 and 24 h HARMONIE-AROME forecasts is shown.The statistics are based on all available data during 2014 from the 5 PROMICE stations: NUK_L, NUK_N, NUK_U, QAS_L and QAS_U.In the rows labelled with "no melt" only data from months where no surface melt occurred are included.T is the 2 m temperature from the model, U is the 10 m wind from the model, and U dir is the 10 m wind direction from the model.A notable exception from this is during April and May at QAS_U, when the temperature deviations at are in the range from −13 to +9 K.During this time the boom height measurements are around 0.5 m, indicating that the measurement station was almost covered with snow.This makes it reasonable to assume that the large temperature deviations in this period arise from measurement uncertainties.The boom heights at the other stations were 1.5 m or more during 2014.The wind direction shown in the right panels of Fig. 5 show a similar increase in offset during periods of surface melt though it is important to note here that there are uncertain-ties due to the difference in height between the model and observed wind data.

Stat
In Fig. 6 the frequency distributions for the temperature differences are shown for all stations.The mean temperature differences (the biases) and the standard deviations are plotted in each subplot.The standard deviation of the Gaussian distribution is used.Assuming that the measurements are unbiased, this standard deviation is the root sum of the squares of the measurement uncertainty and the standard deviation of the model error.Of these the model errors cannot be expected to have a Gaussian distribution, however, we find that the standard deviation is still a useful and practical measure of the model quality.In the right hand column of Fig. 6 only data from months without melt are included.The improved distributions and statistics are clear to see.The temperature bias and standard deviation for the 12 h forecasts for all 2014 data from all 5 stations are +1.31 and 3.61 K, respectively.For the months without melt they are −0.16 and 2.23 K, respectively.In Table 1 the overall model verification statistics are shown for both the 12 and 24 h HARMONIE-AROME forecasts.
In Fig. 7 the frequency distributions for the wind speeds and the wind directions are shown for each station.The mean differences (the biases) and the standard deviations are plotted in each subplot.The wind speed bias and standard deviation for the 12 h forecasts for all 2014 data from all 5 stations are −0.15 and 2.58 m s −1 , respectively.The corresponding wind direction bias and standard deviation are 33.0 and 100.9 • , respectively.For the months without melt they are 10.4 and 89.4 • , respectively.This illustrates the effect of www.adv-sci-res.net/14/323/2017/Adv.Sci.Res., 14, 323-334, 2017  www.adv-sci-res.net/14/323/2017/Adv.Sci.Res., 14, 323-334, 2017  the melting glacier surface on the wind direction.In summer this katabatic wind will follow the surface gradient.When the model does not include the ice melt, the wind direction will follow the synoptic scale winds instead.These verification statistics are summarised in Table 1, where statistics are also shown for the 24 h forecasts.

Modified surface scheme results
We used single column experiments with MUSC at the QAS_L AWS location to test solutions to the surface temperature problem then ran the modified HARMONIE-AROME experimental simulations for July 2015 in the same domain as the operational output shown in Figs. 2 and 3.In Fig. 8, we show the skin and 2 m temperatures at QAS_L in observations and two versions of HARMONIE-AROME, the default version and the version with a modified ISBA scheme in SURFEX.The modified version of the ISBA fluxes scheme shows a significant improvement in skin temperature and 2 m temperature at the glacier surface compared to the unmodified version, though there is still a warm bias during the day and a cold bias, in the evening.The cold bias is likely due to the surface radiatively cooling during the night time, whereas the warm bias, though reduced, likely reflects other surface properties such as albedo which have not been corrected to more realistic values in this version of the model.The applied fix sets only an upper bound on the surface temperature and does not otherwise adjust the surface or sub-surface properties to be more realistic in terms of albedo, heat capacity or conductivity etc.We attribute the problems in unmodified HARMONIE-AROME to unrealistic surface parameter values and the lack of ice melt once the snow has disappeared from the glacier surface.Melting glacier ice is an efficient way to use the positive energy input from radiative and turbulent heat fluxes that would otherwise lead to an increase in surface temperature, the lack of melting thus means that energy that would otherwise be used for the phase change from ice to water instead causes the skin temperature to become much higher than realistic, with a consequent increase in modelled 2 m temperature.
Figure 9a shows that the modelled snow cover at QAS_L disappears on the 20-21 July, which is somewhat later than the observations suggest for this location.The observations (Fig. 9b) from the sonic ranger, mounted on the sensor boom at QAS_L, indicates that snow disappeared at this location on the 21 June, this seems to be confirmed by the broadband albedo at this site which attains a value of 0.6, the value commonly given to clean bare glacier ice on the same date.However, it is important to note that snow depth is not assimilated into HARMONIE-AROME and a full spin up to simulate an accurate snow depth was not performed.Instead, the model was started with a fixed single snow depth across the ice sheet, which is unrealistic compared to the observed snow depth at this site.The observations show however, that the albedo continues to decline throughout the melt season to around 0.2 on the 30 June and even lower by the end of the summer.These low values reflect the formation of cryoconite and the accumulation of dirt and impurities that melt out of the ice as well as microbial activity on the glacier surface (Stibal et al., 2015) but QAS_L lies in the lowest part of the ablation zone and has much lower albedo values than typical for most of the ablation zone in Greenland, (Ryan et al., 2017), as also shown in Fig. 10.The good match between observed and model albedo at this site therefore suggests that for most of the ablation zone the albedo will be underestimated once the surface is snow free.
It is not the aim of this paper to validate the modelled ice sheet albedo, but we compare the HARMONIE-AROME albedo with that derived from MODIS in Fig. 10 and this shows that model albedo is an area where future effort should be directed.In Fig. 10a, areas in the model which still have snow cover, have albedo values of 0.5 or above.These are similar to, though over much of the ice sheet still lower than, those indicated by the MODIS albedo product for the same period.(Fig. 10b).However, the snow free areas in the  model output are uniformly much lower in albedo than in the MODIS albedo output.In MODIS, and confirmed by direct on-ice observations (Moustafa et al., 2017), the slopes of the ice sheet show a clear gradient in albedo values to the very lowest areas of the ablation zone around the edges where the low albedo values in MODIS approach those given in the model.On average the bare ice albedo in the model is therefore too low compared to MODIS, reflecting that bare glacier ice is not a defined surface type within this version of HARMONIE-AROME.As HARMONIE-AROME albedo is given for both cloudy sky and clear sky conditions (Bengtsson et al., 2017), it is not quite the same as the MODIS MOD10A product shown here, where clouds are masked out and this also partly explains the difference in values for albedo over the accumulation zone of the ice sheet.

Conclusions
The NWP model HARMONIE-AROME shows great potential for accurate weather and climate simulation over glaciers in Greenland.We introduced a modification to HARMONIE-AROME that reduces model biases in 2 m temperature.This preliminary work also confirms that the model could be used successfully to model ice sheet surface mass budget at a very high resolution with sufficient modifications to the surface scheme.Planned modifications include the implementation of a glacier surface albedo scheme and glacier surface snowpack hydrology (as in, for example, Langen et al., 2017) and the calculation of surface mass budget based on accumulation of precipitation and ablation processes.Some of these processes will require careful initialisation and careful evaluation, for example, the accumulation of snow on the ice sheet and the new parameterizations and processes that need to be refined further such as biological effects on albedo and parameters that affect modelled sensible and latent heat fluxes (Fausto et al., 2016).The implementation of the recently upgraded full explicit snow scheme (Decharme et al., 2016) in SURFEX v8 will be required in HARMONIE-AROME, in particular to drive an albedo scheme with satellite derived mapped albedo values as in Langen et al. (2017).Nevertheless, the simple modification we show here represents a useful first step in simulating the weather over glaciers more accurately and shows that both HARMONIE-AROME and the spin-off HARMONIE-CLIMATE are viable for regional climate applications in Greenland.

Figure 1 .
Figure 1.(a) HARMONIE-AROME domain for southwest Greenland, GLB refers to the DMI operational model domain.(b) Map showing land (brown) and ice (white) in the domain.Black circles show locations of PROMICE automatic weather stations in the region, labels indicate locations where observational data were used in this analysis as not all stations had data available for the study period.

Figure 2 .
Figure 2. HARMONIE-AROME compared to PROMICE observations at QAS_L for January 2014, clockwise from top left: 2 m temperature, 10 m wind speed, 2 m relative humidity and 10 m wind direction (degrees from North) when compared with QAS_L data for January 2014.

Figure 3 .
Figure 3.As Fig. 2 but for June 2014.Note the divergence from observed 2 m temperature around the 12 June, at around the same time the wind speed becomes much more variable and the wind direction also starts to become more variable and to match less well with observations.

Figure 5 .
Figure 5.Time series plots of differences between 12 h HARMONIE-AROME forecasts of 2 m temperatures (a) and 10 m wind directions (b) relative to the measurements at the 5 PROMICE stations: NUK_L, NUK_N, NUK_U, QAS_L and QAS_U shown from top to bottom.The red points show all the data, while the green points show data for which the surface temperature was measured to be less than 273 K.All available data from 2014 are included.

Figure 6 .
Figure 6.Frequency distribution plots of differences between 12 h HARMONIE-AROME forecasts of 2 m temperatures relative to the measurements at the 5 PROMICE stations: NUK_L, NUK_N, NUK_U, QAS_L and QAS_U shown from top to bottom.In (a) all data from 2014 where the measurements and model output data overlap are included.In (b) only data for which the surface temperature was measured to be less than 273 K are included.

Figure 7 .
Figure 7. Frequency distribution plots of differences between 12 h HARMONIE-AROME forecasts of wind speeds diagnosed at the PROMICE station boom heights (a) and wind directions (b) relative to the measurements at the 5 PROMICE stations: NUK_L, NUK_N, NUK_U, QAS_L and QAS shown from top to bottom.

Figure 8 .
Figure 8.(a) Skin temperature (b) 2 m temperature at QAS_L from HARMONIE-AROME for the default and modified cases.These plots show the days in July after the snow has melted away in the model and bare ice is exposed.The skin temperature is significantly improved by the ISBA fix and the 2 m temperature has the warm bias substantially reduced.

Figure 9 .
Figure 9. (a) Harmonie modelled albedo at QAS_L July 2015; (b) PROMICE Albedo observed at QAS_L, the 21 June is the day when the glacier ice became exposed according to sonic ranger measurements but the glacier albedo continued to decline from here.

Figure 10 .
Figure 10.(a) Total albedo with cloud cover masked out for the 20 July from HARMONIE-AROME (b) 21 July albedo values from the MODIS satellite MOD10A archive.Note that HARMONIE-AROME uses a spectral albedo scheme but the figure shows total albedo, derived from the surface albedo scheme which includes a white sky, black sky and blue sky albedo.The much higher values over the snow-covered part of the ice sheet are largely due to the differences in the way the model and the satellite derived albedos are derived, but also show the importance of an accurate surface snow scheme over glaciers in HARMONIE-AROME.