Numerical simulation of extreme snowmelt observed at the SIGMA-A site, northwest Greenland, during summer 2012

Abstract

The surface energy balance (SEB) from 30 June to 14 July 2012 at site SIGMA (Snow Impurity and Glacial Microbe effects on abrupt warming in the Arctic)-A, (78${}^{\circ }$03${}^{\prime }$ N, 67${}^{\circ }$38${}^{\prime }$ W; 1490 m a.s.l.) on the northwest Greenland Ice Sheet (GrIS) was investigated by using in situ atmospheric and snow measurements as well as numerical modeling with a one-dimensional multi-layered physical snowpack model called SMAP (Snow Metamorphism and Albedo Process). At SIGMA-A, remarkable near-surface snowmelt and continuous heavy rainfall (accumulated precipitation between 10 and 14 July was estimated to be 100 mm) were observed after 10 July 2012. Application of the SMAP model to the GrIS snowpack was evaluated based on the snow temperature profile, snow surface temperature, surface snow grain size, and shortwave albedo, all of which the model simulated reasonably well. Above all, the fact that the SMAP model successfully reproduced frequently observed rapid increases in snow albedo under cloudy conditions highlights the advantage of the physically based snow albedo model (PBSAM) incorporated in the SMAP model. Using such data and model, we estimated the SEB at SIGMA-A from 30 June to 14 July 2012. Radiation-related fluxes were obtained from in situ measurements, whereas other fluxes were calculated with the SMAP model. By examining the components of the SEB, we determined that low-level clouds accompanied by a significant temperature increase played an important role in the melt event observed at SIGMA-A. These conditions induced a remarkable surface heating via cloud radiative forcing in the polar region.

Introduction

Snow and ice on the Greenland ice sheet (GrIS) experienced a record near-surface melt extent in summer 2012 (Nghiem et al., 2012; Tedesco et al., 2013; Hall et al., 2013; Hanna et al., 2014; Bennartz et al., 2013). The physical conditions of the atmosphere and the snow (as well as ice) surface during summer 2012 have been gradually clarified. The most notable feature of the event is that the extent of surface melt was the largest in the satellite era (March 2000 to the present). According to several studies that used satellite data (Nghiem et al., 2012; Tedesco et al., 2013; Hall et al., 2013), at least 95 % of the entire surface of the GrIS melted during the period 30 June to 14 July 2012. At present, it is hypothesized that a significant temperature increase over the GrIS may have been a major cause of the record near-surface melt event. Hall et al. (2013) examined the ice surface temperature (IST) of the GrIS derived from the Moderate-resolution Imaging Spectroradiometer (MODIS) and concluded that the 2012 summer was the warmest (IST: $-$6.38 $\pm$ 3.98 ${}^{\circ }$C) in the MODIS record. Bennartz et al. (2013) have demonstrated that low-level clouds consisting of liquid-water droplets played a key part in the melt event by increasing near-surface temperature via radiative effects. Nghiem et al. (2012) pointed out the existence of an anomalous ridge of warm air that could be identified by a 500 hPa height anomaly, and they concluded that the ridge acted as a strong heat dome that became stagnant over the GrIS during the period. Hanna et al. (2014) concluded that the extreme melt was forced mainly by atmospheric conditions linked with changes in the summer North Atlantic Oscillation (NAO), Greenland Blocking Index, and polar jet stream. The combination of these conditions favored southerly warm air advection along the western coast. According to Fettweis et al. (2013), the frequency of occurrence of anticyclones centered over the GrIS at the surface and at 500 hPa has doubled since the end of the 1990s. They associated this increased frequency with an increase of negative phases of the NAO, which induce more frequent southerly warm air advection along the western Greenland coast and over the neighboring Canadian Arctic Archipelago.

However, a temperature increase cannot induce surface melt if the surface temperature, which is physically controlled by the surface energy balance (SEB), is below 0 ${}^{\circ }$C. The SEB is equal to the sum of the net shortwave radiant flux, net long-wave radiant flux, sensible heat flux, latent heat flux, heat supply by rainfall, and subsurface conductive heat flux (Armstrong and Brun, 2008). In this study, these fluxes are defined to be positive when they are directed towards the snow surface, and a positive sum of these fluxes (net energy flux) induces surface melt only if the surface temperature equals 0 ${}^{\circ }$C. A temperature increase raises the net energy flux mainly by affecting the sensible heat flux, although other energy fluxes, especially net long-wave radiant flux, can of course be affected by a temperature increase. Therefore, it is necessary to pay close attention to other fluxes in addition to the sensible heat flux when surface melt is investigated.

Several attempts that focus on the summer GrIS SEB characteristics have been performed. Presented results show that the net shortwave radiant flux is the main contributor for the surface heating in general; however, detailed characteristics vary from place to place and differ from year to year. Greuell and Konzelmann (1994) unveiled temporal changes in SEB at the ETH Camp (69${}^{\circ }$34${}^{\prime }$ N, 49${}^{\circ }$18${}^{\prime }$ W, 1155 m a.s.l.), west Greenland, during the 1990 summer (June, July, and August). During this summer, average net shortwave radiant flux (82 W m${}^{-\mathrm{2}}$) and sensible heat flux (34 W m${}^{-\mathrm{2}}$) acted to heat the surface, while average net long-wave radiant flux ($-$54 W m${}^{-\mathrm{2}}$) and latent heat flux ($-$28 W m${}^{-\mathrm{2}}$) played a role in cooling the surface. SEB characteristics during the 1991 summer at the same place was presented by Ohmura et al. (1994). According to their results, the absolute values of each dominant SEB component averaged in this summer decreased obviously from the 1990 summer (net shortwave radiant flux was 65 W m${}^{-\mathrm{2}}$; net long-wave radiant flux was $-$44 W m${}^{-\mathrm{2}}$; sensible heat flux was 16 W m${}^{-\mathrm{2}}$; and latent heat flux was $-$6 W m${}^{-\mathrm{2}}$). Summer SEB characteristics at the higher place on the GrIS was described by Cullen and Steffen (2001). They demonstrated that average turbulent (sensible and latent) heat fluxes at Summit (72${}^{\circ }$58${}^{\prime }$ N, 38${}^{\circ }$51${}^{\prime }$ W, 3203 m a.s.l.) during 21 June to 6 July 2000 was small (4 and 3 W m${}^{-\mathrm{2}}$, respectively), while average net shortwave radiant flux (82 W m${}^{-\mathrm{2}}$) and net long-wave radiant flux ($-$68 W m${}^{-\mathrm{2}}$) were comparable with previous results at ETH Camp. van den Broeke et al. (2011) presented long-term records of SEB at three automated weather stations (AWSs) situated along the K-transect, a stake array in southwest Greenland that extends from the ice margin to 1850 m a.s.l. They demonstrated that the temperature and moisture contrasts between ambient atmosphere and (melting) ice surface are less pronounced higher on the ice sheet, resulting in smaller summertime values of turbulent heat fluxes and net long-wave radiant flux at the higher elevations.

In the present study, our goal was to understand the record near-surface melt event that occurred over the GrIS in 2012 by investigating snow–atmosphere interactions from the standpoint of the SEB. For this purpose, we employed data obtained from in situ measurements made at a site named SIGMA (Snow Impurity and Glacial Microbe effects on abrupt warming in the Arctic)-A (Aoki et al., 2014a) in northwest Greenland during the 2012 summer intensive observation period (30 June to 14 July, hereafter IOP). During the latter half of the IOP, which coincided with the period of the record near-surface melt event reported by Nghiem et al. (2012), we observed remarkable melting of the near-surface snowpack and encountered continuous rainfall (Aoki et al., 2014b; Yamaguchi et al., 2014). In the process of calculating the SEB, radiation-related terms could be determined from in situ measurements (described in Sect. 2); however, the sensible heat flux, latent heat flux, and heat supply by rainfall, which were not measured directly at the SIGMA-A site, had to be calculated in some way. In this study, we therefore calculated these fluxes by using the physical snowpack model SMAP (Snow Metamorphism and Albedo Process), which was originally developed and evaluated for seasonal snowpack (Niwano et al., 2012, 2014). Because this study was the first attempt to apply the SMAP model in a polar region, we carried out a preliminary examination of various aspects of the capacity of the SMAP model. This model validation study was expected to increase the reliability both of the SMAP model itself and of the Meteorological Research Institute Earth System Model version1 (MRI-ESM1) (Yukimoto et al., 2011), in which the SMAP model is used to calculate energy exchanges between the snow surface and atmosphere in the global cryosphere.

We begin by presenting an outline of our in situ measurements conducted at SIGMA-A during summer 2012. Then we describe the modification of the SMAP model that enabled its use to simulate polar conditions. We verified the accuracy of the model with respect to the snow temperature profile, snow surface temperature, surface optically equivalent snow grain radius, and shortwave albedo. Finally, we used the SMAP model to explore the SEB during the IOP at the SIGMA-A site by diagnosing the turbulent heat fluxes and heat supplied by rainfall.

