EVALUATION OF DIFFERENT METHODS TO MODEL NEAR-SURFACE TURBULENT FLUXES FOR AN ALPINE GLACIER IN THE CARIBOO MOUNTAINS , BC , CANADA

As part of surface energy balance models used to simulate glacier melting, choosing parameterizations to adequately estimate turbulent heat fluxes is extremely challenging. This study aims to evaluate a set of bulk methods commonly used to estimate turbulent heat fluxes for a sloped glacier surface. The methods differ in their parameterizations of the bulk exchange coefficient that relates the fluxes to the mean meteorological variables measured 2 m above a glacier surface. The performance of 23 bulk approaches in simulating 30-min sensible and latent heat fluxes is evaluated against the measured fluxes from an 5 open path eddy-covariance (OPEC) method. The evaluation is performed at a point scale of an alpine glacier, using one-level meteorological and OPEC observations from a multi-day periods in the 2010 and 2012 summer season. The analysis of the two independent seasons yielded similar findings, listed as following. The bulk method, with or without the commonly used Monin-Obukhov (M-O) stability functions, overestimates the turbulent heat fluxes over the observational period, mainly due to an overestimation of the momentum flux. In the absence of OPEC-derived M-O stability parameter, no method can successfully 10 predict this parameter, which results in poor performances of the M-O stability corrections and consequently the bulk method. The OPEC-derived 30-min momentum flux is linearly related to the measured wind speed, contrary to the proposed quadratic relation by the commonly used bulk methods. An approach based on a katabatic flow model, which assumes a linear relation between the shear stress and the wind speed, outperforms any other bulk approach that we tested in simulating the momentum flux. In agreement with the katabatic flow model, we show that in a more stable atmosphere the bulk exchange coefficient 15 for momentum is smaller. The sensible heat flux can be more successfully modeled if the bulk exchange coefficients for momentum and heat are allowed to follow different parametrization schemes, rather than assuming equal schemes as is the case in the common bulk methods. Further data from different glaciers are needed to investigate any universality of these findings. 1 The Cryosphere Discuss., doi:10.5194/tc-2017-80, 2017 Manuscript under review for journal The Cryosphere Discussion started: 31 May 2017 c © Author(s) 2017. CC-BY 3.0 License.


Introduction
Mountain glaciers of British Columbia are experiencing significant mass losses in response to ongoing climate change and are projected to lose most of their current volume by the end of the century (Radić et al., 2014;Clarke et al., 2015).These projections, however, are highly sensitive to the representation of glacier mass balance in the models, in particular to the model parameterizations used to simulate surface melting.All regional to global scale projections of glacier mass changes rely on a relatively simple representation of glacier melting via so-called 'degree-day' or temperature-index' models (Radić and Hock, 2014).Despite relying on empiricism, these models have relatively successfully dealt with the absence of fine-resolution meteorological input data that are, on the other hand, a prerequisite for a good performance of more physics-based surface energy balance (SEB) models.However, with the ongoing rapid development of computational recourses and global climate models, it is only a matter of time for SEB models to become routinely used in the studies of glacier evolution on regional and global scales.To ensure that future models are firmly routed in theory, new methods are required to test and parameterize turbulent heat fluxes at the glacier surface.
The turbulent sensible and latent heat fluxes are recognized as important components of the SEB over mid-latitude glaciers worldwide.When averaged over a longer period (e.g.weeks, months) these fluxes are generally less than net radiation flux.
Over daily to hourly time scales, however, they can exceed the net radiation (Hock, 2005).Highest melt rates at these short temporal scales often coincide with times of greatest turbulent heat flux (Hay and Fitzharris, 1988).Direct measurements of turbulent fluxes, by eddy-covariance method, are relatively rare because they require sophisticated instrumentation with continuous maintenance which makes them unsuitable for long-term operational purposes.Few studies exist that have directly measured turbulent fluxes over mountain glaciers and all of these studies are restricted to shorter time intervals, e.g.weeks to couple of months (e.g., Munro, 1989;Forrer and Rotach, 1997;van der Avoird and Duynkerke, 1999;Cullen et al., 2007;Conway and Cullen, 2013).A more common approach in the studies of SEB on glaciers is to derive the turbulent fluxes through parameterization schemes that utilize the observations of mean meteorological variables such as near-surface air temperature, wind, and relative humidity (e.g.Hock and Holmgren, 2005;Sicart et al., 2005;Mölg et al., 2008;Reijmer and Hock, 2008;Mölg et al., 2009;MacDougall and Flowers, 2011;Sicart et al., 2011;Huintjes et al., 2015;Prinz et al., 2016).The simplest and most widely applied method is the bulk aerodynamic method (bulk method) (e.g., Braithwaite et al., 1998;Oerlemans, 2000).
The bulk method is an integrated form of a gradient-flux relation, based on the theoretical work of Prandtl (1934) and Lettau (1934), which in the surface boundary layer assumes constant vertical fluxes of momentum, heat and humidity and horizontal homogeneous conditions of the surface.
Despite the common usage of bulk methods in glacier SEB studies, there are many uncertainties in its application.We summarize the main limitations of bulk methods as applied to glacier SEB studies: 1. Surface conditions.The bulk method for sensible heat fluxes relies on the assumption that the difference between near surface air temperature, often measured at 2 m above the surface, and surface temperature well represents the near-4.Stable boundary layers accompanied by katabatic flow.Over a sloping glacier surface, a low-level katabatic jet is a common feature, which means that buoyancy enters the horizontal momentum equation and the turbulence arises independently of the surface roughness (e.g., van der Avoird and Duynkerke, 1999).The Prandtl model for katabatic flow that treats eddy viscosity as a constant value (Prandtl, 1942) is unable to correctly describe the sharp near-surface gradients in wind speed and air temperature that are often observed (Munro, 1989;Oerlemans, 1998).Grisogono and Oerlemans (2001) showed that the Prandtl model can be improved if a varying assigned eddy viscosity profile is used instead of a constant value.Using this 'enhanced ' Prandtl model, Parmhed et al. (2004) found that modeled turbulent fluxes at a glacier surface compare well with the measured fluxes during katabatic flows.An alternative approach to incorporate the model for katabatic flow into the bulk method for turbulent heat fluxes was developed following the classical Prandtl model for slope flows (Oerlemans and Grisogono, 2002).The model is simplistic and therefore agrees only in the overall pattern that shows an increase in sensible heat flux in response to an increase in air-surface temperature difference.Despite its simplicity, application of the model requires calibration with detailed observations of temperature and wind profiles in the surface boundary layer, in addition to the near-surface meteorological observations.Few such measurements exist on glaciers and therefore these models have not been rigorously evaluated so far.
In the light of the above, further research is needed to evaluate and develop methods that are suitable for determining nearsurface turbulent fluxes over a sloping glacier surfaces subject to katabatic flows.In particular, our objective is to evaluate the existing framework of bulk approaches commonly used for parameterizing the turbulent momentum and heat fluxes at a glacier surface.We achieve this by comparing the modeled fluxes, i.e. output of the bulk approaches, to measured fluxes, i.e. the fluxes derived from an open-path eddy-covariance (OPEC) method.OPEC measurements are obtained from an alpine glacier in British Columbia (BC), over a short-term window (weeks) for two summer seasons.In the sections to follow we start with a brief overview of the field site, data collection, eddy-covariance data treatment and overall methodology.This will be followed with a detailed description of the bulk approaches and their evaluation results, and finalized with the discussion and conclusions.
2 Data and Methods

