Forcing and Responses of the Surface Energy Budget at Summit, Greenland

Greenland ice sheet surface temperatures are controlled by an exchange of energy at the surface, which includes radiative, turbulent and ground heat fluxes. Data collected by multiple projects are leveraged to calculate estimates of all surface energy budget (SEB) terms at Summit, Greenland for the full annual cycle from July 2013 June 2014 and extend to longer periods for the radiative and turbulent SEB terms. Radiative fluxes are measured directly by a suite of broadband radiometers. Turbulent 5 sensible heat flux is primarily estimated via the bulk aerodynamic method, and the turbulent latent heat flux is calculated via a two level approach using measurements at 10 and 2 m. The subsurface heat flux is calculated using a string of thermistors buried in the snow pack. Extensive quality-control data processing produced a data set in which all terms of the SEB are present 75% of the full annual cycle, despite the harsh conditions. By including a storage term for a near surface layer, the SEB is balanced in this data set to within the aggregated uncertainties for the individual terms. November and August case studies illustrate that 10 surface radiative forcing is driven by synoptically forced cloud characteristics, especially by low-level, liquid-bearing clouds. The annual cycle and seasonal diurnal cycles of all SEB components indicate that the non-radiative terms are anti-correlated to changes in the total radiative flux, and are hence responding to cloud radiative forcing. Generally, the non-radiative SEB terms and the upwelling longwave radiation component compensate for changes in downwelling radiation, although exact partitioning of energy in the response terms varies with season and near-surface characteristics such as stability and moisture 15 availability. Substantial surface warming from low-level clouds typically leads to a change from a very stable to a weakly stable near-surface regime with no solar radiation or from a weakly stable to neutral/unstable regime with solar radiation. Relationships between forcing terms and responding surface fluxes show that the upwelling longwave radiation produces 5575% (40-50%) of the total response in the winter (summer) and the non-radiative terms compensate for the remaining change in the combined downwelling longwave and net shortwave radiation. Because melt conditions are rarely reached at Summit, 20 these relationships are documented for conditions of surface-temperature below 0◦C, with and without solar radiation. This is the first time that forcing and response term relationships have been investigated in detail for the Greenland SEB. These results should both advance understanding of process relationships over the Greenland icecap and be useful for model validation. 1 The Cryosphere Discuss., doi:10.5194/tc-2016-206, 2016 Manuscript under review for journal The Cryosphere Published: 23 September 2016 c © Author(s) 2016. CC-BY 3.0 License.