Field measurements at the SIGMA-A site

In this section we report the in situ data acquired during the IOP at site SIGMA-A (78${}^{\circ }$03${}^{\prime }$ N, 67${}^{\circ }$38${}^{\prime }$ W, 1490 m a.s.l.) (Aoki et al., 2014a) in northwest Greenland (Fig. 1), and we report the observed meteorological and snow conditions. An AWS was newly installed at site SIGMA-A on 29 June 2012 (Aoki et al., 2014a). In the northwest GrIS, two GC-Net AWS sites exist (Steffen and Box, 2001): one is the Humboldt site (78${}^{\circ }$32${}^{\prime }$ N, 56${}^{\circ }$50${}^{\prime }$ W, 1995 m a.s.l.) and the other is the GITS site (77${}^{\circ }$08${}^{\prime }$ N, 61${}^{\circ }$02${}^{\prime }$ W, 1887 m a.s.l.). Steffen and Box (2001) presented monthly mean temperatures at these stations and also demonstrated that the mean temperature lapse rate over GrIS in summer was 0.4 ${}^{\circ }$C per 100 m. This information allows us to estimate the average near-surface air temperature during July at the SIGMA-A site, resulting in a value of $-$6.5 ${}^{\circ }$C. This region of the northwestern GrIS is now being carefully monitored, because since around late 2005 the ice loss over southern portions of the GrIS has been spreading rapidly northward along the northwest coast (Khan et al., 2010). Dates and times throughout this paper are expressed in terms of the universal time clock (UTC); the difference between local time (LT) at site SIGMA-A and UTC is $-$2 h.

Meteorological observations

The AWS measured air temperature and relative humidity with respect to water with a humidity and temperature probe (HMP155, Vaisala, Finland) protected from direct solar irradiance by a naturally aspirated 14-plate Gill radiation shield (41005, Young, USA) at 3.0 and 6.0 m above the snow surface on 29 June. According to van den Broeke et al. (2009), radiation errors in the unventilated air temperature measurement can be up to 3 ${}^{\circ }$C under conditions of low wind speed and high insolation. However, no corrections were made to measured air temperature in the present study, because low wind speeds ($<$ 1.0 m s${}^{-\mathrm{1}}$) were rarely observed at SIGMA-A during IOP, as noted later in this section. As for relative humidity with respect to water, we converted it into relative humidity with respect to ice when air temperature was below 0 ${}^{\circ }$C and performed the correction presented by Anderson (1994). The AWS also measured wind speed and direction with a wind monitor (05103, Young, USA) at 3.0 and 6.0 m above the surface. The AWS was also equipped with a digital barometer (PTB210, Vaisala, Finland) and an ultrasonic distance sensor (SR50, Campbell, USA) for snow depth monitoring. We used a net radiometer (CNR4, Kipp & Zonen, the Netherlands) to measure downward and upward shortwave (wavelength $\mathit{\lambda }$ $=$ 0.305–2.8 $\mathrm{\mu }$m) and long-wave ($\mathit{\lambda }$ $=$ 4.5–42 $\mathrm{\mu }$m) radiant fluxes. Downward and upward near-infrared ($\mathit{\lambda }$ $=$ 0.715–2.8 $\mathrm{\mu }$m) radiant fluxes were observed by pyranometers (CMP6, Kipp & Zonen, the Netherlands) with RG715 cutoff filter domes. Measured data were sampled and stored in a data logger (CR1000, Campbell, USA) every 1 min. Aoki et al. (2014a) have provided a more detailed description of the AWS. In the present study, we used 5 min average values of AWS data, but to highlight important temporal evolution the SEB data reported after Sect. 5 are 1 h averages.

The location of site SIGMA-A together with ice sheet surface elevation contours (m a.s.l.). Height interval is 100 m.

[Figure omitted. See PDF]

In addition to the automatic measurements, we performed manual spectral albedo measurements using a field spectroradiometer (FieldSpec 3, Analytical Spectral Devices, Inc., USA) ($\mathit{\lambda }$ $=$ 0.35–2.5 $\mathrm{\mu }$m). The details of the measurement procedure were based on the method of Aoki et al. (2000). Use of the obtained spectral albedo in a GSAF (ground-based spectral radiometer system for albedo and flux; Kuchiki et al., 2009) channel with a $\mathit{\lambda }$ of 1.23 $\mathrm{\mu }$m allowed us to retrieve the surface optically equivalent snow grain radius (hereafter simply referred to as the “snow grain size” unless otherwise stated) with the GSAF algorithm developed by Kuchiki et al. (2009). In the latest version of this algorithm, the model of snow grain shape can be arbitrarily chosen to be a sphere, spheroid, or Voronoi aggregate. In this study we selected the Voronoi aggregate model developed by Ishimoto et al. (2012), because this choice led to the best agreement with manually measured in situ surface snow grain size in the case of our data (see Sect. 2.2). We henceforth refer to these data as snow grain size retrieved by ground-based remote sensing (GRS). We used GRS data for model evaluation (Sect. 4.3).

Meteorological conditions measured at SIGMA-A during the IOP (30 June to 14 July 2012): (a) air temperature and relative humidity with respect to water observed at a nominal height of 3.0 m above the snow surface; (b) wind speed and direction observed at 3.0 m above the snow surface on 29 June 2012 (wind direction is denoted in degrees that indicate the direction from which the wind is blowing; the degrees increase from 0 to 360 as the direction is rotated clockwise from north); (c) downward shortwave (SW${}^{\downarrow }$), near-infrared (NIR${}^{\downarrow }$), long-wave (LW${}^{\downarrow }$) radiant fluxes, and the near-infrared fraction of the SW${}^{\downarrow }$ ; and (d) accumulated precipitation obtained from 3-hourly ERA-Interim reanalysis data (ERA-Int), where the black line denotes the original amount of precipitation and the red line reflects the correction based on in situ measurements of accumulated precipitation measured during the period from 20:00 UTC on 12 July 2012 to 11:00 UTC on 14 July 2012.

[Figure omitted. See PDF]

Figure 2 presents time series of meteorological conditions measured with the AWS during the IOP. Until 9 July, air temperature at 3.0 m above the surface was already high and often exceeded 0 ${}^{\circ }$C in the daytime. The time interval from 10 July until the end of the IOP coincided with the record near-surface melt event period reported by Nghiem et al. (2012); during that time air temperature increased slightly and remained above 0 ${}^{\circ }$C, which is much larger than the estimated possible average air temperature at the SIGMA-A site: $-$6.5 ${}^{\circ }$C (Sect. 2) continuously (Fig. 2a). The relative humidity with respect to water recorded at 3.0 m above the surface was related to the high air temperature and tended to be high, the values sometimes approaching 100 % (Fig. 2a). At site SIGMA-A, low wind speeds ($<$ 1.0 m s${}^{-\mathrm{1}}$) were rarely observed at 3.0 m above the surface during the IOP (Fig. 2b). The gradual day-by-day rise of air temperature was accompanied by a strong southeast wind with a speed sometimes over 10 m s${}^{-\mathrm{1}}$. Radiation parameters (Fig. 2c) reflect sky conditions. For example, until 3 July and around 9 July, we observed clear-sky conditions on site. During the clear-sky period, the downward shortwave and near-infrared radiant fluxes showed an obvious diurnal cycle, whereas the downward long-wave radiant flux was almost constant. In contrast, cloudy-sky conditions (subjectively observed prevailing cloud types were Cumulus, Stratocumulus, and sometimes Stratus) on 3, 4, 6, 7, and 8 July and after 10 July were characterized by reduced downward shortwave and near-infrared radiant fluxes and increased downward long-wave radiant flux. During such conditions, the near-infrared fraction of the downward shortwave radiant flux decreased to around 0.3, while it was around 0.5 under clear-sky conditions. Neff et al. (2014) examined synoptic-scale atmospheric conditions over the GrIS during July 2012 from various aspects and summarized notable features as follows: (1) warm air originating from a record North American heat wave (the North American drought of 2012 was the worst since 1895), (2) transitions in the Arctic Oscillation, (3) transport of water vapor via an atmospheric river over the Atlantic to Greenland, and (4) the presence of warm ocean waters south of Greenland. Bonne et al. (2015) clearly showed that moist air mass was advected northward following a narrow band reaching southern Greenland and then it moved northward along the western Greenland coast around 9 July. Observed features of above-mentioned meteorological properties during the IOP at the SIGMA-A site are consistent with these large-scale atmospheric conditions.