Study site and station setup
The study uses data collected from Castle Creek Glacier, an alpine glacier in the Cariboo Mountains, BC (Fig. 1).This 9.5 km 2 mountain glacier flows north for 5.9 km, has an elevation range of 2827-1810 m a.s.l. and contributes meltwater to Castle Creek, a tributary of the Fraser River.The glacier's annual mass balance is monitored since 2009 (Beedle et al., 2014).We utilized the data from an automatic weather station on a glacier (AWS glac ): in summer 2010, the AWS glac was installed at 53 • 3' 2.99" N and 120 • 26' 39.57" W and an altitude of 1967 m a.s.l, while in summer 2012 the AWS glac was installed within ± 10 m of the chosen position in 2010.In the glacier vicinity, at the lateral and terminal moraines, two year-round automatic weather stations have been in operation since 2007/2008 (Déry et al., 2010).The stations are referred to as AWS up (53 • 2' 36" N, 120 • 26' 18" W, 2105 m a.s.l.) and AWS low (53 • 3' 45" N, 120 • 26" 4" W, 1803 m a.s.l.).At AWS up , the mean annual air temperature over 2007-2010 was -2.6 • C, while summer (Jun-Aug) mean air temperature was 6.6 • C.Over the same period, precipitation during summer (Jun-Aug) averages 94 mm at AWS up although rainfall is likely underestimated due to gauge undercatch at the exposed ridge site where mean monthly wind speeds often exceed 5 m s −1 (Déry et al., 2010).The lower part of the glacier (< 2100 m a.s.l.), where the AWS glac is located, is gently sloping with an approximate mean gradient of 7 • .
Meteorological data were sampled at 10-second increments, and 10-minute averages were recorded with a Campbell Scientific CR1000 datalogger.Temperature and relative humidity were measured at a height z t =1.7 m above the ice surface with a Rotronic T/RH sensor, while wind speed and direction were measured at z v =2.0 m (2010 season) and z v =1.9 m (2012 season) 10 with a RM Young wind monitor.In addition to these variables, radiation fluxes (incoming and reflected shortwave-, incoming longwave-, and net radiation) were also recorded.Surface air pressure and liquid precipitation were recorded at the AWS up and AWS down and the measured surface air pressure values were linearly interpolated to represent the surface air pressure at AWS glac .The meteorological sensors at AWS glac , mounted on a floating tripod (Fig. 2), maintained a constant measurement height above the surface over the observational period.At a 3 m distance from the floating tripod, a Campbell Scientific openpath eddy covariance (OPEC) system and sonic ranger (SR50) were installed on a separate tower drilled into the ice (Fig. 2) (Jarosch et al., 2011).As part of the OPEC system, a sonic anemometer (CSAT-3) and a krypton hygrometer (KH2O) were mounted with an initial height of 1.86 m (2010 season) and 2.0 m (2012 season) above the surface.

Eddy-covariance data processing
We applied a series of data processing steps as recommended by Burba (2013) that are part of the Eddy-Pro Software used to process the OPEC data and calculate the turbulent fluxes.The Eddy-Pro processing output yields turbulence statistics and calculated fluxes for each 30-min segment of the raw data.The height of OPEC sensors relative to the ice surface changed due to the surface melting throughout the observational period: the net change of 0.70 m over 12 days in 2010, and 0.75 m over 29 days in 2012.We applied the following corrections to the raw EC data and, as later tested, none of them accounting for differences in the sensor's height produced substantial (>1 %) differences in 30-min averaged turbulent fluxes.Threedimensional (3-D) wind speed, sonic temperature and vapor density fluctuations were recorded at 20 Hz on the datalogger (vapor density fluctuations were later converted to specific humidity fluctuations).At the beginning of each observational season we leveled the instruments to horizontal, and have recored a minimal tilt at the end of each observational period.In both seasons, we set CSAT-3 so that the sensor's arm formed a 90 • angle to the prevailing down-glacier wind direction to minimize flow distortion due to air flow through the station structure.Before calculating the turbulence statistics, we applied a coordinate rotation to the velocity data to align the streamlines into the mean flow using a planar fit method (Wilczak et al., 2001).This method sets the vertical axis of the CSAT-3 perpendicular to a hypothetical plane produced from averaged wind measurements during the observation period.Because the glacier surface undergoes changes with time, we applied sensitivity tests to assess The Cryosphere Discuss., doi:10.5194/tc-2017Discuss., doi:10.5194/tc- -80, 2017 Manuscript under review for journal The Cryosphere Discussion started: 31 May 2017 c Author(s) 2017.CC-BY 3.0 License.whether the alignment of the plane needed adjustment and whether a different coordinate rotation method (e.g.double or triple rotation scheme) produced different output.We recored no substantial changes (>1 % difference) in the 30-min turbulent fluxes in these tests.The latent heat flux data were corrected for oxygen absorption by the KH2O (Tanner et al., 1993) and for the fluctuations in the water vapor density measurements not caused by turbulent eddies (Webb-Pearman-Leuning correction; Webb et al., 1980).Finally, we applied a series of standardized processing steps consisting of corrections for potentially high and low frequency loss, and sensor separation (Ibrom et al., 2007;Moncrieff et al., 2004;Horst and Lenschow, 2009).

The algorithm and quality control for roughness lengths
The covariances of the 3-D wind velocities, temperature and specific humidity, derived by the OPEC system for each 30-min data segment, were transfered to the friction velocity (u * ), temperature scale (θ * ) and specific humidity scale (q * ) using: where u and v are fluctuations of the horizontal wind component around their 30-min mean values, w is fluctuation of the vertical wind component, and T and q are the temperature and specific humidity fluctuations.The OPEC system was also used to indirectly measure the Obukhov length (L): where T v is the 30-min averaged virtual temperature (K), g is the acceleration due to gravity (9.81 m s −2 ), and κ (0.40) is the von Kármán constant.Following Cullen et al. (2007), the turbulent scales and L were then used to derive roughness lengths for each 30-min data segment: where U zv , T zt and q zt are 30-min averages for the wind speed (m s −1 ), air temperature ( • C), and specific humidity (kg kg −1 ) respectively, determined from AWS glac measurements at the senors' heights (z v,t ).In the absence of direct measurements, the surface temperature (T 0 ) was assumed to be at melting point (0 • C) and the surface vapor pressure at saturation (6.13 hPa).
The assumption of consistent melting is corroborated with the sonic ranger measurements showing persistent surface lowering throughout the observational period.In general, assuming that T 0 = 0 • C works well on temperate glaciers during a melting season, and is more accurate than estimating the surface temperature from the longwave radiation measurements (Fairall et al., The Cryosphere Discuss., doi:10.5194/tc-2017-80, 2017 Manuscript under review for journal The Cryosphere Discussion started: 31 May 2017 c Author(s) 2017.CC-BY 3.0 License.
1998) or from a SEB closure (Hock, 2005).We nevertheless quantify errors in our results due to assumed rather than measured surface conditions.Vapor pressure (e) was converted to specific humidity using q = 0.622 e p , where p is observed air pressure (hPa).Ψ v ( zv L ), Ψ t ( zt L ) and Ψ q ( zt L ) (where Ψ t ( zt L ) = Ψ q ( zt L )) are integrated forms of universal functions based on M-O theory.Many such functions have been developed, and the most widely used is the Dyer-Businger flux-profile relationship (Businger et al., 1971;Dyer, 1974).Following Conway and Cullen (2013), for stable stratification ( zv,t L > 0) we use the expressions of Holtslag and de Bruin (1988) in the following forms: where a = 1, b = 2 3 , c = 5, and d = 0.35.For unstable stratification ( zv,t L < 0), which was rarely present during our observational periods, we use the expressions of Dyer (1974).
Turbulence flux data from OPEC systems often contains spurious measurements.The roughness lengths derived from Eq.
(5) -( 7) widely vary because they rely on several mean and turbulence variables.One way to reduce this scatter is to limit analysis on high-quality data (Andreas et al., 2010).In this study we applied a series of filters to the 30-min data segments (points) obtained during the study period following the approach of Conway and Cullen (2013) and Li et al. (2016).Our filtering algorithm proceeds as follows: (a) 'Basic' filter for z 0t and z 0v : removed the points for which < 0, which produce unrealistically large roughness lengths since the exponents in Eq. ( 5) -( 7) become much larger than unity.One of the main reasons for obtaining the negative values is that the bulk approach (finite differences) is a first-order approximation of the gradientflux approach, and the surface values for T 0 and q 0 are assumed and kept constant in time rather than measured directly for each 30-min segment.
(b) 'Stationarity' filter: only steady-state runs were used following the method of Foken (2008)  (e) 'Wind speed' and 'u * ' filters: runs were selected for U zv > 3 m s −1 and u * > 0.1 m s −1 as errors in deriving the roughness lengths become comparatively large for low wind speeds and small friction velocities.In addition, accounting for a relatively strong air flow (U zv > 3 m s −1 ) helps reduce potential errors in T zt due to radiative heating of temperature sensor during periods of strong solar radiation (Huwald et al., 2009).(f) 'Temperature' gradient' filter for z 0t : to ensure that a sufficiently large temperature gradient is detected between the surface and measurement height, we selected the runs for which T zt > 1 • C.
(g) 'Moisture gradient' filter for z 0q : similarly as above, to ensure for sufficiently large moisture gradients between the surface and measurement height, we selected the points that satisfy |e zt − e 0 | > 0.66 hPa, corresponding to a difference between T zt =0 • C and T zt =1 • C.
(h) 'Small values' filter applied to z 0t and z 0q : following Andreas et al. (2010) we assumed that the surface exchange of heat and moisture cannot occur at scales smaller than 10 −7 m.This scale is an approximate mean free path of air molecules at sea level.Considering that this is the minimal scale at which molecular diffusion takes place at the surface, z 0t and z 0q smaller than this are treated as unrealistic (i.e.eddy diffusivity can take place only at larger scales).
(i) 'Large values' filter: no study of turbulent fluxes over glaciers determined roughness lengths larger than 1 m, thus we treated any values of roughness length that exceed 1 m as unrealistic for our site.
We note that it is common to exclude the data during precipitation events as these greatly increase measurement errors associated with the OPEC system.Nevertheless, this filter was redundant in our study because the above filters already removed the points overlapping with the precipitation events (recorded at AWS low ), and because OPEC system (in particular KH20) failed to take measurements during most of these events.