Introduction
The exchange of energy at the Greenland Ice Sheet (GIS) surface must be thoroughly characterized to fully understand the processes that govern surface temperature variability, which is important in monitoring and modeling ice sheet mass balance (Box, 2013). Observations suggest near-surface temperatures are increasing; the GIS is trending toward greater spatial melt extent (McGrath et al., 2013) with increased melt runoff due to atmospheric warming (Hanna et al., 2008). The amalgamated 5 freshwater runoff, in combination with ice discharge, determines how this major reservoir of Northern hemisphere ice affects freshwater input into the North Atlantic and Arctic Oceans, and subsequently, global ocean circulation and sea level rise.
Surface melt processes currently account for approximately half of the total mass loss of the entire GIS (van den , and during prolonged periods of elevated surface temperatures this proportion is even greater (Smith et al., 2014). The melt process occurs in two steps. First, energy flux to the surface is used to increase the surface temperature. Then, after the 10 melting point is reached, excess net surface energy flux is used to convert ice into liquid water. As an increasing area of the interior GIS approaches the melting point of snow in summer, spatial and temporal variations of the net surface energy flux are paramount in determining when the melting point is reached, over what spatial area this occurs, and the amount and rate of melt after this threshold is reached.
The surface energy budget (SEB) is a balance of radiative, turbulent and ground heat fluxes, which are coupled through 15 a variety of processes. Once the surface temperature reaches the melting point of snow, additional energy goes toward melt, limiting the surface temperature to 0 • C. In the absence of phase change, however, a change in one of the SEB terms must be balanced by a change in another term or combination of terms. Importantly, the surface temperature is related to multiple SEB terms including upwelling longwave radiation, turbulent sensible heat, and ground heat fluxes. Over time scales long enough for the surface temperature to adjust, closure of the SEB is achieved and all of the energy exchange at the surface is accounted for. 20 Because of the high emissivity (and hence high longwave absorptivity) of the snowpack, the surface is able to adjust relatively quickly to longwave influences (e.g., whether that is a warm cloud or a cold, clear sky). In contrast to its efficient ability to absorb longwave radiation, the GIS has a high shortwave albedo and reflects much of the incoming solar radiation. Liquidbearing clouds are frequent above the GIS during summer (Shupe et al., 2013b) and have strong implications for increasing melt extent (Bennartz et al., 2013) and meltwater runoff (Van Tricht et al., 2016). In fact, clouds act to radiatively warm the 25 central GIS throughout the year (Miller et al., 2015;Van Tricht et al., 2016), more than would occur via solar radiation acting alone, as a result of the year-round high surface albedo. Thus, the primary radiative influences on raising surface temperatures in this region are the solar zenith angle and occurrence of clouds.
A change in the downwelling radiative flux caused by clouds and/or solar radiation will induce a response of the atmospheric boundary layer and surface. Boundary layer depth and stability are influential for exchange processes, such as sublimation 30 fluxes, which modulate accumulation (Berkelhammer et al., 2016). Miller et al. (2013) shows a degradation of the surfacebased temperature inversion in the presence of liquid-bearing clouds, which impacts the near-surface stability (Hudson and Brandt, 2005) and thus turbulent mixing. A regional modeling case study by Solomon et al. (2016) indicates also that the response of turbulent and conductive heat fluxes to cloud radiative forcing is important when considering surface-atmosphere interactions. Investigating these responses and interactions throughout the year is paramount in discerning the net effect of liquid-bearing clouds on surface temperatures and, consequently, on sub-surface temperatures and melt processes.
The central GIS is a massive reservoir of snow and ice, responding to energy changes at the surface by conducting heat into or out of the subsurface. Thus, the ice sheet damps the effects of either strong radiative warming or cooling at the atmosphere-snow interface. Warmer subsurface temperatures, resulting from warming of surface temperatures, can change the snow morphology 5 and precondition the surface to have less capacity for removing subsequent heat excesses generated by atmospheric processes (Solomon et al., 2016). Proper atmosphere/ice sheet coupling is important to allow for physically realistic radiational cooling at the surface, in order to minimize surface temperature biases in forecast models (Dutra et al., 2015).
Regional and global climate models are a critical tool for understanding the fate of the GIS and attempt to capture the non-trivial interactions between the atmosphere and the GIS. Early studies parameterized the SEB of the GIS using meteoro-10 logical measurements from summer camps in western Greenland and observations of albedo from satellites (van de Wal and Oerlemans, 1994;Konzelmann et al., 1994). More recently, computationally advanced, fully coupled, climate models project enhanced surface melt as GIS surface temperatures increase under future CO 2 forcing scenarios (Vizcaíno et al., 2014). Yet, these state-of-the-art climate models have surface temperature biases over the GIS, likely due to the under representation of liquid-bearing clouds (Kay et al., 2016). To better understand and represent the important processes that currently hinder 15 models, detailed surface-based observations are valuable.
In western Greenland, detailed measurements of the surface mass balance (van de Wal et al., 2005), surface radiation balance (van den Broeke et al., 2008), and surface energy balance (van den Broeke et al., 2011) have been reported, all of which focus on the ablation zone. In central Greenland, the most sophisticated and comprehensive long-term observations of surface energy budget are made at Summit Station. While a majority of the published literature has focused on the summer season (Cullen 20 and Steffen, 2001;Kuipers Munneke et al., 2009) some studies have targeted SEB annual cycles in 2000(Cullen, 2003 and (Hoch, 2005. In addition, various studies have focused on specific components of the SEB, such as surface latent (Box and Steffen, 2001) or sensible heat (Cohen et al., 2007;Cullen et al., 2007;Drüe and Heinemann, 2007) fluxes.
Annual surface radiation fluxes have been reported at Summit by van den Broeke et al. (2008), Cox et al. (2014), andMiller et al. (2015), as well as longwave flux divergence in the boundary layer by Drüe and Heinemann (2007) and Hoch et al. (2007). 25 Yet, prior to May 2010 there have been limited ground-based measurements of the atmospheric state and cloud properties to complement these temporally sporadic SEB investigations and to support process-based understanding of SEB variability on time scales from minutes to seasons.
This study uses comprehensive ground-based measurements to investigate interactions between the atmosphere and the central GIS throughout the year in order to understand how energy exchange drives temporal variability in surface temperature. 30 Summit Station is currently within the accumulation zone, recording only 2 melt events since 1889 (Nghiem et al., 2012). The lack of melt events provides the opportunity to examine relationships between the various surface energy fluxes in all seasons without the energetic influence of phase change at the surface. Using a unique compliment of data/measurements at 30-minute temporal resolution, we present a pair of case studies to illustrate cloud effects on the balance of energy at the surface and, consequently, the subsurface snow in central Greenland. Next, we characterize the annual and diurnal cycles of the radiative, 35 turbulent and conductive heat fluxes for one year and evaluate SEB closure. Finally, we investigate the seasonal responses of the turbulent heat fluxes, subsurface heat flux, and upwelling longwave flux to changes in downwelling longwave and net shortwave fluxes, establishing process-based energy flux relationships.

Measurements and Methods
Near-surface instrumentation at Summit Station (72 • N 38 • W, 3211 m) is used to characterize the surface energy budget in cen-5 tral Greenland. Net radiative (Q), turbulent sensible (SH), turbulent latent (LH), and total subsurface (G) heat fluxes determine the net surface flux (F s ) according to the equation: (1) The total subsurface heat flux (G) considered here is a combination of the conductive heat flux (C) and heat storage in a nearsurface layer (S), detailed in Section 2.5. Each of these four terms is defined such that a positive value sends energy towards the 10 surface and visa versa. For all measurements described here, a 30-minute time window is used; this time window was chosen to fit the constraints set by eddy covariance calculations for sensible turbulent flux (Section 2.3), but is sufficiently brief to capture both the diurnal cycle and the SEB response to atmospheric variability of interest here.
All SEB terms are estimated for 75.3% of an annual cycle, spanning July 2013 -June 2014, although Q, SH and LH are also measured prior to July 2013. The techniques used to calculate each SEB term, the data availability periods, and associated 15 uncertainties are outlined in the following sub-sections. The estimated uncertainty in each SEB term is summarized in Table 1.
While each component of the SEB has it's own uncertainty, at times the various estimates use the same input and are thus not independent. For example the longwave measurements are used to derive the skin temperature, which is input into both the bulk sensible heat flux and conductive heat flux estimates.  smoothed using a 5-day running window to remove erroneous spikes in the snow depth. Realistic longer term discontinuities due to actual snow events were maintained by limiting the period over which data smoothing occurred. Inexplicably, on 27 May 2014 the sonic ranger reported an abrupt 17.8 cm decrease in the surface height. The near surface thermistor variability indicates that this was unrealistic, hence an offset of -17.8 cm was applied to the thermistor depths thereafter through the end of the study period. The standard deviation over 30-minutes of the 1-minute subsurface temperature data indicates that the 30 variability decays as a function of depth because of a decline in the thermal effects of wind ventilation and direct solar heating due to solar penetration. To minimize the impact of these complicating issues a standard deviation threshold of 0.1 is used to determine that the acceptable minimum depth to use for the shallowest subsurface thermistor is about -20 cm.
The specific humidity at 2 and 10 m, which is needed for deriving LH, is calculated from CIBS relative humidity and temperature measurements in combination with NOAA/GMD temperature and pressure measurements. direct measurements of water vapor mixing ratio are obtained via a Picarro model L2120 spectrometer, which was calibrated using a LiCor LI160 dew point generator (Bailey et al., 2015). The instrument directly samples air moisture content once an hour at multiple levels on the 50-m tower using a constrained inlet system to limit large (> 50 µm) hydrometers from being incorporated into the vapor measurements. Comparing meteorologically derived specific humidity values at approximately 1-2 m and 9-10 m above the surface to the highly accurate Picarro measurements reveals a small bias of +0.065 g kg −1 . The percent 5 error, using the Picarro measurements as truth, at the 2 and 10 m levels are 53% and 30%, respectively.
Investigating the surface flux estimates in combination with active and passive cloud property measurements yields a comprehensive understanding of how clouds affect the GIS energy budget. In addition to the aforementioned radiosondes, ICECAPS also measures the cloud properties via a comprehensive suite of instruments, in operation since May 2010. ICECAPS is described in detail by Shupe et al. (2013b). Liquid water path (LWP) and precipitable water vapor (PWV) are estimated using a 10 physical retrieval via a pair of microwave radiometers, similar to Turner et al. (2007). In a dry environment, such as Summit, it is advantageous to use a total of 3 channels (23.84, 31.40, 90.0 GHz) to increase sensitivity and effectively reduce uncertainty (LWP ≈ 5 g m −2 , PWV ≈ 0.35 mm) (Crewell and Löhnert, 2003). The primary changes to the LWP values estimated in Miller et al. (2015) are an improved liquid-water model (TKC, Turner et al., 2016) and the use of three channels in the retrieval instead of four. By excluding the 150.0 GHz channel, biases in LWP retrievals due to precipitating ice hydrometers will not impact 15 the overall statistical results (Pettersen et al., 2016). Vertically resolved cloud presence is determined by a 35-GHz Millimeter Cloud Radar (MMCR).

Radiative Flux
Four broadband radiation components comprise the net radiation at the surface (Q):

25
Data processing for radiation measurements used here is similar to Miller et al. (2015), including corrections to the LW↓ components based on the net longwave radiation and comparison to colocated broadband radiation measurements operated by NOAA-GMD. An estimated Gaussian longwave radiation measurement uncertainty of 4-5 W m −2 (Gröbner et al., 2014) for the 1-minute data translates to about 1 W m −2 uncertainty in the 30-minute average data. Assuming an emissivity uncertainty of 0.005 a LW-derived surface temperature has an approximate uncertainty of 0.6 • C, which is derived from the radiation 30 measurements thusly: where surface emissivity( ) = 0.985 and σ = is the Stefan-Boltzmann constant. Comparing LW↑ to similar, proximate NOAA/GMD radiation measurements indicates that there is general agreement within the estimated 4-5 W m −2 uncertainty of the longwave radiative components. Yet, for very cold surface temperatures (i.e., < -46 • C) differences between the NOAA/GMD and ETH LW↑ are more pronounced. Hence, a third degree polynomial was used to fit the difference between the ETH and NOAA/GMD LW↑ as a function of the ETH LW↑. A correction factor (y) was applied based on the measured ETH LW↑ (x) value according 5 to the equation: y = −14.99+0.1715x−0.000668x 2 +8.579e−7x 3 , which assumes the more recently calibrated NOAA/GMD pyrgeometers are accurate. After applying the adjustments to LW↑ and LW↓ (Miller et al., 2015) the 1-minute LW data are consistent with a total uncertainty of 4-5 W m −2 .
The surface albedo (SW↑/SW↓) is affected by the solar zenith angle, and for clear-sky days should have a minimum at solar noon. During 2014 an asymmetry in the diurnal cycle is observed in the measured albedo, where the albedo in the 10 morning is up to 10% lower than in the evening. The NOAA/GMD measurements, which are mounted on the same fixed arm, indicate the same issue (possibly a gradual slope to the surface due to snow drifts). There is good agreement between the SW↓ measurements and the total direct plus diffuse SW↓ values when available, suggesting that this issue is unlikely a leveling problem in the SW↓ component. Hence, the SW↑ value is estimated in 2014 using the SW↓ value according to the equation:

25
The net surface flux is influenced by the temperature of the overlying air, i.e. warmer near-surface air will increase the sensible heat transferred to the surface. Direct heat transfer, via conduction, from the atmosphere to the snowpack is only prominent very close to the surface, thus heat is primarily transferred via turbulent eddies. These eddies act to mix the air within the surface layer, reducing the vertical temperature gradient. Estimates of the sensible heat flux are calculated using two independent methods: eddy correlation method and the bulk aerodynamic method. The eddy correlation (EC) method (e.g., Oke, 1987) calculates the covariance between the anomalies in the vertical wind (w ) and temperature (θ ) to determine the turbulent sensible heat flux according to the equation: where the constants are the density (ρ) and heat capacity (c p ) of air. By using direct measurements of windspeed and temperature from a three-dimensional sonic anemometer, an accurate calculation of the heat exchange at ≈ 2 m is obtained.

5
A 30-minute averaging period is a short enough time-window to exclude issues of non-stationarity while still long enough to include low frequency contributions to the turbulent heat flux. Various quality-control (QC) measures are implemented to ensure the data is representative of the entire sensible heat flux during the 30-minute window. QC measures exclude large changes in windspeed or wind direction, upwind contamination by the experimental apparatus, and ± 30% deviations from characteristic -5/3 slope in the inertial subrange (Kaimal et. al. 1972). Applying the QC criteria flags 75% of the available data, 10 spanning September 2011 -June 2014. Thus, for the 85% of this period that either have instrument downtime or where the data are QC flagged, an alternative approach is used.
Due to the limited data set available from the EC method, a bulk aerodynamic method is also used in order to fill in data gaps for the time period June 2011 -June 2014. The bulk transfer method uses Monin-Obukhov similarity theory to estimate turbulent sensible heat flux at the surface: where U is the mean horizontal wind speed at 2 m, T surf is the skin temperature, T 2m is the temperature at 2 m, and C h is the sensible heat transfer coefficient for the 2 m reference height (Persson et al., 2002;Fairall et al., 1996). NOAA/GMD meteorological data is the primary source of the 2 m temperature measurements and data gaps are filled with CIBS temperature data. Cup anemometer measurements fill in data gaps of the sonic anemometer-derived 2 m windspeed measurements. measurements and snow emissivity assumptions, or roughness length values. Sensible heat flux discrepancies could also be due to measurement height differences between the EC and bulk methods. While the bulk method uses the measured surface skin temperature the EC values are measured at 2 m, which could differ from the sensible heat flux directly at the surface under very stable conditions. This suggests that the true SH uncertainty is smaller than estimated here. The covariance u * and bulk u * are well correlated (0.84) with a RMS difference of 0.55 m s −1 and the bulk values are biased low (-0.026 m s −1 ). Changing the 5 velocity roughness length to 4.5e-4 m, which is that determined for snow-covered multi-year sea ice (i.e., Persson et al., 2002) increases the RMS differences for the sensible heat flux by 1.4 W m −2 , suggesting that variability in the roughness of the surface could contribute to error in the bulk parameterization. A majority of the 8.7 W m −2 uncertainty in the bulk estimates is likely due to uncertainties in the skin temperature as estimated from a constant surface emissivity. From June 2011 to June 2014 the bulk estimates are available for 78% of the time period. Thus, filling in EC data gaps with the bulk values vastly 10 improves the temporal coverage of the sensible heat estimates.

Turbulent Latent Heat Flux and Stability
Turbulent eddies also affect the surface energy budget by transferring latent heat toward or away from the surface. Frequently the specific humidity increases with height above the surface, resulting in a transfer of latent energy toward the surface possibly resulting in deposition. The bulk method used by Persson et al. (2002) assumes saturation conditions at the surface, which is 15 not always a valid assumption for dry snow (Albert and McGilvary, 1992). In central Greenland the two-level profile method has been shown to be superior to the bulk method (Box and Steffen, 2001) as it can account for sublimation and deposition to the surface.
The profile method used here is similar to Steffen and DeMaria (1996) such that the latent heat flux is calculated from near-surface horizontal wind (U) and mixing ratio (q) gradients (∆ = value at 10 m -value at 2 m) according to the equation: where ρ is the density of air, L v is the latent heat of vaporization, k is the von Kármán constant (0.4), and z r is the log mean height ( ∆z ln(z2z1 −1 ) ). The stability functions for the transfer of momentum (φ m ) and water vapor (φ e ) are corrections based upon the stability of the boundary layer and will either increase (unstable conditions) or decrease (stable conditions) the surface flux.
A measure of boundary layer stability is attained via calculation of the bulk Richardson number (Ri). The sign of Ri indicates 25 whether mechanical mixing (positive) or buoyancy (negative) is more important in producing turbulence. Ri is dependent on the gradient in virtual potential temperature (∆θ v ),wind speed (∆u) and respective measurement heights (∆z) according to the equation: where g is the acceleration due to gravity (9.81 m s −2 ) andθ v is the average virtual temperature (K) between the two levels. 30 In accordance with Steffen and DeMaria (1996), Ri is used to calculate the stability corrections. Coefficients for relating Ri to the stability factors are obtained from a study conducted in 2000, which used eddy correlation turbulence measurements to Table 2. Stability functions for unstable and stable regimes from Cullen (2003).
obtain the relationships in Table 2 (Cullen, 2003). LH is the data set most susceptible to data gaps because there must be input values of specific humidity, wind speed, and temperature at both the 2 m and 10 m levels. Yet by using the best available meteorological data from NOAA/GMD and/or the 10 CIBS project we estimate LH for 81% of the time period from March 2012 -June 2014. The main driver of uncertainty is the estimation of the mixing ratios with uncertainties of 53% and 30% at 2 m and 10 m, respectively, as compared to the Picarro measurements. The resultant error contribution (60%) to the LH estimate dominates the contribution from uncertainty in the wind speeds.

15
The energy flux from the overlying atmosphere to the subsurface includes direct radiative heating of the snowpack due to solar penetration (Kuipers Munneke et al., 2009), the thermal effects of wind ventilation (Albert and McGilvary, 1992), and conduction. To minimize the complications in calculating sub-surface heat flux caused by the other factors, an estimation of the conductive heat flux (C) at a depth below the solar penetration depth (at least 20 cm) combined with a heat storage (S) in the snow above this level is used to provide an estimation of the total subsurface heat flux (G), such that In this study we calculate the storage heat flux across the uppermost layer and assume the heat flux to the subsurface below is equivalent to C.
The conductive heat flux (C) represents the diffusion of heat between the subsurface and the overlying surface. The effectiveness of the heat transfer is a function of the thermal conductivity of the snow (K) and the vertical temperature gradient 25 (∆T /∆z): The temperature gradient for the uppermost subsurface layer (∆T 01 ) is estimated as the difference between the surface tem- Equation 3) and the temperature measured by the shallowest, sub-surface sensor. To estimate C, at ≈20 cm depth, the conductive heat flux at the two levels bracketing this depth are calculated and averaged, according to the equation: The thermal conductivity of the upper most layers of snow is estimated from average density profile measurements taken layer (≈20 cm), S is calculated by the layer averaged temperature difference (δT ) between chronologically adjacent time steps (δt = 30 minutes), where T 1 is the temperature of the shallowest thermistor at a depth z 1 (simliar to Hoch (2005)); and c ice is the specific heat of ice and ρ is the average density of the layer. The large uncertainty in the skin temperature measurements (0.6 • C) are close to the average temperature change from one time step to the next (0.76 • C) resulting in an The temperature variability at and below the ice sheet surface is important for understanding the flow of heat through this interface and can influence processes like snow compaction and melt. Figure 1  In the spring, fall and winter, surface-based temperature inversions are prevalent (Miller et al., 2013) and the warmest layers of the atmosphere occur between 100-1000 m above ground level as can be seen in Figure 1a. In fact, the minimum temperature

Liquid-bearing cloud with insolation
A case study on 6 August 2013 also illustrates the longwave warming effect of liquid-bearing clouds and investigates the additional influence of shortwave radiation. Similar to the first case study, surface temperature variability is driven by the downwelling radiative flux, which in this case is a combination of longwave and shortwave influences.
MMCR measurements (Figure 3a) indicate a clear-sky scene from 2 to 6 UTC, a low-level cloud from 6 -13.5 UTC, clear-5 sky from 13.5 to 16 UTC, a deep cloud from 18 to 22 UTC, and finally a low-level cloud during the last hour of the case study period. The low-level cloud is mixed phase from 6 to 13.5 UTC and LWP values ranging from 0 -15 g m −2 (Figure 3b) indicate that it is optically thin. LWP values ranging from 0 -20 g m −2 also indicate that the deep cloud later in the day is mixed phase from 18 to 21 UTC although after ≈19 UTC LWP values are low due to competition from falling ice into the mixed phase layer from above. The low-level cloud from 23 to 24 UTC is optically thicker then the previous low-level cloud 10 with LWP ranging from 5 -30 g m −2 .
The presence of the optically thin liquid-bearing cloud (6 -13.5 UTC) produces an approximate increase of 70 W m −2 of LW↓ compared to the preceding clear-sky scene. Over this period shortwave radiation increases the net radiation at the surface by an additional 5 -75 W m −2 . In response, LW↑ radiation increases by 50 W m −2 . The combination of thin liquid-bearing clouds and insolation produces positive net radiation at the surface from 9.5 to 13 UTC (Figure 3c). During the daytime clear- 15 sky period the net radiation is near zero, indicating that shortwave warming is offset by the longwave cooling at the surface. Net radiation again goes positive in the presence of liquid-bearing clouds that occur after 18 UTC. After 18 UTC the net radiation declines as the solar radiation diminishes.
The compensating response of the non-radiative terms to changes in the downwelling radiation, shortwave and/or longwave, are similar to the November case study. The sensible heat flux decreases from 29 W m −2 at 5 UTC to -9 at 12.5 UTC (Fig-20 ure 3d). The fact that the SH is negative during the presence of the liquid-bearing cloud indicates that the surface temperature is warmer than the 2 m temperature, thus the near-surface atmospheric layer is unstable. The conductive heat flux decreases from 9 W m −2 at 5 UTC to 0 W m −2 at 12.5 UTC, indicating the subsurface temperature gradient as been reduced (Figure 3e). The  Over the entire year the SEB residual, or the sum of all the SEB terms, when available (75.3% of the time), is 0.9 W m −2 .
The monthly residuals (top numbers in Figure 4) indicate that there are times of the year when the residuals are larger but there is no apparent seasonality in the combined SEB terms. Generally, the monthly mean residuals could be due to energy imbalances, under sampling, measurement biases, and/or measurement uncertainties. Each monthly residual is below the total SEB uncertainty (excluding the S term) of 12.4 W m −2 .

Diurnal Cycles
The magnitudes of the monthly mean SEB terms are small from May -August (< 10 W m −2 ), yet the diurnal variability peaks during this period, driven largely by the solar cycle. The net radiative flux increases during times of peak insolation (Figure 5a), although the high surface albedo limits the maximum Q to 40 W m −2 . The maximum values of the net radiative flux occur in July, when the sun still rises more than 30 • above the horizon and liquid-bearing clouds are frequent (Figure 6a, b), which act 20 to radiatively warm the surface at Summit Station year round (Miller et al., 2015).
Counteracting the net radiative flux, the sensible heat flux is negative for large sun angles and warms the surface by approximately 20 W m −2 when the sun is below the horizon (Figure 5b). The diurnal variability for this term is largest in summer due to an enhanced diurnal cycle of the near-surface temperature gradient (Miller et al., 2013). The cooling effect of the conductive heat flux (Figure 5c) is most prominent when the sun is above the horizon and is maximized at solar noon. In agreement with 25 the results in Figure 4, more conductive surface cooling occurs in the spring compared to the fall due to the lag in subsurface response, which results in relatively colder subsurface temperatures in the spring. The diurnal variability of the latent heat flux is largest in June ranging from hourly-average values of -33 to 12 W m −2 (Figure 5d) due to an increase in available moisture ( Figure 6c).
Sun angle, and the associated change to the net radiative flux, is a main driver of energy fluxes at the surface ( Figure 5).

Responses to Surface Radiative Forcing
The surface energy budget at Summit Station is largely driven by changes in the downwelling radiation. In general, the LW↑, turbulent, latent, and subsurface heat fluxes (response terms) respond to changes in the LW↓ and net SW flux (forcing terms).
The response terms are not always governed by the forcing terms, as, for instance, under high wind conditions the turbulent heat fluxes can operate independently as the Ri in these cases is dominated by the wind shear (see Equation 7). Cloud presence 5 influences the radiational balance at the surface by modulating the downwelling radiation; increasing LW↓ and decreasing SW↓. Miller et al. (2015) show that clouds increase the net surface radiation compared to an equivalent clear-sky scene, because the high year-round surface albedo limits the magnitude of the cloud SW cooling effect to less than that of the LW warming effect. Statistical relationships for the current study reinforce the fact that liquid-bearing clouds increase the forcing terms during two distinct periods; with and without solar insolation ( Figure 7a). Hence, the occurrence of liquid-bearing clouds 10 correspond to warmer surface temperatures in both circumstances ( Figure 7b) and consequently greater LW↑ (Figure 7c), which is proportional to the surface temperature to the fourth power.
LW↑ has less variability (all cases in Figure 7c) than the variability of the forcing terms (all cases in Figure 7a). In addition, the differences between the cloudy and non-cloudy states are more pronounced in Figure 7a, compared to Figure 7c

Boundary-Layer Stability Response
The degree to which the overlying atmosphere can dynamically interact with the surface is important for determining the 20 turbulent heat exchange. Atmosphere/ice sheet interaction is modulated by low-level stability, which can be influenced by both thermodynamic and dynamic processes. Mechanical mixing, via high wind speeds, is one way to decrease near-surface stability and increase turbulence near the surface. The 10 m wind speed is greater than 8 ms −1 for 16% of 32130 stability estimates.
The median Richardson number decreases from 0.19 for all cases to 0.06 for the cases that report higher wind speeds (>8 ms −1 ), showing the expected decrease of stability. In addition, cloud-driven atmospheric mixing can also affect the low-level 25 atmospheric structure (Shupe et al., 2013a) and liquid-bearing cloud presence, especially in combination with enhanced solar radiation, decrease the near-surface temperature gradient (Hudson and Brandt, 2005;Miller et al., 2013).
This study explicitly shows that the radiative influences of liquid-bearing clouds and/or insolation create neutral or unstable boundary-layer conditions. When the sun is below the horizon, as for the first case study (Section 3.2.1), the presence of liquidbearing clouds decreases the stability such that a majority of the cases are weakly stable (0 < Ri < 0.25) (Figure 7d). In the stable. However, when optically thick liquid-bearing clouds (LWP > 30 gm −2 ) are present the boundary-layer is near-neutral on average. Interestingly, optically thin liquid-bearing clouds (5 gm −2 > LWP > 30 gm −2 ) lead to more frequent occurrence of more unstable conditions in the presence of insolation, because these clouds emit significant longwave radiation while also allowing significant penetration of solar radiation, thus producing the maximum surface heating. Our results that liquid-bearing clouds of intermediate thickness lead to higher instability agree with studies that show these clouds produce the maximum cloud 5 radiative forcing for elevated sun angles (Bennartz et al., 2013;Miller et al., 2015). Hence, liquid-bearing clouds and/or solar insolation enhance turbulent mixing, facilitating sensible and latent heat exchange, although instability (negative Ri) requires SW↓.

SEB Responses
Process-based relationships distill our understanding of the underlying physical processes into a succinct form that is infor-10 mative, yet practical. While clouds, the solar cycle, and other processes can influence the downwelling radiation, process relationships between response terms and forcing terms reveal how variability in downwelling radiation affects the other SEB terms. A linear regression of the relationship between the forcing and response terms yields a slope of -1.01 (Figure 8a), indicating that the SEB is largely radiatively driven, the response terms account for all of the forcing energy flux, and there is approximate closure for the SEB terms calculated here. The scatter in this relationship is due to measurement uncertainties, 15 mismatches of response times in different terms, and the spatial distribution of the instrumentation. The annual evolution of this slope (Figure 9) shows that the SEB response terms balance the forcing terms to within ≈10% in all months of the year.
Thus, any change in forcing terms elicits an approximately equal change in flux through the combination of response terms.
The response to the radiative forcing can be evaluated for each term independently (Figure 8b-e), and as a function of month, showing that each term responds differently throughout the annual cycle ( Figure 9). The slope of the linear fit provides 20 an estimate of the relative magnitude (percentage) of the response of each term. The RMS error of the monthly response estimates in Figure 9 are calculated by comparing the estimated values, using the linear fit, to the measured values ( Figure 10).
The RMS error includes the uncertainty of the measurements involved, any delay in response time greater than 30 minutes, and variability in the physical response not represented by the linear fit. Generally, the RMS error of the linear fits of all response terms to the driving terms are on the same order of magnitude as the combined uncertainty of the SEB components. 25 The annual response in the LW↑ term (Figure 8b to the increased response of the latent heat flux for this specific month (Figure 9). Any increase (decrease) of response of 30 an individual term will effectively decrease (increase) the change in surface temperature, and hence the response of LW↑, to radiative forcing.
The response of the sensible heat flux is fairly constant at ≈ -0.11 throughout the annual cycle ( Figure 9) due to its dependence on both the near-surface temperature gradient and stability (heat transfer coefficient -see Equation 5). For weakly stable conditions, the former term dominates decreasing (increasing) the heat flux for surface warming (cooling), while for very stable conditions the latter term dominates limiting turbulent exchange and increasing (decreasing) the sensible heat flux for surface warming (cooling) (e.g., Grachev et al., 2005). Since this Summit data generally shows a decrease in sensible heat flux for an increase in the forcing terms (surface warming), this is consistent with weakly stable conditions on the unstable side of the stability transition shown by Grachev et al. (2005). Therefore, the response of the sensible heat flux to changes in the 5 surface temperature is similar throughout the year and does not show an annual cycle. However, the RMS error of the linear fit ( Figure 10) during winter (9.7 Wm −2 ) is greater than during summer (6.0 Wm −2 ) (i.e., there is more scatter in the sensible heat response in winter), suggesting that conditions in winter are at times very stable and that the sensible heat flux response to radiative forcing is then different. In summer, conditions are rarely very stable so the response in sensible heat flux is more strongly correlated with the change in the forcing terms.

10
The response of the latent heat flux increases in summer compared to other months of the year (Figure 9). The amount of available moisture (Figure 6c) peaks in summer and average PWV values for non-summer (winter) months are below 2 mm (1 mm). Thus, changes to near-surface stability due to changes in the forcing terms produce a smaller response when moisture gradients are small in magnitude.
The response of the conductive heat flux to radiative forcing is greatest in winter (Dec -Feb) at 22% compared to ≈10% 15 in summer (June-August). Seasonal changes in the conductive heat response are due to changes in snow density, thermal conductivity, and subsurface temperatures. Warmer subsurface temperatures resulting from prior warm surface temperatures precondition the snowpack, reducing its ability to remove heat from the surface. Decreased density in the summer decreases the thermal conductivity of the near-surface snow pack, also limiting the ability of the subsurface to remove energy from the surface. The RMS error of the linear fit of the conductive heat flux to the forcing terms is relatively low with an annual mean 20 of 3.2 Wm −2 .
The response of the heat storage in the upper subsurface layer is important to consider when accounting for all the energy responses at the surface. Even though the annual mean of S is less than 1 Wm −2 (i.e., there is effectively no annual net change in temperature in the near-surface snow), it is highly variable (annual standard deviation = 62.5 Wm −2 ) as this layer can warm or cool rapidly from one half hour period to the next. The heat storage response to the forcing terms also accounts for 25 subsurface heating due to solar penetration. Over the annual cycle the response of S is 7%. The response of S is at a maximum in April -May at 23-25% (Figure 9), indicating relatively cold near-surface snow is able to store larger amounts of energy originating from radiative sources.
Since scatter in S in response to forcing is so large, we first examine the scatter of all the other terms jointly. The RMS error of the linear fit of (LH + SH + C -LW↑) vs. (LW↓ + SWnet) is maximum in July (19.6 Wm −2 ) and has an annual mean value S -LW↑ to the forcing terms, shown by the error bars in Figure 9, is primarily due to the variability and associated uncertainty in S. However, correctly accounting for the ground heat flux in the upper most layer provides near closure of the surface energy balance, a critical accomplishment of the synthesis of comprehensive datasets given here.
At the ice sheet/atmosphere interface surface temperature is the linchpin that connects the subsurface to the atmospheric boundary layer, responding to changes in the net flux at the surface. The variability in the surface temperature is controlled 5 by changes in the forcing terms and modulated by the response terms. An increase in radiative forcing warms the snowpack; increasing the surface temperature and decreasing the near-surface atmospheric stability. Not surprisingly, the response terms are all associated with surface temperature; either directly proportional, or a function of the near-surface temperature gradient.
Latent heat flux is also dependent on the near-surface moisture gradient and the ground heat flux is dependent on the thermal conductivity of the snow pack leading to seasonal differences in their responses. This study highlights the importance of the 10 seasonal changes in the non-radiative responses, which determine the annual cycle of the LW↑ response.
The seasonal response of the SEB to cloud presence is estimated by combining the radiative effects of clouds with the observationally based and statistically derived relationships between the forcing and response terms. Cloud radiative forcing (CRF) at the surface, as detailed in Miller et al. (2015), is the instantaneous net radiative effect of clouds. Furthermore, changes in the forcing terms elicit a response of the surface temperature and the non-radiative SEB terms. Thus, we combine the monthly 15 CRF values reported in Miller et al. (2015) and monthly responses, calculated from the maximum available data (Figure 9), to estimate the corresponding increase in LW↑ and decreases in SH, LH and G attributed to cloud presence. Figure 11a shows LW↑ has the smallest increase due to CRF in May (9.8 Wm −2 ), the largest increase in October ( Over the annual cycle atmospheric temperatures in the free troposphere (> 1 km) are well correlated with surface tempera- 30 tures, although energy exchange processes at the surface enhance surface temperature variability. In general, time-series data, monthly mean values, and monthly diurnal cycles all show that the non-radiative SEB terms oppose the increase or decrease of the net radiation. Liquid-bearing clouds and solar insolation strongly modulate the radiative flux that reaches the surface, which affects subsurface temperatures, the stability of the boundary-layer, and the near-surface temperature gradients. A pair of case studies illustrate how all the pieces fit together to depict how an increase in surface radiation elicits a response in the surface temperature, while also indicating that the increase in temperature is lessened by a decrease in sensible and conductive heat fluxes. The resultant compensation of the non-radiative SEB terms thereafter affects the net amount of surface warming that occurs due to cloud radiative forcing and/or insolation. Similar compensation is apparent when looking at longer-term 5 averages.
To examine these relationships in more detail, radiative forcing terms (LW↓ + net SW) were related to the response terms (SH, LH, C, S and -LW↑) throughout the annual cycle. Linear regression analysis, for the year-round dataset relating the response terms as a function of the forcing terms, resulted in a -1.01 slope, indicating general closure in the calculated SEB terms. On average LW↑, which is directly linked to surface temperature, responds by about 70% of a perturbation in incident 10 radiation, with a somewhat diminished response in summer. Quantifying how each non-radiative response changes throughout the year provides insight into how much SH, LH and/or G limit the surface temperature increase due to the occurrence of liquid-bearing clouds and/or insolation: -Latent heat flux response is near-zero for much of the year, with an increased response in summer.
-Sensible heat flux response is fairly constant throughout the annual cycle (11%).

15
-Ground heat flux, consisting of both heat storage in the upper most ≈ 20 cm of snow and the conductive flux below this layer, is the largest non-radiative response for most of the year, with a decreased response in summer.
The enhanced summer latent heat flux response is due to an increase in available moisture and an increase in turbulence during relatively frequent periods of neutral/unstable near-surface conditions. In winter the effect of the stable boundary-layer is to dampen the response of the turbulent sensible heat flux, yet this dampening effect is offset by the enhanced near-surface 20 temperature gradient. The consequence of a limited sensible heat exchange during periods of strong radiational cooling is that the sensible heat flux response is relatively constant throughout the annual cycle. Finally, the ground heat flux response decreases in the summer due to decreases in near-surface snow density and warmer subsurface temperatures.
Previous studies by Cullen (2003)  annual cycle of the SEB variations, due to the atmosphere being relatively warm and cloudy, the monthly mean values of this study in the winter are similar to the earlier studies. However, the summer monthly total radiation is lower in our annual cycle (July 2013 -June 2014) compared to previous studies, thus the recent monthly averages for non-radiative terms are greater than the 2000 values. The differences in the annual cycles could be due to possible decreases in cloud cover (Comiso and Hall, 2014). 30 In central Greenland, cloud presence in winter (longwave forcing) is unable to produce a neutral stratification. It is only with insolation that neutral and unstable conditions exist. In contrast, over Arctic sea ice, wintertime conditions are near-neutral or even slightly unstable nearly 25% of the time (Persson et al., 2002). More instability over sea ice compared to the Greenland ice sheet may be due to warming of the surface from below due to oceanic influences. Springtime/summertime near-neutral and slightly unstable conditions with shortwave forcing observed here is similar to that observed over sea ice (e.g., Ruffieux et al., 1995;Persson et al., 1997Persson et al., , 2002. Also in agreement with our findings are process diagrams obtained via a modeling study over sea ice (Sterk et al., 2013) that found the non-radiative SEB terms lessen the change in surface temperature due to changes in downwelling radiation. Moreover, observational studies over sea-ice (Persson et al., 2002) and in the Greenland ablation 5 zone (van den Broeke et al., 2011) suggest that if/when Summit Station more frequently experiences melt the non-radiative compensation, detailed in this study, may be significantly diminished as energy goes towards surface melting.
These central Greenland results can be used to evaluate how well the annual and diurnal cycles of the SEB terms are represented in climate models and reanalyses, and specifically the relationship among key terms. It is known that global climate models underestimate the occurrence of liquid-bearing clouds above Greenland (Kay et al., 2016). We estimate that 10 the underrepresentation of clouds, especially liquid-bearing clouds, should produce annual surface temperature biases ranging from 0 to -6.9 • C. If the representation of liquid-bearing clouds were to improve then the modeled downwelling radiation would likely also improve, but it is unclear if the other SEB terms would realistically adjust. A regional or global climate model's modus operandi is to achieve absolute closure of the SEB, hence this study will be useful in future studies as a valuable tool for pinpointing the processes responsible for possible model surface temperature bias over Greenland and for evaluating model 15 representation of physical processes at the ice sheet/atmosphere interface.  Distributions are represented by box-and-whisker plots (the box indicates the 25th and 75th percentiles, the whiskers indicate 5th and 95th percentiles, the middle line is the median, and the * is the mean).