At the SIGMA-A site, no rainfall was observed on site until 9 July. After 10 July, we encountered continuous rainfall until 14 July (Aoki et al., 2014b; Yamaguchi et al., 2014). Because precipitation is one of the critical parameters for a snowpack model that is used to calculate mass balance of a snowpack, the total amount of precipitation during 10–14 July was estimated on the basis of precipitation collected with a bucket during the period from 12 July, 20:00 UTC, to 14 July, 11:00 UTC, and on the 3-hourly European Centre for Medium-Range Weather Forecasts (ECMWF) Interim reanalysis (ERA-Interim) data (Dee et al., 2011). At first, the precipitation that accumulated in the bucket from 12 July, 20:00 UTC, to 14 July, 11:00 UTC, and the 3-hourly ERA-Interim reanalysis data were compared. For this comparison, ERA-Interim 1 h accumulated precipitation was simply estimated to be one-third of the 3-hourly accumulated precipitation. The comparison indicated that the accumulated precipitation obtained from the ERA-Interim reanalysis data was lower by a factor of $\mathrm{1}/4.9$. The most prevailing reason for this discrepancy was not a misrepresentation of the true area of rainfall but just underestimation by the ERA-Interim reanalysis (Fig. 3). As demonstrated by Chen et al. (2011), ERA-Interim data originally tend to underestimate annual accumulation in the area north of 68${}^{\circ }$ N, even though it shows close spatial pattern of accumulation to the observations over the whole area of GrIS. In addition, it is possible that insufficient horizontal resolution (0.75${}^{\circ }$) and a hydrostatic atmospheric model, which cannot reproduce a short-time mesoscale convective system realistically in general, might have caused the large discrepancy. Finally, the precipitation that accumulated between 10 and 14 July, 100 mm, was obtained by multiplying ERA-Interim reanalysis data by a factor of 4.9 (Fig. 2d). According to Ohmura and Reeh (1991), annual total precipitation near the SIGMA-A site is extrapolated to be around 200–300 mm w.e. The estimated total precipitation during this event can account for more than 30–50 % of the annual total precipitation. The procedure to input this information into the SMAP model is described in Sect. 3.3.

Accumulated precipitation from 00:00 UTC on 10 July 2012 to 12:00 UTC on 14 July 2012 obtained from the 3-hourly ERA-Interim reanalysis data. Black colored circle indicates the location of site SIGMA-A.

[Figure omitted. See PDF]

Snow pit measurements

Over the course of the IOP, we performed snow pit observations in the local morning (around 09:30 LT) in order to characterize the physical profiles of the snow (Yamaguchi et al., 2014). We also collected near-surface snow samples to measure the mass concentrations of light-absorbing snow impurities such as black carbon (BC) and dust (Aoki et al., 2014b).

The near-surface layer (NSL) of the snow was recognizable by quite a thick bottom ice layer. We focused on this near-surface layer in the snow pit observations, which were performed every day from 30 June to 13 July except for 11 and 12 July, when there was heavy rain at the SIGMA-A site. The thickness of the layer (above the ice layer) was 88 cm on 30 June. At present, the NSL has not been determined to be the latest annual layer because of the lack of justification. Measured properties included profiles of snow grain shape, snow temperature (including bottom ice layer temperature), snow density, volumetric water content, and geometric snow grain size. In addition, the snow grain size of the top 10 cm was measured manually after 1 July. In this snow grain size measurement process, we followed the definition of “$r_{\mathrm{2}}$” presented by Aoki et al. (2000, 2003, 2007): half the branch width of dendrites or half the width of the narrower portion of broken crystals. Figure 4 shows the temporal evolution of the observed snow grain shape profile in the NSL at the SIGMA-A site from 30 June to 13 July 2012. On the whole, until 9 July this layer could be roughly divided into three parts: a bottom layer composed of depth hoar, a middle layer of rounded grains (sometimes including faceted crystals), and a top layer composed of melt forms, including an ice layer of variable thickness. After 10 July, a wet snow layer expanded from the top to the bottom as the air temperature increased slightly (Fig. 2a), and the rainfall event occurred (Fig. 2d). The measured snow temperature profile is presented in Sect. 4.1 and compared to the results simulated by the SMAP model. In this comparison, a total of 221 profiles (after rejecting unphysical data) were available. As for the other snow properties, we used only the data that we obtained on 30 June to specify the initial physical conditions of the snowpack at the SIGMA-A site for the SMAP model simulations. The procedure is briefly explained in Sect. 3.3.

Temporal evolution of observed main snow grain shape profile within the NSL, which was defined by a thick bottom ice layer, at SIGMA-A during the IOP. Vertical axis (relative snow depth) denotes the height above the bottom thick ice layer. Characters and colors indicating snow grain shape follow the definition by Fierz et al. (2009). In sequence from the left they denote precipitation particles, decomposing and fragmented precipitation particles, rounded grains, faceted crystals, depth hoar, surface hoar, melt forms, and ice layer.

[Figure omitted. See PDF]

Snow samples to measure the mass concentrations of snow impurities were collected from the top 2 cm and the 2–10 cm layer every other day from 30 June to 12 July (Aoki et al., 2014b). To obtain the mass concentrations of elemental carbon (EC), the samples were melted and filtered on site, and the filters were weighed and assayed with a carbon analyzer (Lab OC-EC Aerosol Analyzer, Sunset Laboratory Inc., USA). In this study we assumed, following Aoki et al. (2011), that EC was equivalent to BC. The dust concentrations were determined by subtracting total carbon concentrations from the total impurity concentrations, which were estimated from the difference between the dry weights of the filters before and after filtering (Aoki et al., 2003, 2007). The detailed analytical procedure is precisely described by Kuchiki et al. (2009) and Aoki et al. (2011, 2014b). Figure 5 shows the temporal changes in the mass concentrations of EC and dust in the near-surface snowpack at the SIGMA-A site during the IOP; the EC concentrations increased rapidly after 6 July. Aoki et al. (2014b) have discussed the reasons for this obvious increase, and they concluded that effects of sublimation/evaporation and snowmelt amplification due to the low scavenging efficiency of meltwater played an important role.

Numerical simulation with the SMAP model

SMAP model overview

The multilayered physical snowpack model SMAP calculates mass and energy balances of the snowpack by taking snow settlement, phase changes, water percolation, and snow metamorphism into account. The most distinguishing characteristic of the SMAP model is the fact that it incorporates the physically based snow albedo model (PBSAM) developed by Aoki et al. (2011). PBSAM calculates snow albedo and the solar heating profile in the snowpack by explicitly considering effects of snow grain size and light-absorbing snow impurities such as BC and mineral dust. In default configuration, the SMAP model requires precipitation (partitioned in the model into snow and rain by using the algorithm to calculate snow : rain ratios as a function of wet bulb temperature; Yamazaki, 2001), air pressure, wind speed, air temperature, relative humidity, downward ultraviolet (UV)–visible and near-infrared radiant fluxes, the diffuse components of UV–visible and near-infrared radiant fluxes, downward long-wave radiant flux, subsurface heat flux, and the mass concentrations of snow impurities (BC and dust) (Niwano et al., 2012). In the present study, the diffuse components of the UV–visible and near-infrared radiant fluxes were not available, as explained in Sect. 2.1. Under these circumstances, SMAP calculates these fluxes by using the scheme of Goudriaan (1977) together with internally diagnosed cloud fraction calculated from air temperature and simulated net long-wave radiant fluxes at the snow surface. The procedure for obtaining cloud fraction is based on the approach of van den Broeke et al. (2004, 2006). Niwano et al. (2014) used data obtained during the winters of 2007–2009 at Sapporo to investigate the effectiveness of the process, and they demonstrated that the accuracy of the simulated snow depth and snow surface temperature were comparable in magnitude to the accuracy of the default configuration.

Measured mass concentrations of elemental carbon (EC) and dust in the top 2 cm of snow layers and the subsurface (2–10 cm) layers at the SIGMA-A site during the IOP. Dashed lines indicate periods when there were no measurements.

[Figure omitted. See PDF]

The main governing equation of the SMAP model is a one-dimensional energy balance equation that takes solar heating of the snowpack and melt–freeze cycles into account (Niwano et al., 2012). This equation is approximated and solved with the Crank–Nicolson finite difference implicit method. In this numerical solution, the SMAP model assumes each model snow layer to have a thickness $d$ that is allowed to range between $d_{\text{min}}$ and $d_{\text{max}}$. The values of $d_{\text{min}}$ and $d_{\text{max}}$ used in the present study are discussed in Sect. 3.3. Solar heating of the snowpack was calculated with a component of PBSAM. Input requirements for this component include internally calculated profiles of snow grain size and snow water equivalent, as well as the mass concentrations of snow impurities given externally from in situ measurements or host global or regional circulation models. Snow grain size was calculated by employing a model geometry that envisions two spherical ice particles connected by a neck (Lehning et al., 2002). By calculating the specific surface area of snow per unit volume with the nonspherical model geometry, SMAP obtains snow grain size. Temporal evolution of snow grain size is governed by the following three types of metamorphism regimes: equi-temperature metamorphism, temperature-gradient metamorphism, and wet snow metamorphism. The formulation of these metamorphism laws is based on results from the SNOWPACK model (Lehning et al., 2002). Furthermore, the SMAP model implicitly takes into account the effects of snow metamorphism under alternating temperature gradients (small rounded grains can be formed even when the temperature gradient is large when the sign of the temperature gradient changes with a 24 h cycle) (Pinzer and Schneebeli, 2009) by forcing temperature gradient metamorphism, which generally induces rapid grain growth, not to occur in the top 20 cm of the each model layer.