Bulk methods
All bulk methods used in this study are rooted in the gradient-flux relations, in which the turbulent fluxes of momentum (τ ), sensible heat (Q H ) and latent heat (Q E ) in the surface boundary layer are proportional to the time averaged gradients of wind speed (U ), potential temperature (θ) and specific humidity (q), expressed by: where ρ a is the air density, c p is the specific heat of air at constant pressure (1005 J kg −1 K −1 ), L v is the latent heat of vaporization (2.514 MJ kg −1 ), z is the height above the surface, K M is the eddy viscosity, and K H and K E are the eddy diffusivities for heat and vapor exchange, respectively.Eq. ( 10) also defines the friction velocity u * which we use henceforth as a surrogate for surface shear stress.We note that by convention in SEB on glaciers, the heat fluxes are positive if they transport heat towards the surface, and negative if they transport heat away from the surface.

Methodology outline
Before we go into details, we provide an overview of the bulk approaches, whose performance in simulating the turbulent fluxes of momentum, heat and humidity, we evaluate in the study.In total we evaluate 23 bulk schemes, and the details of each scheme are provided together with the results, in order to facilitate easier readability.The modeled 30-min turbulent L .These schemes require estimates of the roughness lengths (z 0v,t,q ) which we assume to be constant in time and equal to the mean log value of the OPEC-derived 30-min roughness lengths for neutral condition only (i.e. when all filters are applied to the OPEC data).In addition to these three schemes, we use C M −O parametrization with the scalar roughness lengths derived from the surface renewal model of Andreas (1987), yielding (4) C SR parameterization scheme.Second, we derive the turbulent fluxes using the K-approach method, but instead of modeled friction velocity (u * ) and modeled zv,t L we use their OPEC-derived values.This yields the following schemes: (5)  Oerlemans and Grisogono (2002), and an approach with an integrated varying eddy viscosity profile K Int (Grisogono and Oerlemans, 2001;Parmhed et al., 2004).In the K Int approach where we analyze whether two parameters that determine K Int (later introduced as K max and H K ) can be expressed as functions of OPEC-derived zv,t L .This yields two schemes: Finally we test a 'hybrid' approach which combines the K-approach with the K Int method (with OPEC-derived or modeled

Estimates of the roughness lengths
We apply the filters discussed in Section 2.3 to the OPEC data to obtain the high-quality estimates of 30-min turbulent fluxes, which are then used to determine the roughness lengths using the Eq.( 5)-( 7).For the 12-day observational period in 2010, prior to the filtering, there are 515 values for z 0v,0t while 307 values for z 0q .The fewer initial data points for z 0q is due to high sensitivity of KH2O sensor to the build-up of droplets on the sensor window during precipitation events.During these events KH2O mainly produced spurious values or failed to record the data, while CSAT-3 recorded values, albeit with questionable quality.After filtering, usable data reduced by 98 %, giving nine high-quality estimates for z 0v,0t , but none for z 0q (Fig. 3).
Similarly, initial number of data points in the 29-day observational period in 2012 is 1161 for z 0v,0t and 621 for z 0q , while the post-filtering number of data points are 30 for z 0v , 25 for z 0t and six for z 0q .Among all the filters, the 'neutrality filter' (| zv,t L | < 0.1) eliminates the largest percent of the data points, followed-by the 'stationarity' filter.Despite the significant loss of viable data, these filters are essential to achieve a more accurate determination of roughness lengths with the bulk method.
OPEC-derived fluxes are assumed valid only during the steady-state or stationary conditions (i.e. in the absence of intermittent turbulence and/or gravity waves).For the neutral conditions, the choice of the stability correction functions has negligible effect on the calculated fluxes (the correction function approaches zero as | zv,t L | approaches zero).Furthermore, the neutral stability conditions are more commonly observed during overcast and windy conditions rather than during clear-sky conditions with prevailing katabatic flows.This fact minimizes the chance of selecting 30-min intervals with katabatic flows for which the bulk method, used to determine the roughness lengths, might be deficient.If we remove the 'neutrality' filter, and instead select all values for which | zv,t L | < 2, we obtain 196, 191, and 68 values for z 0v , z 0t , and z 0q , respectively in 2010, while for 2012 season we obtain 276, 259, and 60 values for z 0v , z 0t , and z 0q , respectively.The threshold of zv,t L = 2 is chosen because the universal stability functions for stable stratification are commonly defined up to z L = 2 (Foken, 2008), which represents strongly stratified stable regime.
We determine the time windows of prevailing katabatic flows (shaded areas in  3 m s −1 ), preferably reaching maximum 30-min mean values > 5 m s −1 .With this rough identification criteria we may miss some shorter time intervals with katabatic flows, but a more restrictive approach reduces the uncertainty that originates from the selection of katabatic flow periods in our further analysis.The identified prevailing katabatic periods (Fig. 4) overlap with clear sky conditions and with relatively high near-surface air temperatures (T zt ), as is commonly expected for katabatic flows at glacier surfaces (e.g., Oerlemans et al., 1999).For each roughness length estimate from 30-min segments (Eq.( 5) -( 7)) we also estimate its error (error bars in Fig. 4) derived through an error propagation method for a multi-variable function (e.g., Bevington, 1969), where the measurement error for wind speed and assumed errors for surface temperature and water vapor pressure are propagated (keeping just the linear terms in Taylor expansion) into the final error for the roughness length.The following error values are used: measurement error for wind speed δU zv =0.11 ms −1 , assumed error for surface temperature δT 0 =0.5 • C, and assumed error for surface water vapor pressure δe 0 = 0.23 hPa.To provide an expected (mean) value for  roughness lengths during neutral conditions we took the average of their logged values given the measurement's log-normal distribution.To incorporate the individual error of each estimate into the mean value, we derive a weighted mean roughness length, where the weights are inversely proportional to the individual errors.The results (Table 1) show that the weighted mean values do not substantially differ from the mean values and we thus use the non-weighted mean values in the remaining analysis.For both seasons, our results show that mean log z 0v is of the order of 10 −3 m and two orders of magnitude larger than mean log z 0t .The overall error of these mean values, expressed as one sigma (± a standard deviation), shows uncertainty of one order of magnitude around the mean value (Table 1).To determine z 0q for 2010 season, which had no high-quality measurements for z 0q ; we thus assume z 0q ≈ z 0t since the equality holds for the mean log values in the 2012 season.

K-approach with common parameterizations for C
Following Conway and Cullen (2013) we use the bulk methods derived by integrating the flux-gradient equations with three different assumptions: (1) C log method -assumes that wind speed, temperature, and humidity have logarithmic profiles with height (2) C Rib method -assumed that the logarithmic profiles are modified with stability corrections based on the bulk Richardson number (R ib ), (3) C M −O method -assumes that the logarithmic profiles are modified with M-O universal stability functions of zv,t L , and (4) C SR -same as C M −O but with modeled scalar roughness lengths.In the following sections we first provide a brief theoretical background and derivation of each method, which is then followed by an evaluation of those approaches.

K-approach with logarithmic profiles for wind and temperature (C log )
According to the mixing-length theory (Prandtl, 1934) for a neutral surface layer, eddy viscosity can be parametrized as (Kapproach; Stull, 1988): The Cryosphere Discuss., doi:10.5194/tc-2017-80,2017 Manuscript under review for journal The Cryosphere Discussion started: 31 May 2017 c Author(s) 2017.CC-BY 3.0 License.
Assuming K M = K H = K E , inserting these parameterizations in the gradient flux relations (Eq.( 10) -( 12)) and integrating the equations with respect to height, while treating u * , Q H , Q E as constants, gives the bulk aerodynamic expressions: where dimensionless exchange coefficients for the sensible and latent heat flux, respectively.Wind speed (U zv ), air temperature (T zt ) and water vapor pressure (e zt ) represent the time averaged values at their measurement heights z v,t , while zero subscript in these variables denotes their surface values (note that z t = z q ).The air density at AWS glac (ρ a ) is derived as the air density at standard sea-level pressure (ρ 0 =1.29 kg m −3 at 0 • C) multiplied by the ratio between the air pressure estimated at AWS glac (p) and the standard sea-level pressure (p 0 ; 1013 hPa).The dimensionless exchange coefficients for the neutral atmosphere assume logarithmic profile of wind, temperature and humidity with height, and take the following form: C q,log = κ ln zt z0q . (19) 3.2.2K-approach with R ib -based stability corrections (C Rib ) The method assumes the same parameterization for K M,H,E (Eq.( 13)) as C log method, but it allows a flux reduction in a stable stratification via the bulk Richardson number (R ib ) defined as: where temperature is given in Kelvin.For stable conditions (0 < R ib < 0.2) which prevail over a melting glacier surface, the coefficient C for the property y (either v, t or q) is (Webb, 1970): Starting from the parametrization for neutral conditions (Eq. ( 13))) but allowing K M,H,E to vary in response to atmospheric stability (generally, K for statically unstable surface layer > K for neutral > K for statically stable), K M (and equivalently K H and K E ) can be expressed as: ( We use the integrated form of M-O stability functions from Holtslag and de Bruin (1988) and Dyer (1974) as expressed earlier in Section 2.3.The calculation of Ψ requires an estimate of L, which requires an estimate of Q H and u * .This thus becomes a root-finding problem for the three coupled equations, which is either solved by Newton's method or a fixed-point iteration method.We chose the latter since the method is shown to successfully work for this set of equations (e.g., Berkowicz and Prahm, 1982;Lee, 1986) and has been widely used in glacier studies (e.g., Munro, 1989;Conway and Cullen, 2013).In this iterative method each 30-min Q H and u * are initially derived assuming the neutral case ( zv,t L = 0).This approach allows an estimate of L and Ψ, which in the next iteration yields a new estimate for Q H and u * .These steps are repeated until no significant change in Q H occurs (e.g.> ± 1 W m −2 difference), a condition which is satisfied within first five iterations (Munro, 1989).

K-approach with C M −O parameterization and surface renewal model (C SR )
In all the parameterizations for C we use a constant value for roughness lengths derived as the mean log value from the 30min values (Eq.( 5) -( 7)) calculated for the neutral conditions only.An alternative approach in deriving the scalar roughness lengths (z 0s ≡ z 0t,0q ) is through the surface renewal model by Andreas (1987).The model uses a simple similarity arguments considering the structure of the viscous sublayer to derive z 0s from known z 0v as a function of Reynolds roughness number ν , where ν is the kinematic viscosity of air, equal to 1.46 × 10 −5 m 2 s −1 ): where b 0,1,2 are the polynomial coefficients (Table 2).The model has been successfully tested for relatively smooth snow and ice surfaces, while for a rough, hummocky glacier surface Smeets and van den Broeke (2008) found that somewhat different values for the polynomial coefficients (Table 2) give better performance.Here we apply the model for the scalar roughness lengths (Eq.( 24)) and input the modeled values for z 0t,0q into C M −O method, which gives us the fourth bulk scheme in this section: (4) C SR method.