Recently, the SMAP model was updated with respect to water movement in the snowpack, snow settlement, and turbulent heat exchanges under very stable conditions (Niwano et al., 2014). Water movement in the snowpack is now governed by the Richards equation (Richards, 1931), which takes into consideration Darcy's law; hydraulic diffusivity and hydraulic conductivity are calculated by the van Genuchten model (van Genuchten, 1980) adapted to snow (Shimizu, 1970; Hirashima et al., 2010; Yamaguchi et al., 2010, 2012). As for snow settlement, the updated version of the SMAP model calculates the viscosity coefficient of snow using the scheme presented by Vionnet et al. (2012), which was developed to improve the performance of the Crocus model (Brun et al., 1989, 1992; Vionnet et al., 2012; Carmagnola et al., 2014). To ensure small but non-zero turbulent heat exchanges under very stable conditions, we set an upper bound of 0.1 on the Richardson number. During IOP, low wind speeds were rarely observed at the SIGMA-A site, as mentioned in Sect. 2.1. As a result, preliminary numerical simulation revealed that the impact of the upper bound was not clear during IOP; however, it is still set in this study.

Adaption of the SMAP model to the Greenland snowpack

In this section we discuss the adaptation of the model that enabled the SMAP model to be suitable for polar (especially Greenland) snowpack simulations. To simulate the temporal evolution of snow temperature accurately by solving the one-dimensional energy balance equation in the snowpack, the thermal conductivity of snow $k_{\text{eff}}$ (W m${}^{-\mathrm{1}}$ k${}^{-\mathrm{1}}$) should be assigned an appropriate value. In this study, we employed the parameterization of Anderson (1976) that describes $k_{\text{eff}}$ as a function of snow density ${\mathit{\rho }}_{\text{s}}$ (kg m${}^{-\mathrm{3}}$): $$k_{\text{eff}}=0.021+2.5{\left(\frac{{\mathit{\rho }}_{\text{s}}}{1000}\right)}^{\mathrm{2}}.$$ This equation is widely used in polar regions (e.g., van den Broeke et al., 2005; Kuipers Munneke et al., 2009; van As, 2011; Brun et al., 2011).