Evaluation results: K-approach with common parameterizations for C
We summarize the evaluation results of the bulk schemes:   1).The modeled fluxes (y-axis in Fig. 5) are the mean values from the Monte Carlo ensemble, whereas the error bars are represented with ± 1 standard deviation (σ) of the ensemble.As illustrated in the scatter plots (Fig. 5), each bulk scheme overestimates the turbulent fluxes (MBE > 0) over the entire period.The scheme with the lowest RMSE for the friction velocity is C Rib , while for the sensible and latent heat fluxes is C SR .The variability in friction velocity (r=0.26,p-value > 0.05) is least successfully modeled, while statistically significant correlations (p < 0.05) are found for both Q H and Q E fluxes, with the highest correlation coefficients found for C log and C Rib method (Table 3).The scatter plots for the 2010 season (not shown) are similar to these for the 2012 season, while the correlation values are smaller in 2010 due to smaller sampling period (Table 3).We also calculated the total model RMSE from the individual 30-min error estimates (error bars in Fig. 5) over the entire period, and found that it is smaller than the RMSE between modeled and observed fluxes for each of the bulk methods.
As illustrated in Fig. 5, the performance of C SR method does not significantly differ from C M −O because the OPEC-derived values for scalar roughness lengths (neutral conditions only) mostly agree with the estimated values from the surface renewal model (Fig. 6).Large roughness Reynolds numbers (Re * ), estimated for our site, indicate a rough flow regime for which the surface renewal model (Andreas, 1987) predicts the scalar roughness lengths two orders of magnitude smaller than z 0v .This prediction agrees well with our observed values, especially when we consider the mean log values of roughness lengths rather than their individual 30-min estimates.Despite the large scatter of 30-min OPEC-derived z0t z0v values around the predicted values from the surface renewal model, more than 50 % of the observed values fall between the predictions of Andreas (1987) and those of Smeets and van den Broeke (2008) derived for hummocky ice with z 0v > 1 × 10 −3 m (Fig. 6).

3.2.6
Evaluation results: K-approach with observed variables (C u * , C z L ) In the gradient-flux relation, used to derive the K-approach bulk methods, the eddy viscosity is parameterized as a function of  variables may influence the performance of the bulk methods.Our goal in this section is to investigate the influence of two components (u * and zv,t L ) on the method performance.Therefore, we assess the turbulent fluxes using the (1)-(4) bulk schemes with OPEC-derived u * and OPEC-derived Obukhov length.The following methods are tested with the OPEC-derived u * only:   1) is used to calculate Re * .
. The findings for Q H from the 2012 season, which has more data points, agree well with the 2010 season: C log u * , and C SR u * methods produced the lowest RMSE and the highest correlation coefficients among all the K-approach methods with observed variables (Table 3).These results suggest that in order to better simulate Q H with the bulk method, getting correct estimates of u * is more important than getting correct estimates of zv L .The variance error (scatter of residuals around 1:1 line) in the C log u * method is smaller than in the C M −O z L method (Fig. 7).The OPEC-derived u * underestimates Q H over the entire period (MBE < 0).Relative to the common K-approach with modeled u * and modeled zv L , the K-approach with the observed variables yields smaller RMSE, especially for C log method where the RMSE is more than halved.The C log u * method outperforms any parameterization with the M-O stability corrections, even in the case when both u * and z . The results for Q E , for both years accord with the results for Q H : using measured instead of modeled u * improves the performance of the bulk method.Despite being less affected by the M-O stability corrections than the modeled sensible heat fluxes are, the modeled latent heat fluxes more closely resemble the observed ones if the OPECderived rather than modeled zv L are used.In summary, C log with measured u * outperforms the other parameterizations but underestimates the sensible heat flux (MBE = -13.4W m −2 for 2012 season, Table 3).In contrast, C log with modeled u * significantly overestimates the sensible heat fluxes (MBE = 23.2W m −2 for the 2012 season, Table 3).