The snow surface roughness length for momentum $z_{\mathrm{0}}$ affects turbulent heat exchanges between the snow surface and the atmosphere. Brock et al. (2006) compiled published data of $z_{\mathrm{0}}$ and demonstrated that the values were quite dependent on surface conditions and places. Considering this, we assumed the following constant values for $z_{\mathrm{0}}$ (mm) presented by Greuell and Konzelmann (1994): $$z_{\mathrm{0}}=\left\{\begin{array}{ll}0.12& \text{for}\text{snow}\text{before}\text{melting}\\ 1.3& \text{for}\text{snow}\text{after}\text{melting}.\end{array}\right.$$ They simulated the mass balance and englacial temperature at the ETH Camp, west Greenland, during the 1990 summer melting period with these values. During IOP, the surface condition at SIGMA-A was smooth. van den Broeke et al. (2009) presented time series of $z_{\mathrm{0}}$ at sites S5 (490 m a.s.l.), S6 (1020 m a.s.l.), and S9 (1520 m a.s.l.) on the K-transect in southwest Greenland during August 2003 to August 2007. The above-mentioned values for $z_{\mathrm{0}}$ (Eq. 2) fall in the range of their result at S9, where surface condition was relatively smooth. The uncertainty of the SMAP model calculated SEB caused by the choice of $z_{\mathrm{0}}$ is investigated in Sect. 5. In the SMAP model, roughness lengths for heat and moisture are calculated following Andreas (1987) as explained by Niwano et al. (2012).

Model configuration

In the present study we calculated the temporal evolution of snow physical states in the NSL (see Sect. 2.2) at the SIGMA-A site from 30 June 2012, 16:45 UTC, to 14 July 2012, 12:00 UTC. To perform detailed model simulations we divided the NSL into several thin model layers. According to Dadic et al. (2008), the recommended model vertical resolution for an Arctic snowpack is on the order of 10 mm, because there is a strong feedback between small-scale snow structure and snow temperature that should be simulated accurately by snowpack models. Based on this consideration, the SMAP model vertical layer thickness (see Sect. 3.1) used in this study ranged between $d_{\text{min}}$ $=$ 2 mm and $d_{\text{max}}$ $=$ 6 mm. Accordingly, it became necessary for the time step of the numerical integration $\mathrm{\Delta }t$ to be shortened. Thus, $\mathrm{\Delta }t$ was set to 30 s in the present study.

Input parameters used to drive the SMAP model in the present study were air pressure; air temperature and relative humidity at 3.0 m above the snow surface (Fig. 2a); wind speed at 3.0 m above the snow surface (Fig. 2b); downward shortwave, near-infrared, and long-wave radiant fluxes (Fig. 2c); accumulated precipitation (Fig. 2d); the temperature of the thick bottom ice layer in the NSL (lower boundary condition); and the mass concentrations of snow impurities (Fig. 5). In addition, the emissivity of the snow surface ${\mathit{\epsilon }}_{\text{s}}$ was assumed to be 0.98 (Armstrong and Brun, 2008; van As, 2011) throughout this study.

The time interval for meteorological properties was selected to be 5 min, whereas daily values were used for other properties. The corrected 3-hourly accumulated precipitation data were divided equally into each time interval. Regarding the measured mass concentrations of snow impurities (Fig. 5), we followed the same method used by Niwano et al. (2012), who equated the values in the top 2 cm of the model layers of the snowpack to the corresponding observed values and assigned the observed values in the 2–10 cm depth interval equally to the lower model layers. During the period when measurements of the mass concentrations of snow impurities were unavailable (dashed lines in Fig. 5), we used the values measured at the nearest point in time. The initial physical states of the snowpack in the NSL on 30 June (profiles of snow grain shape, snow temperature, snow density, volumetric water content, and geometric snow grain size) were taken from the snow pit observations conducted on 30 June 2012, 16:45 UTC (Yamaguchi et al., 2014). These observations were distributed within the model layers, and the layer thickness $d$ was in all cases set to 5 mm in the NSL at that time. Because the SMAP model calculates shortwave albedo as a function of snow grain size and the mass concentrations of snow impurities, information regarding an optically equivalent snow grain size profile was necessary. In the present study, an optically equivalent profile was obtained by multiplying the geometric profile by a factor that produced agreement between the calculated shortwave albedo and observations on 30 June 2012, 16:45 UTC. The assigned factor was 0.88.

Comparison of snow temperature profiles in the NSL at the SIGMA-A site between snow pit observations (OBS) and profiles simulated with the SMAP model at 11:30 UTC on 8 July 2012. Relative snow depth denotes the height of the NSL above the thick bottom layer of ice. The snow depths simulated by the SMAP model were adjusted to the measurements.

[Figure omitted. See PDF]

Model evaluation using the data at the SIGMA-A site

In this section we evaluate the SMAP model using the data measured at the SIGMA-A site during the IOP to adapt it to the GrIS snowpack. The following parameters were validated quantitatively: snow temperature profile, snow surface temperature, surface snow grain size, and shortwave albedo. When a measured snow temperature profile was compared against simulation results, the depth of the NSL simulated by the model (the SMAP model tended to underestimate the NSL's depth by $-$2.0 cm compared to the snow pit measurements during the IOP) was adjusted to the measured depth as a post-process in which model simulated internal properties were not modified at all. The measured snow surface temperature $T_{{\text{s}}_{\mathrm{0}}}$ was calculated from the following relationship between the observed downward and upward long-wave radiant fluxes ($L^{\downarrow }$ and $L^{\uparrow }$): $$L^{\uparrow }={\mathit{\epsilon }}_{\text{s}}\mathit{\sigma }{T}_{{\text{s}}_{\mathrm{0}}}^{\mathrm{4}}+\left(\mathrm{1},-,{\mathit{\epsilon }}_{\text{s}}\right)L^{\downarrow },$$ where $\mathit{\sigma }$ is the Stefan–Boltzmann constant. The model performance was assessed in terms of the root mean square error (RMSE) and mean error (ME) (Table 1). In this paper ME is defined as the average difference between simulated and observed values.

Comparison of SMAP simulation results with in situ measurements. RMSE and ME are the root mean square error and mean error (the average of the difference between simulated values and observed values), respectively. Figures in parentheses indicate scores when the observed surface temperature was negative (i.e., dry snow conditions).

Parameters	RMSE	ME
Snow temperature profile (${}^{\circ }$C)	0.60 (0.81)	$-$0.16 ($-$0.24)
Snow surface temperature (${}^{\circ }$C)	0.58 (0.94)	0.25 (0.68)
Surface snow grain size (mm)	0.21	0.17
Shortwave albedo	0.022	0.008

Snow surface temperature at the SIGMA-A site during the IOP observed with the AWS and simulated by the SMAP model. Observed snow surface temperature was calculated from observed downward and upward long-wave radiant fluxes (Sect. 4).

[Figure omitted. See PDF]

Snow temperature profile

We first examined whether the SMAP model could reproduce the internal physical states of the NSL in terms of the snow temperature profile. Model performance during the IOP (Table 1) indicated that the model simulated the temperatures reasonably. An example profile comparison on 8 July showed that the top 40 cm of the observed snowpack was wet (Fig. 6), and the condition was reproduced by the SMAP model. Even when such observed wet snow conditions were excluded from the statistical assessment of model performance, the order of magnitude of the RMSE and ME did not change (Table 1).

Snow surface temperature

Snow surface temperature affects all energy flux components of the SEB except for the net shortwave radiant flux. A comparison of observed and simulated snow surface temperatures during the IOP (Fig. 7) as well as the obtained MEs (Table 1) indicated that the SMAP model tended to overestimate snow surface temperature; however, the RMSE obtained in this study (0.58 ${}^{\circ }$C) was an improvement compared to the previous study by Niwano et al. (2012) at Sapporo during the winters of 2007–2009 (2.45 and 2.3 ${}^{\circ }$C for the winters of 2007–2008 and 2008–2009, respectively). In this case as well, scores under dry snow surface conditions were still reasonable (RMSE $=$ 0.94 ${}^{\circ }$C and ME $=$ 0.68 ${}^{\circ }$C). These foregoing statistics are almost the same order of magnitude as the analogous statistics from previous, detailed snow modeling studies performed in polar regions (Kuipers Munneke et al., 2009, 2012; Brun et al., 2011).

Surface snow grain size

For the SMAP model, snow grain size is one of the key parameters to be simulated accurately, because the SMAP model calculates the snow albedo and solar heating profile in the snowpack using the PBSAM, for which snow grain size is an input parameter. Figure 8a depicts the simulated snow grain size profile in the NSL, and Fig. 8b compares the surface snow grain size determined from in situ manual measurements and GRS against the simulation result. Because the GRS estimates agreed well with in situ measurements, and because we could obtain more data from GRS than from in situ manual measurements, model validation was performed against GRS data in the present study. The simulated surface snow grain size sometimes decreased abruptly (Fig. 8b), although a new snowfall event was not observed during the IOP. The abrupt decrease in grain size was caused by rapid surface melting and subsequent continuous exposure of layers beneath that were simulated by the SMAP model (Fig. 8a). The order of magnitude of the simulated snow grain size was almost the same with the SMAP model and GRS, and the RMSE (0.21 mm) and ME (0.17 mm) (Table 1) were almost the same order of magnitude as the analogous errors reported by Niwano et al. (2012) for seasonal snow simulations at Sapporo, Japan, during the winters of 2007–2009 (RMSE $=$ 0.31 and 0.15 mm, and ME $=$ $-$0.04 and $-$0.02 mm for the winters of 2007–2008 and 2008–2009, respectively). During 3–5 July, when we observed surface hoar with small grains of snow on the surface at the SIGMA-A site (Yamaguchi et al., 2014), discrepancies between observations and simulations stood out (Fig. 8b). The reason for this failure is that the SMAP model could not simulate surface hoar formation during this period. In the SMAP model surface hoar is created when the latent heat flux is positive and the wind speed is less than 3 m s${}^{-\mathrm{1}}$ (Föhn, 2001); however, the simulated $H_{\text{L}}$ from the night of 2 July to the morning of 5 July was continuously negative (more detailed discussion follows in Sect. 6.2).

In the SMAP model, the latent heat flux is calculated on the assumption that the snow surface is saturated (Niwano et al., 2012), an assumption that is widely made by many physical snowpack models. However, as Box and Steffen (2001) have pointed out, this method can detect surface hoar deposition only in cases of extreme temperature inversion, and this limitation leads to underestimation of surface hoar. Because this result suggests that reconsideration of the method of calculating the latent heat flux in the SMAP model is necessary to obtain an accurate SEB, we discuss this topic further in Sect. 6.2.

Temporal evolution of (a) simulated snow grain size profile in the NSL and (b) surface snow grain size from in situ snow pit observations (black diamond patterns), ground-based remote sensing (GRS) using spectral albedo data measured with a spectrometer (purple dots), and simulated by the SMAP model (solid green curve) at the SIGMA-A site during the IOP. Vertical axis in the upper panel (relative snow depth) denotes the height above the thick bottom ice layer in the NSL.

[Figure omitted. See PDF]

Shortwave albedo

Now we focus on the shortwave albedo, which is another important parameter for estimating snow–atmosphere interactions as well as snow surface temperature. Figure 9 compares observed shortwave albedos with albedos simulated with the SMAP model during the IOP. The comparison shows that the SMAP model successfully reproduced the observed diurnal variations. The RMSE and ME in Table 1 support this assessment. Furthermore, the SMAP model simulated observed spikes when snow albedo rapidly increased, especially around 4 July and after 10 July. The cause of these spikes is the fact that the snow albedo is generally higher under cloudy-sky than under clear-sky conditions (Liljequist, 1956; Yamanouchi, 1983; Aoki et al., 1999). These results imply that the component of the PBSAM driven by the observed near-infrared and UV–visible fractions of the downward shortwave radiant flux, as well as the diffuse fractions calculated by the procedure described in Sect. 3.1, played an important role in improving the precision of the SMAP model simulations. Therefore, once the SMAP model or the PBSAM are coupled with atmospheric models, it is necessary for such host atmospheric models to simulate the presence or absence of cloud realistically. King et al. (2015) also argued that efforts to improve model simulations of surface energy balance and melt in the polar region should concentrate initially on reducing biases in modeled shortwave and long-wave radiant fluxes which are caused by deficiencies in the representation of cloud properties.

Shortwave albedos at the SIGMA-A site during the IOP observed with AWS and simulated by the SMAP model.

[Figure omitted. See PDF]

In Sect. 4.3, however, we found clear discrepancies between observed and simulated surface snow grain size, especially during 3–5 July, when simulated shortwave albedo agreed well with observations. This paradox can be explained as follows: as demonstrated by Yamanouchi (1983), the difference between the downward near-infrared radiant fluxes under clear-sky and cloudy-sky conditions is larger than the analogous difference in the downward visible radiant fluxes. The explanation is that most of the absorption of solar radiation by clouds occurs in the near-infrared region, and the difference in multiple reflection effects due to the spectral differences of surface albedo magnifies the impact of this preferential near-infrared absorption. Figure 2c depicts the observed near-infrared fraction of the downward shortwave radiant flux. Actually, the near-infrared fraction was depleted during the cloudy-sky conditions observed on 3, 4, 6, 7, and 8 July and after 10 July (Sect. 2.1). In the PBSAM, the shortwave albedo is calculated from a weighted summation of visible albedo and near-infrared albedo. The weights for these albedos are the visible and near-infrared fractions of the downward shortwave radiant flux (Aoki et al., 2011). A decrease of the near-infrared fraction therefore increases the influence of the visible albedo on the calculated shortwave albedo, and the influence of the near-infrared albedo, which is mainly affected by snow grain size (Wiscombe and Warren, 1980), decreases. The simulated shortwave albedo therefore agreed with observations, even during 3–5 July, when the SMAP model could not reproduce surface hoar and the associated small near-surface snow grain size.

SEB during the IOP at SIGMA-A

Now we look into the SEB during the IOP at the SIGMA-A site to elucidate the physical conditions of both the snowpack and the atmosphere that led to the remarkable melting around 12 July, when the record near-surface melt occurred over more than 95 % of the entire surface of the GrIS (Nghiem et al., 2012). The SEB equation (Armstrong and Brun, 2008) can be written as follows: $$S_{\text{net}}+L_{\text{net}}+H_{\text{S}}+H_{\text{L}}+H_{\text{R}}+H_G=Q_{\text{net}},$$ where $S_{\text{net}}$ is the net shortwave radiant flux, $L_{\text{net}}$ is the net long-wave radiant flux, $H_{\text{S}}$ is the sensible heat flux, $H_{\text{L}}$ is the latent heat flux, $H_{\text{R}}$ is the heat flux associated with rainfall calculated as a function of rainfall rate and a difference in rain temperature (wet bulb temperature is assumed) and surface temperature (Niwano et al., 2012), $H_G$ is the subsurface conductive heat flux, and $Q_{\text{net}}$ is the net energy flux at the snow surface. As already mentioned, these fluxes are defined to be positive when they are directed into the snow surface. The snow surface is heated when the flux is positive and surface melt occurs if the surface temperature is 0 ${}^{\circ }$C, whereas it is cooled if the flux is negative. In this study, $S_{\text{net}}$ and $L_{\text{net}}$ were calculated from in situ measurements. Other fluxes on the left-hand side of Eq. (6) were estimated as a function of measured snow surface temperature by using the SMAP model. Figure 10 shows the temporal changes of the 1 h averaged simulated SEB during the IOP. $S_{\text{net}}$ remained positive throughout the IOP. $L_{\text{net}}$ was negative for much of the time prior to 9 July but was positive most of the time after 10 July. $H_{\text{S}}$ was close to 0 W m${}^{-\mathrm{2}}$ until 9 July, but it gradually increased after 10 July and sometimes reached about 50 W m${}^{-\mathrm{2}}$. The other turbulent heat flux, $H_{\text{L}}$, was negative most of the time during the first half of the IOP, but it became generally positive after 10 July as well. After 10 July there was heavy rainfall frequently, but its impact on the SEB was quite small ($H_{\text{R}}$ was less than 10 W m${}^{-\mathrm{2}}$ most of the time, although it sometimes reached more than 15 W m${}^{-\mathrm{2}}$). Finally, $H_G$ showed clear diurnal variation: it heated the surface especially during the night time, while it was almost 0 W m${}^{-\mathrm{2}}$ in the daytime as a result of isothermal profile in the near-surface snowpack caused by meltwater percolation. As a result, until 9 July the total surface energy flux, $Q_{\text{net}}$, clearly varied diurnally, being negative during the night and positive during the day. However, after 10 July it remained positive at all times.

One-hour averaged surface energy balance at the SIGMA-A site during the IOP. The upper panel (a) shows net shortwave radiant flux ($S_{\text{net}}$), net long-wave radiant flux ($L_{\text{net}}$), sensible heat flux ($H_{\text{S}}$), latent heat flux ($H_{\text{L}}$), heat supply by rainfall ($H_{\text{R}}$), and subsurface conductive heat flux ($H_G$). These fluxes are defined to be positive when they are directed into the snow surface. The lower panel (b) denotes net energy flux ($Q_{\text{net}}$). The snow surface is heated when $Q_{\text{net}}$ is positive, whereas it is cooled if $Q_{\text{net}}$ is negative.

[Figure omitted. See PDF]

Because there was a clear contrast in the calculated SEB characteristics during the period from 30 June to 9 July (“Period-1”) and the time interval 10–14 July (“Period-2”), we compared the average SEB components between Period-1 and Period-2 (Fig. 11) to characterize the SEB at the SIGMA-A site around 12 July, when continuous melting was observed. Figure 11 reveals dramatic modulations of both $S_{\text{net}}$ and $L_{\text{net}}$, and changes in $H_{\text{S}}$ and $H_{\text{L}}$ were also remarkable. Quantitatively, there were significant positive increments in $H_{\text{S}}$ (4.0 to 15.3 W m${}^{-\mathrm{2}}$; $+$11.3 W m${}^{-\mathrm{2}}$), $H_{\text{L}}$ ($-$18.0 to 7.2 W m${}^{-\mathrm{2}}$; $+$25.2 W m${}^{-\mathrm{2}}$), and $L_{\text{net}}$ ($-$54.9 to 2.1 W m${}^{-\mathrm{2}}$; $+$57.1 W m${}^{-\mathrm{2}}$), the total positive increment being 96.7 W m${}^{-\mathrm{2}}$. There was a remarkable decrease only in $S_{\text{net}}$ (83.0 to 40.7 W m${}^{-\mathrm{2}}$; $-$42.3 W m${}^{-\mathrm{2}}$). As a result, $Q_{\text{net}}$ increased by 49.1 W m${}^{-\mathrm{2}}$ (20.1 to 69.2 W m${}^{-\mathrm{2}}$). These values, calculated during Period-1, are almost equal to the surface fluxes from June to August averaged during the summers of 2000–2011 over the GrIS accumulation area based on the MAR regional climate model (Fettweis et al., 2011) and MODIS data presented by Box et al. (2012). However, the SEB characteristics during Period-2 (signs and orders) resemble those obtained at Langfjordjøkelen, Norway (Giesen et al., 2014).

Comparison of average surface energy balance components ($S_{\text{net}}$ is net shortwave radiant flux, $L_{\text{net}}$ is net long-wave radiant flux, $H_{\text{S}}$ is sensible heat flux, $H_{\text{L}}$ is latent heat flux, $H_{\text{R}}$ is heat supply by rainfall, $H_G$ is subsurface conductive heat flux, and $Q_{\text{net}}$ is net energy flux) between Period-1 and Period-2 at the SIGMA-A site.

[Figure omitted. See PDF]

The decrease of $S_{\text{net}}$ was due to both a reduction of downward shortwave radiant flux (352.7 to 203.9 W m$-\mathrm{2}$) and a slight increase in shortwave albedo from 0.775 to 0.810. The latter value is as high as the 2000–2011 summer (June to August) average albedo over the GrIS accumulation area (0.809) reported by Box et al. (2012) and the cause of increase was mainly the appearance of clouds (discussed in Sect. 4.4). Modulation of $H_{\text{S}}$ was caused mainly by an increase of temperature, which is clearly apparent in Fig. 2a. The cause of the increase in $H_{\text{L}}$ was basically the same: increases in relative humidity and air temperature. The latter reached more than 0 ${}^{\circ }$C. The abrupt transition of radiative properties is understandable from the perspective of cloud radiative forcing (defined as the difference in net surface radiant fluxes under cloudy-sky and clear-sky conditions) in the polar region: at the snow surface shortwave cloud radiative forcing is negative and long-wave cloud radiative forcing is positive (Aoki and Yamanouchi, 1992). We thus conclude that the appearance of low-level clouds (Sect. 2.1) accompanied by a remarkable increase of temperature (Fig. 2a) played an important role and induced surface heating during Period-2. The heating occurred via a large increase in $L_{\text{net}}$ which was able to substantially compensate for the reduction of $S_{\text{net}}$. Bennartz et al. (2013) have investigated the effect of clouds on the temporal development of surface temperature by performing numerical sensitivity tests with a prognostic surface energy balance model and by parameterizing downward radiant fluxes at the surface. They have argued that low-level liquid clouds played a critical role in the enhancement of surface melting at Summit, Greenland, by increasing near-surface temperature through their radiative effects. The present results, which highlight the importance of low-level clouds, are consistent with their conclusion.

Discussion

In this section, we investigate uncertainties in the calculated SEB that can be induced by model settings and calculation methods of the SMAP model. The purpose of this work is to confirm the validity of SEB characteristics presented in Sect. 5. Because $S_{\text{net}}$ and $L_{\text{net}}$ are directly given from in-situ measurements (mentioned in Sect. 5), only uncertainties related to turbulent heat fluxes are discussed here.

Effects of model settings on the calculated SEB

As mentioned in Sect. 3.2, we assumed constant values for $z_{\mathrm{0}}$ (0.12 and 1.3 mm for snow before and after melting, respectively) following Greuell and Konzelmann (1994) in this study. However, as Braithwaite (1995) pointed out, the uncertainty in $z_{\mathrm{0}}$ causes significant error in turbulent heat flux calculations. Therefore, firstly we performed numerical sensitivity tests based on the SEB calculation performed in Sect. 5, where the value of $z_{\mathrm{0}}$ was perturbed between the possible maximum and minimum values. These values were determined as follows: from the compilation by Brock et al. (2006) the maximum value for $z_{\mathrm{0}}$ of the GrIS snow surface was found to be around 10 mm, while we assumed the minimum value for $z_{\mathrm{0}}$ of the GrIS snow surface to be around 0.01 mm based on the results by Smeets and van den Broeke (2008). Compared to the original run ($H_{\text{S}}$ and $H_{\text{L}}$ during IOP were estimated to be 7.6 and $-$10.0 W m${}^{-\mathrm{2}}$, respectively), results from sensitivity studies revealed that absolute values of the turbulent heat fluxes were increased with the condition that $z_{\mathrm{0}}$ was set to be 10 mm (mean differences against the original SEB calculation were 0.8 W m${}^{-\mathrm{2}}$ and $-$1.4 W m${}^{-\mathrm{2}}$ for $H_{\text{S}}$ and $H_{\text{L}}$, respectively; the same tendency was also reported by van As, 2011), while they were reduced if $z_{\mathrm{0}}$ was set to be 0.01 mm (mean differences against the original SEB calculation were $-$2.9 and 3.8 W m${}^{-\mathrm{2}}$ for $H_{\text{S}}$ and $H_{\text{L}}$, respectively). However, these uncertainties still do not affect characteristics of temporal changes in SEB during IOP at the SIGMA-A site discussed in Sect. 5 so much.

Secondly, we investigated effects of ${\mathit{\epsilon }}_{\text{s}}$ introduced in Sect. 3.3 on SEB calculations. In this sensitivity test, ${\mathit{\epsilon }}_{\text{s}}$ was set to be 1.0 and surface temperature (to be input to the SMAP model) was calculated only from observed $L^{\uparrow }$. The result indicated that mean differences of turbulent heat fluxes against the original SEB calculation were 1.1 and 0.9 W m${}^{-\mathrm{2}}$ for $H_{\text{S}}$ and $H_{\text{L}}$, respectively. This result implies that the sensitivity of SEB calculation to the choice of ${\mathit{\epsilon }}_{\text{s}}$ is small, and SEB characteristics during IOP at the SIGMA-A site presented in Sect. 5 is still valid. van As (2011) also performed this type of sensitivity test and demonstrated ${\mathit{\epsilon }}_{\text{s}}$'s small impact on the SEB calculation.

Uncertainties caused by the latent heat flux calculation method

As discussed in Sect. 4.3, calculation of the latent heat flux $H_{\text{L}}$ based on the assumption that air at the surface is saturated with respect to ice at the snow surface temperature might lead to underestimation of surface hoar deposition, which implies that internally diagnosed $H_{\text{L}}$ in the SMAP model simulation presented in Sect. 4 can be underestimated. Here we refer to this method as 1LM, in accordance with Box and Steffen (2001). According to Box and Steffen (2001), employing two-level atmospheric data to calculate $H_{\text{L}}$ can solve this problem, an approach we designate as the 2LM method. As noted in Sect. 2.1, the SIGMA-A AWS measured temperature, relative humidity, and wind speed at heights of 3 and 6 m above the surface. In this section we calculate the latent heat flux using the AWS data at these two heights, and we investigate whether the SMAP model can simulate surface hoar formation with small snow grain sizes by using the 2LM method.

The original formulation of $H_{\text{L}}$ employed by the SMAP model (1LM) is based on the bulk method and is expressed as follows (Niwano et al., 2012): $$H_{\text{L}}=\frac{{\mathit{\rho }}_{\text{a}}{L}_v{\mathit{\kappa }}^{\mathrm{2}}{u}_{\mathrm{1}}\left(q_{\mathrm{1}},-,q_{{\text{s}}_{\mathrm{0}}}\right)}{\left[\ln ,\left(\frac{z_{\mathrm{1}}}{z_{\mathrm{0}}}\right),-,{\mathrm{\Psi }}_{\text{M}},\left(\frac{z_{\mathrm{1}}}{L}\right)\right]\left[\ln ,\left(\frac{z_{\mathrm{1}}}{z_Q}\right),-,{\mathrm{\Psi }}_{\text{H}},\left(\frac{z_{\mathrm{1}}}{L}\right)\right]},$$ where ${\mathit{\rho }}_{\text{a}}$ is the density of air, $L_v$ is the latent heat of sublimation or evaporation, $\mathit{\kappa }$ is the von Kármán constant, $u_{\mathrm{1}}$ is the wind speed at a measurement height $z_{\mathrm{1}}$ (in this study, $z_{\mathrm{1}}$ is 3 m), $q_{\mathrm{1}}$ and $q_{{\text{s}}_{\mathrm{0}}}$ are the specific humidity at $z_{\mathrm{1}}$ and the snow surface, respectively, ${\mathrm{\Psi }}_{\text{M}}$ and ${\mathrm{\Psi }}_{\text{H}}$ are profile functions for momentum and heat, respectively, $z_{\mathrm{0}}$, and $z_Q$ are roughness lengths for momentum and moisture, respectively, and $L$ is the Obukhov length. When we now use atmospheric data at two heights ($z_{\mathrm{1}}$ and $z_{\mathrm{2}}$; $z_{\mathrm{2}}$ is 6 m in this study), the formulation of the bulk method with the 2LM model can be rewritten as follows: $$H_{\text{L}}=\frac{{\mathit{\rho }}_{\text{a}}{L}_v{\mathit{\kappa }}^{\mathrm{2}}\left(u_{\mathrm{2}},-,u_{\mathrm{1}}\right)\left(q_{\mathrm{2}},-,q_{\mathrm{1}}\right)}{\left[\ln ,\left(\frac{z_{\mathrm{2}}}{z_{\mathrm{1}}}\right),-,{\mathrm{\Psi }}_{\text{M}},\left(\frac{z_{\mathrm{2}}}{L}\right),+,{\mathrm{\Psi }}_{\text{M}},\left(\frac{z_{\mathrm{1}}}{L}\right)\right]\left[\ln ,\left(\frac{z_{\mathrm{2}}}{z_{\mathrm{1}}}\right),-,{\mathrm{\Psi }}_{\text{H}},\left(\frac{z_{\mathrm{2}}}{L}\right),+,{\mathrm{\Psi }}_{\text{H}},\left(\frac{z_{\mathrm{1}}}{L}\right)\right]},$$ where $u_{\mathrm{2}}$ and $q_{\mathrm{2}}$ are the wind speed and specific humidity at $z_{\mathrm{2}}$, respectively. Other parameters are calculated by the same method used by Niwano et al. (2012). The choice of ${\mathrm{\Psi }}_{\text{M}}$ and ${\mathrm{\Psi }}_{\text{H}}$ depends on stability conditions in the atmospheric boundary layer. When the atmosphere is stable, the SMAP model assumes that ${\mathrm{\Psi }}_{\text{M}}$ $=$ ${\mathrm{\Psi }}_{\text{H}}$ and calculates the profile functions according to Holtslag and De Bruin (1988), whereas the SMAP model carries out the calculations with functions determined by Paulson (1970) under unstable conditions.

(a) Temporal evolution of 1 h averaged latent heat fluxes calculated from atmospheric measurements at 6 and 3 m above the surface (OBS_2LM) and measurements at 3 m above the surface (OBS_1LM), as well as the SMAP simulation forced by one-level (3 m above the surface) atmospheric measurements (SMAP_1LM). (b) Comparison of surface snow grain size at the SIGMA-A site during the IOP with surface snow grain size from ground-based remote sensing (GRS) using spectral albedo data measured with a spectrometer and those from the SMAP simulations driven by two-level (6 and 3 m above the surface) atmospheric measurements (SMAP_2LM) or one-level (3 m above the surface) atmospheric measurements (SMAP_1LM).

[Figure omitted. See PDF]

Figure 12a shows the temporal changes of the 1 h averaged latent heat fluxes calculated from two-level atmospheric measurements using the 2LM (OBS_2LM) compared to the result presented in Sect. 5 (OBS_1LM) as well as original SMAP simulation result (SMAP_1LM) presented in Sect. 4. Obviously, the result indicated that the latent heat flux from OBS_1LM was almost compatible with the heat flux from SMAP_1LM, implying that the OBS_1LM is also likely to underestimate the latent heat flux. Comparison between OBS_2LM and OBS_1LM shows that the former obviously tends to be higher than the latter, and sometimes the signs of the fluxes are different. According to previous studies (Box and Steffen, 2001; Cullen et al., 2014), the sign and magnitude of the latent heat fluxes from the 1LM and 2LM methods agree reasonably at low elevations on the GrIS, whereas they often differ from each other at the higher elevations. Measurements conducted previously in the northwest GrIS (the Humboldt and GITS sites) showed that the net annual sublimation from the 1LM and 2LM methods did not agree sufficiently at both sites (Box and Steffen, 2001). In the former site, the sign was contrasting, while the magnitude was remarkably different at the latter site.

We next performed another numerical simulation with the SMAP model (SMAP_2LM), where shortwave albedo and snow surface temperature were simulated with the same manner as performed in Sect. 4. In this case, the latent heat flux determined directly from the 2LM method was used only for the assessment of surface hoar formation, while that from the 1LM method was employed for a boundary condition of the governing one-dimensional diffusion equation in the SMAP model (Niwano et al., 2012). Figure 12b compares the surface snow grain size from the SMAP_2LM case with the surface snow grain size from the SMAP_1LM case and GRS. A difference between results of the SMAP_2LM and SMAP_1LM cases is apparent on 4 July (Fig. 12b), when snow grain size is lower for the SMAP_2LM case than for the SMAP_1LM case. The low grain size estimated with SMAP_2LM ($<$ 0.4 mm) can be attributed to modeled surface hoar formation. Actually, the latent heat flux simulated with OBS_2LM on 4 July was positive (Fig. 12a), and the winds at 3 m above the surface around 4 July were calm (Fig. 2b). These conditions are suitable for modeled surface hoar formation, as explained in Sect. 4.3. Although SMAP_2LM succeeded in reproducing the surface hoar detected at the SIGMA-A site during the IOP, which implies that the 2LM method calculated latent heat flux might be more probable than that from the 1LM method, there were still obvious discrepancies with GRS results around 2–5 July (Fig. 12b). To resolve this problem, it might be necessary to consider the possibility of a new snow grain size for surface hoar, which is now calculated as a function of air temperature in the same manner as new snowfall (Niwano et al., 2012).

According to Box and Steffen (2001), the uncertainty of the 2LM method increases as the temperature, humidity, and wind speed differences between two measurement heights decrease. This motivated us to investigate the significance of the latent heat flux calculated by the 2LM method during IOP. In this inquiry, gradients (positive downward) of wind speed, temperature, and vapor pressure between the surface and the lower measurement height and between the lower and upper measurement height were investigated at first. Averaged gradients between the surface and the lower measurement height during the IOP were 1.6 s${}^{-\mathrm{1}}$, 0.3 ${}^{\circ }$C m${}^{-\mathrm{1}}$, and $-$0.15 hPa m${}^{-\mathrm{1}}$, respectively. The value for vapor pressure is very close to that obtained at Summit during 2000–2002 as reported by Cullen et al. (2014). However, averaged gradients between the lower and upper measurement heights were nearly 0, except for the case of wind speed, for which it was 0.2 s${}^{-\mathrm{1}}$. Focusing on the period from 00:10 to 00:20 UTC on 4 July when SMAP_2LM detected the surface hoar, vapor pressure gradients showed opposite signs: $-$0.13 hPa m${}^{-\mathrm{1}}$ for the 1LM method and 0.01 hPa m${}^{-\mathrm{1}}$ for the 2LM method. Although this result explains the reason why only the 2LM method succeeded in the surface hoar detection, the latter value is still very small. These make it difficult to assess uncertainties of the 2LM method caused by each sensor as expected. In fact, numerical sensitivity studies with perturbed input parameters considering absolute accuracy of temperature, relative humidity, and wind speed sensors ($\pm$0.2 ${}^{\circ }$C, $\pm$2 %, and $\pm$0.3 m s${}^{-\mathrm{1}}$, respectively) in the 2LM calculation modified the picture of calculated turbulent heat fluxes drastically in any calculations. Even when relative differences in the accuracy of two sensors at the lower and upper measurement heights were considered (according to our relative calibration of the instruments performed in advance, air temperature and wind speed sensors at two levels showed no significant difference; however, as for relative humidity, the upper sensor tended to be lower by 1.2 % compared to the lower sensor), differences in calculated latent heat flux with perturbed input parameters were quite large as well. Therefore, we should conclude that underestimation of the latent heat flux calculated with the 1LM method is plausible, although the exact order of underestimation was quite hard to detect during this study period.

Conclusions

In this study, we investigated the record near-surface melt event that occurred over the GrIS during the IOP (30 June to 14 July 2012) from the standpoint of the SEB. We used data measured in situ at the SIGMA-A site, where significant increases of air temperature, relative humidity, and downward long-wave radiant flux, as well as heavy rainfall and abrupt near-surface snowmelt, were observed beginning on 10 July. Although radiation-related components of the SEB could be determined from AWS data, other fluxes were not measured directly. We therefore employed the physical snowpack model SMAP to calculate the $H_{\text{S}}$, $H_{\text{L}}$, and $H_{\text{R}}$. Because this was the first attempt to adapt the SMAP model to a polar region, we carried out a preliminary analysis of various aspects of the performance of the SMAP model. We calculated the snow temperature profile in the NSL, snow surface temperature, surface snow grain size, and shortwave albedo; we compared these calculated values with in situ measurements. In the numerical simulation, the initial conditions of the snow were specified from the snow pit measurements conducted on 30 June 2012, 16:45 UTC. The SMAP model was subsequently driven by 5 min averaged meteorological data, including air pressure, air temperature, relative humidity, wind speed, and downward shortwave, near-infrared, and long-wave radiant fluxes. The SMAP model was also driven by the daily temperature at the bottom of the thick ice layer in the NSL, the mass concentrations of snow impurities obtained every other day, and 5 min averaged accumulated precipitation based on 3-hourly ERA-Interim reanalysis data corrected by in situ bucket measurements.

Validation results revealed that the RMSE for the snow temperature profile and snow surface temperature were reasonable. Regarding surface snow grain size, simulation results were compared against GRS data obtained from spectral albedo measurements post-processed with a GSAF algorithm. Although the RMSE and ME that we obtained were comparable in magnitude to those reported from previous model validation studies performed at Sapporo, Japan, the small snow grain size associated with the surface hoar observed during 3–5 July could not be simulated by the SMAP model. In the SMAP model, surface hoar is created when $H_{\text{L}}$ is positive and the wind speed is less than 3 m s${}^{-\mathrm{1}}$. However, the simulated $H_{\text{L}}$ from the night of 2 July to the morning of 5 July was continuously negative. Despite these discrepancies, the simulated shortwave albedo was in reasonable agreement with observations throughout the IOP (RMSE $=$ 0.023 and ME $=$ 0.008). The cause of the agreement between the simulated shortwave albedo and observations, even during the 3–5 July period when the SMAP model could not reproduce surface hoar and associated small near-surface snow grain size, was a decrease of the near-infrared fraction of the downward shortwave radiant flux caused by the appearance of low-level clouds during the IOP. This change in the shortwave radiant flux increased the influence of visible albedo and in turn decreased the effect of near-infrared albedo on the shortwave albedo. These physical processes are explicitly taken into account by the PBSAM, an important component of the SMAP model, which is an advantage of PBSAM.

Using the SMAP model, we investigated the SEB at the SIGMA-A site during the IOP. Radiation-related components of the SEB were directly given from the AWS measurements, whereas other components were calculated by the SMAP model as a function of measured snow surface temperature. The calculated SEB was clearly different between Period-1 (30 June to 9 July) and Period-2 (10 to 14 July): $L_{\text{net}}$ increased dramatically by $+$57.1 W m${}^{-\mathrm{2}}$ ($H_{\text{S}}$ and $H_{\text{L}}$ also increased by $+$11.3 and $+$25.2 W m${}^{-\mathrm{2}}$, respectively) after 10 July, whereas $S_{\text{net}}$ decreased significantly by $-$42.3 W m${}^{-\mathrm{2}}$. Consequently, $Q_{\text{net}}$ clearly varied diurnally (negative during the night and positive during the day) until 9 July. However, the fact that it remained continuously positive after 10 July explains the continuous melt event observed at the SIGMA-A site. We discussed the reason for this remarkable transition of radiative properties, and we concluded that it was caused by the appearance of low-level clouds accompanied by a significant temperature increase. The result was surface heating during Period-2 via the process of cloud radiative forcing in the polar region.

In order to confirm the validity of SEB characteristics during IOP, additional error analyses were conducted. During this process, it turned out that the sign of latent heat fluxes from the 1LM and 2LM methods differed especially when the surface hoar was observed (around 4 July). The former showed negative, while the latter turned positive and designated the surface hoar formation. Therefore, the latent heat flux calculated by the 2LM method seemed to be plausible; however, uncertainty involved in the 2LM method was so large that we could not confirm its significance.

M. Niwano developed the SMAP model code and performed the numerical calculations. T Aoki designed the SIGMA AWS system. T. Aoki, S. Matoba, S. Yamaguchi, T. Tanikawa, M. Niwano, K. Kuchiki, and H. Motoyama installed the SIGMA-A AWS on the GrIS and constructed data-acquisition system. S. Matoba and S. Yamaguchi performed snow pit measurements. T. Aoki and K. Kuchiki processed data of mass concentrations of snow impurities. M. Niwano, T. Aoki, and T. Tanikawa conducted spectral albedo measurements. K. Kuchiki retrieved GRS from the spectral albedo data. M. Niwano prepared the manuscript with contributions from all co-authors.

Acknowledgements

We sincerely thank Tetsuhide Yamasaki for logistical and field support of our expedition and Sakiko Daorana for her help during our stay at Qaanaaq, Greenland. We thank Climatec Inc. (Japan) for manufacturing the AWS installed at the SIGMA-A site and Masae Igosaki for supporting the laboratory measurements of snow impurities. We would like to thank two reviewers, John King and an anonymous referee, and the scientific editor of this paper, Michiel van den Broeke, for their helpful comments and suggestions. This study was supported in part by (1) the Japan Society for the Promotion of Science (JSPS) grant-in-aid for scientific research (S) number 23221004, (2) the Global Change Observation Mission – Climate (GCOM-C)/the Second-generation GLobal Imager (SGLI) Mission, the Japan Aerospace Exploration Agency (JAXA), (3) the Experimental Research Fund for Global Environment Conservation, Ministry of the Environment of Japan, and (4) the Grant for Joint Research Program, Institute of Low Temperature Science, Hokkaido University. The map showing the location of site SIGMA-A (Fig. 1) was created by NunaGIS (http://en.nunagis.gl/) operated by Asiaq, Greenland Survey. Edited by: M. van den Broeke