Evaluation results: K-approach with new M-O stability functions (C M −O new)
In the C M −O method with modeled or OPEC-derived (measured) zv,t L , we assumed that the most commonly used M-O stability functions in glacier studies (Holtslag and de Bruin, 1988;Dyer, 1974) are applicable to our site.Our next step is to derive the new integrated stability functions, Ψ v ( zv L ) and Ψ t ( zt L ) to examine how well the optimization of the functions to specifically  fit our data can improve the performance of the C M −O method.To do so, we first apply Eq. ( 5) to derive Ψ v using measured u * , mean log z 0v (Table 1), and measured U zv .We fit a second-order polynomial (Fig. 8) which represents the new integrated stability function for momentum flux, Ψ v ( zv L ).Similarly, we derive Ψ t via Eq.( 6) using measured Θ * , mean log z 0t (Table 1), and measured T zt .The best fit polynomial in the plot, showing a relation between derived Ψ t and zt L (Fig. 8), represents the new integrated stability function for the sensible heat flux, Ψ t ( zt L ).In this fitting exercise we used an additional condition that Ψ v,t ( Due to few data points available during the very stable conditions ( zv,t L > 1), the z-less scaling is assumed following the common approaches for a very stable atmospheric boundary layer (Foken, 2008).Therefore the new stability functions give: .51, and Ψ t ( zt L > 1) = Ψ t ( zt L = 1) =4.50.The goodness-of-fit is calculated by regressing the modeled versus the OPEC-derived Ψ v and Ψ t values, giving statistically significant correlations (p-values < 0.01) with r=0.48 and r=0.35, respectively.The fitted polynomials, i.e. the coefficients in Eq. ( 25) and ( 26) are derived from the 2012 data only and then validated with 2010 data.The statistically significant correlations (p-values < 0.01) between modeled and observed values for Ψ v ( zv L ) and Ψ t ( zt L ) for the validation data (r=0.43 and r=0.41, respectively) demonstrate that the new stability corrections are successfully applicable to both seasons.Errors due to self-correlation (e.g., Klipp and Mahrt, 2004) are probably present in our dataset, however, we assume them negligible since we eliminated the non-stationary turbulence during which these errors can be large.Due to insufficiently large data sample for measured Q E , we do not derive a new stability function for humidity (Ψ q ( zt L )), nor assume that it is equal to the newly derived Ψ t .Instead, we set Ψ q ( zt L ) to zero since the stability seems to play a less important role in the Q E flux modification relative to Q H flux modification, as shown for the K-approach with observed variables (Section 3.2.6).L via the fixed-point iteration method of Munro (1989).As expected, the new stability corrections provide a better match between modeled and observed sensible heat fluxes than in the case with the stability corrections from literature (Fig. 7).If the modeled zv,t L is used, however, the improvement of the bulk method performance is almost negligible.Despite the fact that the fixed-point iteration method successfully converges to a solution, the predicted solutions for

Methods based on katabatic flow models
Here we describe in more details our two types of bulk methods that are rooted in the katabatic flow models, and present the evaluation results: Rather than integrating the flux-gradient relations with a chosen parametrization for K, Oerlemans and Grisogono (2002) derived a bulk approach for surface heat fluxes through a simplified scaling of the governing equations for heat and momentum balance in a 1-D katabatic flow model (Prandtl, 1942).Their basic assumption was that the katabatic flow is characterized by a well defined wind maximum, setting the exchange coefficient for heat proportional to the maximum wind speed and to the height of the wind maximum.With this scaling approach they determined a 'katabatic bulk exchange parameter' for heat flux (C kat , in m s −1 ) as: where k and k 2 are dimensionless empirical constants, Θ s is the temperature deficit at the glacier surface (negative value; , γ is the background potential temperature lapse rate, and P r denotes the eddy Prandtl number Defining the exchange coefficient in this way makes the sensible heat flux to increase quadratically with the surface temperature deficit: Θ s can be replaced by 2 m air-surface temperature difference (∆T ), but it gives only a crude estimate of the actual temperature deficit.The latent heat flux is similarly expressed as: Neither Q H nor Q E explicitly depend on u * because the momentum flux and wind speed are obtained from the katabatic flow model in which the eddy viscosity is parametrized through the background variables (Θ s and γ).C kat becomes smaller as the atmospheric boundary layer becomes more stable (larger γ) because the katabatic flow weakens with increasing γ.As noted in Oerlemans and Grisogono (2002), this is a dynamic effect which should not be confused with the effect of stratification on the exchange coefficients as predicted by the M-O theory.

Evaluation results: C kat method
Our data are not suitable to adequately test this bulk method since we do not have the vertical profile observations of potential temperature needed to estimate Θ s and γ.Nevertheless, our data are sufficient to investigate the validity of the quadratic relation between Q H and ∆T .In particular, assuming that the quadratic fit holds, we will optimize C kat so that the modeled Q H gives the best fit to the observed Q H .Although this optimization can be performed in multiple different ways, we chose to set all the parameters constant except one used for model tuning.The values of constant parameters are taken from the field data from a 10 km long valley glacier in the Austrian Alps (PASTEX-94; Greuell and Struijk, 1994): k 2 =1, γ=0.005K m −1 , P r=2 and T 0 =273 K.The chosen optimization parameter k is then evaluated by minimizing the RMSE between modeled and observed 30-min Q H for each season separately, using the air-surface temperature difference (T zt − T 0 ) as a substitute for Θ s .
The optimization of the k parameter in this method, as our (15) C kat method, yields k = 4.12 × 10 −4 and k = 4.00 × 10 −4 for 2012 and 2012 season, respectively.When this optimized C kat is used to derive sensible and latent heat fluxes (Eq.( 28) and ( 29)) the modeled fluxes simulate the observed ones better (smaller RMSE and MBE) than any of the K-approach methods with common parameterizations (1)-(4) (Table 3).There is a quadratic relation between OPEC-derived Q H and ∆T (Fig. 9), but this relation only explains 58 % of the variance in the observations.Our results do not guarantee the validity of Eq. ( 28), but demonstrate that the simple katabatic model of Oerlemans and Grisogono (2002) reproduces the basic characteristics of katabatic flow at our site; this model can simulate the turbulent heat exchange more successfully than the K-approach method.
Keeping all the parameters, except k, constant in Eq. ( 27) and using the observed 30-min values for Q H and ∆T in Eq. ( 28), we find that k (and therefore C kat ) decreases as the static stability ( zv,t L ) increases, however, the regression between k and zv,t L is not statistically significant.The decrease in C kat as the atmosphere becomes more stable is expected according to the katabatic flow model (Oerlemans and Grisogono, 2002).
Equivalently to the C kat method, the second method assumes that the katabatic flow dominates at the site throughout the observational period.Following Grisogono and Oerlemans (2001) and Parmhed et al. (2004) we use the assumed linear-Gaussian profile for eddy viscosity K(z) as where K max is the maximum value of K(z) (in m 2 s −1 ) reached at the height H K (in m) above the surface.Near the surface, this parameterization for K(z) is similar to O'Brien's cubic polynomial approximation applicable to neutral and stable boundary layer (O'Brien, 1970;Stull, 1988).In the original model of Grisogono and Oerlemans (2001), K max and H K are related to the background variables (γ, Θ s ), slope of the surface, and P r.To adequately determine K max and H K , detailed observations of the wind speed and temperature profiles with height are required.Setting P r=1 (thus setting K M =K H =K E ), inserting the parameterizations for K(z) into the bulk-gradient relations (Eq.( 10) -( 11)) and integrating them with respect to z (treating K max and H K as constants), we obtain the following bulk expressions for the turbulent fluxes: where which is valid for z v < H K .In the absence of vertical profile observations we make use of the OPEC data (i.e.30-min values for u * ), which when inserted in Eq. ( 31) together with observed U zv can give the best estimate for measured K Int (in s m −1 ).
Using the measured K Int in Eq. ( 34), and the log mean value (for neutral conditions only) of z 0v we derive: (a) measured values for K max when H K is set to a constant value of H k =20 m; and (b) measured values for H K when K max = 0.8 m 2 s −1 .The assumed constant values for K max and H K are taken from their frequently observed values from a field study at a glacier in Iceland (Parmhed et al., 2004) as this was the only study we found that directly measured the two parameters.We also investigate whether the measured K max and H K can be expressed as functions of the OPEC-derived zv L .If there is a strong empirical relation, this will indicate that from the local stability metric one can constrain the parameters in the katabatic model.

Evaluation results: K
The relation between the measured 30-min K max (H K ) and the OPEC-derived zv L is strongly non linear (Fig. 10).The best fitting models (Fig. 10) have the following power-law relations: where K max0 =0.22 m 2 s −1 , and H K0 =71.52 m.Note that these empirical fits are assumed valid within the range of observed values for K max and H K taken from Parmhed et al. (2004): 0.03 m 2 s −1 < K max < 2.1 m 2 s −1 , and 5 m < H K < 90 m.The power-law models are optimized on the 2012 data only (Fig. 10), while the 2010 data is used for the model validation.The validation yields high correlation coefficients between the 30-min modeled and observed values: r=0.89 for the K max values, and r=0.83 for the H K values.This result indicates that, for our study site, the same empirical model is equally applicable to both field seasons.
The static stability drives the strength of the proportionality between measured friction velocity and wind speed (Fig. 10), so that the higher the stability, the smaller the constant of proportionality between u * and U zv .The modeled u * derived from Eq. ( 31), where K Int is evaluated using H k =20 m and the K max empirical model (Eq.( 35 this non-linear relation on zv L .For comparison, we also show modeled u * versus U zv where the C log method (K-approach) is used to determine the friction velocity (Eq.( 14)).As expected, for the near-neutral conditions (| zv L | < 0.1), the K-approach works relatively well, because of the near-linear dependence between the OPEC-derived u * and wind speed.However, as soon as the stability increases, the linearity is replaced by a relation that more closely resembles u * ∝ U zv (i.e.momentum flux proportional to wind speed) where the strength of the proportionality depends on the static stability.The same results are found for 2010 season.Finally, we compare the modeled versus observed turbulent fluxes incorporating: (16) K Int K max z L method, which takes the empirical relation for K max (Eq.( 35)) and the OPEC-derived zv L , and ( 17) L method, which takes the empirical relation for H K (Eq.( 36)) and the OPEC-derived zv L .The results show that both methods have almost the same performance, and they give the best estimate of the friction velocity or momentum flux (Fig. 11) among all the bulk schemes we tested so far.The K Int -approach, however, does not perform as well for the sensible and latent heat fluxes as it does for the friction velocity: correlation between observed and modeled fluxes are the lowest among the methods tested (Table 3).We conclude that, while the linear-Gaussian parametrization of K(z) works well for the momentum flux, it works poorly for the heat flux, indicating a shortcoming in the commonly used assumption that K M and K H can share the same parametrization.

Hybrid methods with K
Our final group of bulk methods reflect a 'hybrid' approach where the friction velocity is assessed from the K Int method (Eq.( 34)), while the heat fluxes are calculated from the K-approach C log method (or C M −O new method) following Eq.( 18) and ( 19): In this way, we assume that the linear-Gaussian profile of K(z) dictates the solution for the wind speed profile with height (Grisogono and Oerlemans, 2001), while the sensible and latent turbulent heat fluxes are calculated assuming the validity of the logarithmic vertical profiles for temperature and specific humidity (Prandtl, 1942) respectively, but use the fixed-point iterative method of Munro (1989) to derive zv L .

Evaluation results
The hybrid methods outperform the K Int methods in simulating the Q H fluxes, and outperform the K-approach with the common parameterizations in simulating all the fluxes (Fig. 11).The best performance for Q H fluxes arises from the method, but the application of the new stability corrections only make sense if the OPEC-derived   3).

Discussion and Uncertainty Analysis
We summarize the inter-comparison results, i.e. error estimates, between modeled and observed fluxes, for a selection of the most relevant bulk approaches used in the study (Fig. 12).The error estimates are shown for the 2012 season only, and agree well with the estimates for the 2010 season.By the most relevant bulk approaches we consider those that use OPEC-derived roughness lengths and mean meteorological variables (e.g.2-m temperature, wind speed, specific humidity averaged over 30min time intervals) and modeled or OPEC-derived M-O stability parameter ( zv,t L ).Note that the schemes with OPEC-derived zv,t L do not strictly qualify as the bulk approaches, however, we include them in the analysis in order to distinguish between the actual and the empirical bulk model performance, where the empirical model uses OPEC-derived rather then modeled zv,t L .In implementing the bulk schemes, we assumed that the glacier surface has constant roughness lengths (z 0v,0t,0q ), equal to the mean log values from filtered OPEC data for neutral conditions only (Table 1), throughout the observational period.
The error estimates for each bulk scheme (Fig. 12) are presented as a Mean Square Error (MSE, which is equal to RMSE 2 from Table 3) decomposed into a squared MBE and a variance error (VE), satisfying MSE = MBE 2 + VE.We also show a squared model RMSE calculated from the individual 30-min error estimates (error bars in Fig. 5) over the whole observational period, reflecting the sensitivity of the bulk scheme to the uncertainty in the roughness lengths.The hybrid approach with the K Int method for assessing u * outperforms any other actual bulk model in the study (Fig. 12).The OPEC-derived zv,t L significantly improves the simulations of the turbulent fluxes, in some cases reducing the MSE by more than 50 %.Again, the empirical K Int -approach performs the best in simulating the shear stress among all the empirical bulk approaches, while the empirical hybrid approach (K Int in combination with C M −O new), with the new stability correction Ψ t ( zt L ), performs the best in simulating Q H .For modeled u * , the total error is dominated by the contribution from the mean bias error (MBS 2 > VE), while for the modeled Q H and Q E it is the variance error that dominates.The better performing bulk methods are those that are able to reduce MBE, while VE stays roughly unchanged across the methods.The VE gets reduced only when stability corrections are acting on the OPEC-derived zv,t L , reveling that the remaining variability in the data, which is not explained by 10 the variability in the mean meteorological variables (e.g.U z , ∆T ), is likely explained by the variability in the assumption that γ is inversely proportional to the incoming longwave radiation at the surface (L down ), we performed a correlation analysis between ( zv L ) and L down timeseries.While there is no statistically significant negative correlation for the 30-min values, there is a significant negative correlation when the smoothed timeseries of zv L are used.In other words, the long term (multi-day) fluctuations in the static stability are shown to be negatively correlated with the long term L down fluctuations (results not shown).Considering the successful performance of the K Int -approach based on the assumed linear-Gaussian profile for K(z), a possible extension of our work would be to investigate whether a more exact parametrization for K (as a function of z, ∂U ∂z , and ∂T ∂z ) in the surface layer can be determined as a solution to the governing equations (heat and momentum balance) for a katabatic flow.On the other hand, the poor performance of K Int -approach is assessing the sensible heat fluxes, reveals a shortcoming in the assumption that the eddy diffusivity (K H ) and eddy viscosity (K M ) share the same parametrization.One way to improve for this shortcoming would be to introduce a vertically varying eddy Prandtl number, P r(z).

Conclusions
The main objective of the study was to evaluate commonly used bulk approaches for simulating turbulent heat fluxes at a sloped glacier surface.In particular, we investigated the stability components of the bulk methods, and attempted to improve upon these, assuming site-specific roughness length values are available.In addition to the K-approach with different parameterizations for the bulk exchange coefficient (C), we included a set of less commonly used bulk approaches developed from katabatic flow models (Grisogono and Oerlemans, 2001;Oerlemans and Grisogono, 2002).These 23 bulk schemes, for simulating 30-min turbulent fluxes of momentum, sensible and latent heat, were evaluated against the equivalent fluxes obtained by an open path eddy-covariance (OPEC) method.The evaluation was performed at a point scale of an alpine glacier in BC, using one-level meteorological and OPEC observations from a multi-day period in the 2010 and 2012 summer seasons.Prior to the inter-comparisons, the OPEC data were subjected to a set of quality control corrections and filters.Few 30-min data segments met the criteria for the near-neutral static stability and the steady-state turbulence conditions, revealing a challenge to perform this type of analysis at sites dominated by the stable conditions and drainage (katabatic) flows.Despite the small OPEC data sample available for determining the surface roughness lengths (z 0v,0t,0q ) and the high-quality 30-min turbulent fluxes, the evaluation results derived from the two independent field seasons agree.A summary of our main findings are: Bulk exchange coefficients could be derived from the OPEC-derived roughness lengths for neutral stratification, where the mean log values for z 0v,0t,0q over the entire period are assumed representative of the actual surface roughness and constant in time.In other words, the widely used K-approach facilitated by the standard meteorological observations on glaciers, provides a good approximation of the OPEC-derived 30-min turbulent fluxes during the neutral stability conditions.However, as the stability increases, the K-approach with commonly used parameterizations of the bulk exchange coefficients performs relatively poorly in comparison to the bulk approaches based on a katabatic flow model.
According to the Monin-Obukhov (M-O) stability theory, under stable stratification, the bulk exchange coefficient should be modified by the universal stability functions.However, the stability functions widely used in glacier studies perform The Cryosphere Discuss., doi:10.5194/tc-2017Discuss., doi:10.5194/tc- -80, 2017 Manuscript under review for journal The Cryosphere Discussion started: 31 May 2017 c Author(s) 2017.CC-BY 3.0 License.relatively poorly on our site.The K-approach with commonly used parameterizations for C, with or without the stability corrections, overestimates the sensible heat fluxes (Q H ) over the observational period.The overestimation is due to a positive mean bias in the modeled friction velocities (u * ).Thus, without an adequate model for u * , the use of the stability corrections for Q H will probably fail.Considering that similar findings are derived from a different glacier in this region (Fitzpatrick et al., 2017) these conclusions might apply to other alpine glaciers in similar climatic settings.
The inability of M-O stability functions to adequately correct for the overestimation of u * and Q H during stable conditions was partly resolved by developing the new empirical stability corrections tuned to our data.These new stability functions work well when OPEC-derived zv,t L is used as an input to the bulk method.If OPEC data is not available, however, zv,t L cannot be successfully predicted, and therefore the new stability functions fail to improve the bulk model performance.The proposed fixed-point iterative method for predicting zv,t L , which has been previously argued to perform successfully (Berkowicz and Prahm, 1982;Lee, 1986;Munro, 1989), fails to provide a good approximation of OPEC-derived zv,t L .
Getting a correct estimate of u * , rather than a correct estimate of zv,t L , is shown to be more important for improving the simulation of Q H .We found that, contrary to the K-approach predictions, the relation between the OPEC-derived 30-min u * and observed wind speed (U z ) is not well represented by a linear fit.Our data, instead, reveals u * ∝ √ U z , where the constant of proportionality is a nonlinear function of zv,t L .
The bulk approach derived from the katabatic flow model, abbreviated as K Int approach, outperforms any other bulk scheme in simulating u * , because the K Int approach is rooted in u * ∝ √ U z .The approach performs poorly in simulating Q H , however, revealing a shortcoming in the commonly used assumption that eddy viscosity and diffusivity can use the same parametrization.Applying the K Int approach to assess u * , which is then used in the K-approach with the new M-O stability function to assess Q H gives the best performance across all the bulk methods we tested.In short, Q H can be more successfully modeled with Q H ∝ u * ∆T , where the constant of proportionality is a function of the M-O stability parameter, rather than using the K-approach that relies on Q H ∝ U z ∆T .
The bulk exchange coefficient in the K Int approach is found to be a function of zv,t L , in a very similar way as the bulk exchange coefficient is a function of the background potential temperature lapse rate (γ) in the original model of Grisogono and Oerlemans (2001).As expected by the original model and confirmed by our results, the atmospheric stability increases (larger γ in the model or larger zv,t L in our data) as the bulk exchange coefficient decreases.Provided that γ, taken from observations or regional climate models, is shown to be a good proxy for zv,t L , the K Int approach could be a predictive tool for simulating the turbulent fluxes during prevailing katabatic flows.Clearly more observations are needed to adequately evaluate the model, in particular, multi-level observations of wind speed and temperature within the boundary layer above the glacier surface.
As initially stated, this study provides a foundation towards developing a better parameterization of turbulent fluxes on

Figure 1 .
Figure 1.Map of Castle Creek Glacier with locations of automatic weather station on the glacier (AWS glac ) and two AWS in the glacier vicinity.
are compared to the equivalent OPEC-derived fluxes with the use of standard evaluation metrics: rootmean-square-error (RMSE), mean-bias-error (MBE) and Pearson correlation coefficient (r).First, we derive 30-min turbulent fluxes for each season using the K-approach bulk method with three different parameterization schemes for the bulk exchange coefficients (C): (1) C log -which assumes logarithmic vertical profiles of wind, temperature and specific humidity, (2) C Rib -which assumes stability corrections to the fluxes using R ib , and (3) C M −O -which assumes stability corrections using the universal stability functions of zv,t 8) C SR u * , (9) C M −O z L , (10) C SR z L , (11) C M −0 u * z L , (12) C SR u * z L .In the schemes with C M −O , we apply the most commonly used integrated stability functions in glacier studies.In addition, we derive our own empirical integrated stability functions from the relations between OPEC-derived fluxes and zv,t L , which yields the parameterizations: (13) C M −O new z L and (14) C M −O new.Third, we use bulk methods derived from a simple katabatic flow model: (15) C kat method -following

Figure 3 .
Figure 3. Initial number of observations for estimating the roughness lengths (from 30-min averages of OPEC data) and remaining number of observations after each filter is applied (in the order from top down) for 2010 and 2012 observational periods.Dotted line represents the same but without applying the filter for neutral stability conditions, i.e. all points that satisfy | zv L | < 2 are used in the calculation of roughness lengths.
Fig. 4) by identifying longer periods (several hours to days) with almost stationary down-glacier wind direction (200 • for our site) and relatively high mean wind speeds (> The Cryosphere Discuss., doi:10.5194/tc-2017-80,2017 Manuscript under review for journal The Cryosphere Discussion started: 31 May 2017 c Author(s) 2017.CC-BY 3.0 License.

Figure 4 .
Figure 4. 30-min values of meteorological variables and SEB fluxes measured at the AWS glac in 2010 (left panel) and 2012 (right panel): near surface air temperature (Tz), accumulated precipitation (P ; measured at AWS down ), wind speed (Uz), wind direction (U dir ), incoming (K down ) and reflected shortwave radiation (Kup), OPEC-derived sensible (QH ) and latent (QE) heat flux.Bottom panel shows the OPECderived roughness lengths for momentum (z0v), temperature (z0t) and humidity (z0q), after filtering OPEC data, and their estimated errors.Shaded in blue are the identified time intervals of prevailing katabatic flow (see text).

The
Cryosphere Discuss., doi:10.5194/tc-2017-80,2017 Manuscript under review for journal The Cryosphere Discussion started: 31 May 2017 c Author(s) 2017.CC-BY 3.0 License.whereψ v is a dimensionless shear and L is Obukhov length.Using this parametrization for K and integrating the flux-gradient equations, the expression for C (for the property y) is: (1) C log , (2) C Rib , (3) C M −O , and (4) C SR .First we compare 30min turbulent fluxes of momentum, sensible and latent heat flux, from the 2012 season, derived from the OPEC data (observed fluxes) with the modeled fluxes from each bulk scheme (Fig.5).Only the filtered 30-min segments, except the 'neutrality' filter (Section 2.3), are used for this comparison.To estimate an error for each modeled 30-min flux (error bars in Fig.5) we use a Monte Carlo approach where each bulk method is run 1000 times with randomly perturbed roughness length values.The roughness length values for each Monte Carlo run are picked randomly from an assumed normal distribution of their log The Cryosphere Discuss., doi:10.5194/tc-2017-80,2017 Manuscript under review for journal The Cryosphere Discussion started: 31 May 2017 c Author(s) 2017.CC-BY 3.0 License.Table 2. Values for the polynomial coefficients in Eq. (24) for temperature (bt) and humidity (bq) for three different aerodynamic flow regimes according to (a) Andreas (1987), and (b) Smeets and van den Broeke (2008) Re * ≤ 0.135 0.135 < Re * < 2.5 0.135 ≤ Re * ≤ 1000 0.135 ≤ Re * ≤ 1000 z0v < 10 −3 m z0v > 10 −3 m z, u * and M-O stability parameter ( z L ).If u * and z L are modeled rather than directly measured, the modeling error of the input The Cryosphere Discuss., doi:10.5194/tc-2017-80,2017 Manuscript under review for journal The Cryosphere Discussion started: 31 May 2017 c Author(s) 2017.CC-BY 3.0 License.

Figure 5 .
Figure 5.Comparison of 30-min values of observed (OPEC derived) and modeled (bulk method) shear stress (friction velocity, u * ), sensible (QH ) and latent (QE) heat flux, using four parameterization schemes in the bulk method (see text): C log , C Rib , CM−O, CSR.Dots and error bars show the mean and standard deviation of Monte Carlo ensembles with ∼1000 simulations per point.Red points show values during the neutral conditions (| z v,t L | <0.1) only, and blue points indicate values during the prevailing katabatic flow.Root-mean-square-error (RMSE; W m −2 ), mean bias error (MBE; W m −2 ), and Pearson correlation coefficient (r) are provided for each scatter plot.

Figure 6 .
Figure 6.The ratio of roughness lengths for temperature (z0t) and momentum (z0v) (panel a) and the ratio of roughness lengths for humidity (z0q) and z0v (panel b) versus roughness Reynolds number (Re * ).Model 1 (solid line) is the theoretical prediction ofAndreas (1987), while model 2 (dashed line) is the prediction ofSmeets and van den Broeke (2008).Mean log value of z0v specific to each season (Table1) is used

The
Cryosphere Discuss., doi:10.5194/tc-2017-80,2017 Manuscript under review for journal The Cryosphere Discussion started: 31 May 2017 c Author(s) 2017.CC-BY 3.0 License.Table 3. Results of the comparison between modeled and OPEC-derived sensible (QH ) and latent (QE) heat fluxes, expressed as root-meansquare-error (RMSE), mean bias error (MBE; modeled minus observed), and Pearson correlation coefficient (r) for 2012 and 2010 (values in parenthesis) observational period, given for a range of methods (1) to (23) used to model the fluxes (see text).
zv,t L = 0) = 0.The results yield the following empirical expressions for the new stability corrections applicable to The Cryosphere Discuss., doi:10.5194/tc-2017-80,2017 Manuscript under review for journal The Cryosphere Discussion started: 31 May 2017 c Author(s) 2017.CC-BY 3.0 License.

Figure 7 .
Figure 7.Comparison of 30-min values of observed (OPEC derived) and modeled (bulk method) shear stress (friction velocity, u * ) and sensible heat flux (QH ) using the parameterization schemes in the bulk method (see text): C log u * , CM−O z L , CM−O new z L , and CM−O new.Red points show values during the neutral conditions (| z v,t L | <0.1) only, and blue points indicate values during the prevailing katabatic flow.Root-mean-square-error (RMSE; W m −2 ), mean bias error (MBE; W m −2 ), and Pearson correlation coefficient (r) are provided for each scatter plot.

Figure 8 .
Figure 8. Panels (a) and (b): OPEC-derived 30-min values for stability corrections Ψv and Ψv versus OPEC-derived zv ,t L , and the best fit polynomials representing the new integrated stability functions Ψv( zv L ) and Ψt( z t L ) (solid line), in comparison to the stability functions from the literature (dashed line).Red points show values during the neutral conditions (| z v,t L | <0.1) only, and blue points indicate values during the prevailing katabatic flow.In (c): 30-min modeled z v,t L , via the fixed-point iteration method applied in CM−O, versus 30-min OPEC-derived

Finally
, we test the performance of the C M −O method with the new stability functions, yielding: (13) C M −O new z L method 5 that considers the OPEC-derived zv,t L , and (14) C M −O new method that takes modeled zv,t

L
poorly resemble the OPEC-derived values.The modeled 10 zv L applied in the C M −O and C M −O new methods are significantly lower than the OPEC-derived zv L (Fig.7).Because the static stability zv,t L is consistently underestimated in the bulk methods with the fixed-point iterative scheme, the stability corrections have a very small effect in modifying the modeled turbulent heat fluxes.The Cryosphere Discuss., doi:10.5194/tc-2017-80,2017 Manuscript under review for journal The Cryosphere Discussion started: 31 May 2017 c Author(s) 2017.CC-BY 3.0 License.

Figure 9 .
Figure 9. Observed 30-min sensible heat fluxes (QH ) versus observed air-surface temperature difference (∆T ), and modeled 30-min QH fluxes (black dotted line) according to the C kat method, for (a) 2012 and (b) 2010 season.Each observation point is colored according to its 30-min value of z t L (colorbar).Panels (c) and (d) show the parameter k, derived from 30-min observations of QH and ∆T , versus OPEC-derived z t L for the 2012 and 2010 seasons, respectively.

Figure 10 .
Figure 10.Kmax (panel a) and H k (panel b), calculated from 30-min OPEC-derived u * and observed Uz v , plotted against OPEC-derived zv L for 2012 season.Black circles represent all observations, blues circles are observed values during prevailing katabatic conditions, and red dotted line is the empirical function for Kmax( zv L ) (panel a) and HK ( zv L ) (panel b).Panel (c): OPEC-derived 30-min friction velocity, u * , versus observed wind speed Uz v , for 2012 season, where each observation point is colored according to its 30-min value of zv L (colorbar).Solid lines are modeled u * , using KInt parameterization with H k =20 m and modeled Kmax, for four different stability conditions ( zv L =0.1, 0.5, 1.0, 1.5).Dashed line is modeled u * using the bulk method with C log parameterization.
)), successfully represents two observational features: (i) a non-linear relation between friction velocity and wind speed; and (ii) a functional dependence of The Cryosphere Discuss., doi:10.5194/tc-2017-80,2017 Manuscript under review for journal The Cryosphere Discussion started: 31 May 2017 c Author(s) 2017.CC-BY 3.0 License.

L
is used.Without a priori knowledge of zv,t L , but using the proposed fixed-point iterative method to derive L, the skill in simulating u * and Q H with the The Cryosphere Discuss., doi:10.5194/tc-2017-80,2017 Manuscript under review for journal The Cryosphere Discussion started: 31 May 2017 c Author(s) 2017.CC-BY 3.0 License.

Figure 11 .
Figure 11.Comparison of 30-min values of observed (OPEC derived) and modeled friction velocity (u * ) and sensible heat flux (QH ) using different bulk approaches (see text): C kat , KInt Kmax z L , C log Kmax z L , CM−O new Kmax z L , and C log Kmax.Red points show values during the neutral conditions (| z v,t L | <0.1) only, and blue points indicate values during the prevailing katabatic flow.Root-mean-square-error (RMSE; W m −2 ), mean bias error (MBE; W m −2 ), and Pearson correlation coefficient (r) are provided for each scatter plot.

Figure 12 .
Figure 12.Panel (a): Mean Square Error (MSE), squared Mean Bias Error (MBE 2 ), Variance Error (VE) and squared model RMSE (δ 2 ; derived from the Monte Carlo runs with perturbed roughness lengths -see text) for a selection of most relevant bulk approaches in the study, for 2012 season.Note that MSE = MBE 2 + VE.Panels (b) and (c): OPEC-derived 30-min sensible heat flux (QH ) versus measured products u * ∆T and Uz v ∆T , respectively, for 2012 season.Each point is colored according to its 30-min OPEC-derived zv L .Panels (d) and (e): Measured 30-min values of u * Uz v and u * √ Uz v, respectively, versus OPEC-derived zv L for the two years.
assessment of Q H with the C M −O methods has the largest sensitivity to the perturbations in the roughness lengths, probably because the uncertainties in both z 0v and z 0t propagate into the modeled zv,t L and then, via the fixed-point iterative method, back-propagate into the Q H estimates.The best overall performing model (C M −O new K max z L ) for Q H has the highest sensitivity to the The Cryosphere Discuss., doi:10.5194/tc-2017-80,2017 Manuscript under review for journal The Cryosphere Discussion started: 31 May 2017 c Author(s) 2017.CC-BY 3.0 License.
sloping glacier surfaces subject to katabatic flows.OPEC data from an expanded SEB monitoring program on glaciers in the The Cryosphere Discuss., doi:10.5194/tc-2017-80,2017 Manuscript under review for journal The Cryosphere Discussion started: 31 May 2017 c Author(s) 2017.CC-BY 3.0 License.region, commenced in 2014, will allow for further analysis and help us to address the outstanding questions that arose in this study.One of the biggest challenges, as we see it, for both measuring and parameterizing the turbulent heat fluxes, is the presence of the intermittent turbulence that often accompanies katabatic flows.To tackle the challenge from the measurement side one would need to use the detection methods of intermittency through a spectrum and/or wavelet analysis on the eddycovariance data.On the modeling side, one would need to invest into new ways of parameterizing the eddy viscosity that go The Cryosphere Discuss., doi:10.5194/tc-2017-80,2017 Manuscript under review for journal The Cryosphere Discussion started: 31 May 2017 c Author(s) 2017.CC-BY 3.0 License.The Cryosphere Discuss., doi:10.5194/tc-2017-80,2017 Manuscript under review for journal The Cryosphere Discussion started: 31 May 2017 c Author(s) 2017.CC-BY 3.0 License.

Table 1 .
Mean and ± standard deviation (σ) of log roughness lengths of momentum (z0v), temperature (z0t) and humidity (z0q) derived for neutral stability conditions (| z L | < 0.1) during 2010 and 2012 observational period.Also shown is weighted mean ± weighted standard deviation (σw; see text) and number of points (N ) used to assess these values, i.e. number of points remaining after all the filters have been applied to the original data.
method with the new stability functionΨ t ( zt L ) in the parameterization for C M −O,t .The final two are (22) C log K max method and (23) C log H K method, which share the same formulation as (18) and (19),