Spatial and temporal variations in glacier aerodynamic surface roughness during the melting season, as estimated at the August-one ice cap, Qilian mountains, China

The aerodynamic roughness of glacier surfaces is an important factor governing turbulent heat transfer. Previous studies rarely estimated spatial and temporal variation in aerodynamic surface roughness (z0) over a whole glacier and whole melting season. Such observations can do much to help us understand variation in z0 and thus variations in turbulent heat transfer. This study, at the August-one ice cap in the Qilian mountains, collected three-dimensional ice surface data at plot scale, using both automatic and manual close-range digital photogrammetry. Data were collected from sampling sites spanning the whole ice cap for the whole of the melting season. The automatic site collected daily photogrammetric measurements from July to September of 2018 for a plot near the center of the ice cap. During this time, snow cover gave way to ice and then returned to snow. z0 was estimated based on micro-topographic methods from automatic and manual photogrammetric data. Manual measurements were taken at sites from the terminals to the top of the ice cap; they showed that z0 was larger at the snow and ice transition zone than in areas that are fully snow or ice covered. This zone moved up the ice cap during the melting season. It is clear that persistent snowfall and rainfall both reduce z0. Using data from a meteorological station near the automatic photogrammetry site, we were able to calculate surface energy balances over the course of the melting season. We found that high or rising turbulent heat, as a component of surface energy balance, tended to produce a smooth ice surface and a smaller z0 and that low or decreasing turbulent heat tended to produce a rougher surface and larger z0.


Introduction
The roughness of ice surfaces is an important control on air-ice heat transfer, on the ice surface albedo, and thus on the surface energy balance (Greuell and Smeets, 2001;Hock and Holmgren, 2005;Irvine-Fynn et al., 2014;Steiner et al., 2018). The snow and ice surface roughness at centimeter and millimeter scales is also an important parameter in studies of wind transport, snowdrifts, snowfall, snow grain size, and ice surface melt (Bruce and Smeets, 2000;Brock et al., 2006;McClung and Schaerer, 2006;Fassnacht et al., 2009a;Fassnacht et al., 2009b). Radar sensor signals, such as SAR (Oveisgharan and Zebker, 2007), 35 altimeters and scatter meters, are also affected by ice and snow surface roughness (Lacroix et al., 2007;Lacroix et al., 2008).
One of the most important of these influences is the aerodynamic roughness of z0, which is related to ice surface topographic roughness in a complex way (Andreas, 2002;Lehning et al., 2002;Smith et al., 2014;Smith et al., 2016). Determination of z0 based on topographic roughness is therefore of great interest for energy-balance studies (Greuell and Smeets, 2001).
Glacier surface z0 has been widely studied through such methods as eddy covariance (Munro, 1989;Smeets et al., 2000;Smeets 40 and Van den Broeke, 2008;Fitzpatrick et al., 2019), or wind profile (Wendler and Streten, 1969;Greuell and Smeets, 2001;Denby and Snellen, 2002;Miles et al., 2017;Quincey et al., 2017). However, micro-topographic estimated z0 shows some advantages, such as lower scatter, rather than profile measurements over slush and ice (Brock et al., 2006), and ease of application at different locations (Smith et al., 2016). Current research has increasingly used micro-topographic method to estimate z0. It has also become clear that it is important to estimate z0 over the entire course of the melting season and at many 45 points on the glacier surface, as z0 is prone to large spatial and temporal variation (Brock et al., 2006;Smeets and Van den Broeke, 2008). This variation is due to variations in weather and snowfall (Albert and Hawley, 2002). The micro-topographic estimated z0 allows repeated measurement at many points on the glacier surface, which is not possible with wind profile or eddy covariance methods.
Photogrammetry has been increasingly popular as a method to measure the aerodynamic surface roughness of snow and 50 ice (Irvine-Fynn et al., 2014;Smith et al., 2016;Miles et al., 2017;Quincey et al., 2017;Fitzpatrick et al., 2019).. Initially, the Micro-topographic method was developed as snow digital photos were taken against a dark background plate. The contrast between the surface photo and the plate could then be quantified as a measure of glacier roughness (Rees, 1998). This methods still widely applied to quantify glacier surface roughness (Rees and Arnold, 2006;Fassnacht et al., 2009a;Fassnacht et al., 2009b;Manninen et al., 2012). A more recent method, as described by Irvine-Fynn et al. (2014), uses modern consumer-grade 55 digital cameras to do close-range photogrammetry at plot scale (small plots of only a few square meters). Appropriate image settings and acquisition geometry allow the collection of high-resolution data (Irvine-Fynn et al., 2014;Rounce et al., 2015;Smith et al., 2016;Miles et al., 2017;Quincey et al., 2017). Such data facilitates the distributed parameterization of aerodynamic surface roughness over glacier surfaces (Smith et al., 2016;Miles et al., 2017;Fitzpatrick et al., 2019). Precision of microtopographic estimated z0 also became a major concern, and many comparative studies with the aerodynamic method 60 (eddy covariance or wind towers measurements) were carried out over debris-covered or non-debris covered glaciers. The difference was within an order of magnitude for some studies (Fitzpatrick et al., 2019) or strongly correlated (Miles et al., 2017).
Previous researchers have performed some long-term, systematic studies of glacier surfaces (Smeets et al., 1999;Brock et al., 2006;Smeets and Van den Broeke, 2008;Smith et al., 2016). The current study applied such methods to the study of snow and ice aerodynamic surface roughness during melting season at the August-one ice cap. We used both automatic digital photogrammetry and manual photogrammetry. Automatic methods allowed us to monitor daily variations in aerodynamic surface roughness; manual methods allowed us to characterize aerodynamic surface roughness variation along the main glacial flow line. We also recorded meteorological observations, so as to study the impact of weather conditions (e.g. snowfall or rainfall) on aerodynamic surface roughness. This data allowed a further effort to characterize variation of plot-scale z0 from an 70 energy balance perspective.

Study area and meteorological data
The August-one glacier ice cap is located in the middle of Qilian Mountains on the northeastern edge of the Tibetan Plateau ( Figure 1a, 1b). The glacier is a flat-topped ice cap that is approximately 2.3 km long and 2.4 km 2 in area. It ranges in elevation 80 from 4550 to 4820m a.s.l. (Guo et al., 2015). This study was conducted during the melting season of 2018, a season characterized by high precipitation. Energy balance analysis indicated that net radiation contribute 86% and turbulent heat fluxes contribute about 14% to the energy budget in the melting season. A sustained period of positive turbulent latent flux exists on the August-one ice cap in August, causing faster melt rate in this period (Qing et al., 2018).

85
Researchers had access to meteorological data that had been recorded continuously since September 2015, when an automatic weather station (AWS) was sited at the top of the ice cap (Table 1). The AWS measures air temperature, relative humidity, and wind speed at 2 and 4 m above the surface. Air pressure, incoming and reflected solar radiation, incoming and outgoing long wave radiation, glacial surface temperature (using an infrared thermometer) are measured at 2 m height. Mass balance is measured by a Campbell Scientific ultrasonic depth gauge (UDG) close to the AWS. An all-weather precipitation gauge 90 adjacent to the AWS measures solid and liquid precipitation. All sensors sample data every 15 seconds. Half-hourly means are stored on a data logger (CR1000, Campbell, USA). Throughout the entire melting season (from June to September) researchers periodically checked the AWS station, to make sure that it remained horizontal and in good working order. During the entire study period, precipitation total was 261.3mm as measured at the AWS. Of that, 172.1 mm was snow or sleet and 89.2 mm was rainfall ( Figure.   This was done on July 10, 2018. A wooden frame, 1.5 m wide, and 2 m long, was put on the ice surface. This frame served as a geo-reference control field ( Figure. 3a). Four feature points demarcated the control field; three additional points served as 105 check points. A Canon EOS 1300D cameras with an image size of 5184×3456 pixels was connected to the frame. The camera lens was set in wide-angle mode (focal length of 27mm). The f-stop was fixed at f 25 with an exposure time of 1/320s. The camera was programmed to automatically take seven pictures over a period of ten minutes. The photography was repeated at three-hour intervals from 9:00 AM to 18:00 AM, Beijing time. During the ten-minute photography periods, the camera moved along a 1.5 m long slider rail. The camera was 1.7m above ice surface and moved along the control frame. The seven pictures 110 taken during this period were merged to produce a picture of ice surface topography at millimeter scale ( Figure 3b). This apparatus took pictures over a period of three months (July 12 to September 15, the melting season). Sixty-four days of data were recorded. Each daily photography series produced four sets of pictures (twelve hours, three hour intervals). The bestexposed photo sets were manually selected and used as that day's data. We also set up instrumentation to record incoming and reflected solar radiation. Samples were taken every 15 seconds; 10-minute means were stored on a data logger (CR800,

115
Campbell, USA) located at a height of 1.5m. Surface elevation changes caused by accumulation and ablation was measured by a digital infrared hunting-video camera, which took pictures of ice-surface gauge stakes located near the automatic photogrammetry site.

Manual photogrammetry
Manual close-range photogrammetry was used to survey glacier surfaces at several different locations of the ice cap.
Observations were made on four days: July 12 and 25, and later on August 3 and 28. It should be noted that when the July Channels account for only a small portion of the glacier surface area. These surfaces show extreme variability of z0 (Rippin et al., 2015;Smith et al., 2016). For that reason, we distributed the manual photogrammetry study sites over the glacier surface in such a way as to cover most surface types and topographic regions without including any channels ( Figure 1b). We photographed a total of thirty-six sites over the four days of observation.

130
Study plots were demarcated with a 1.1×1.1m portable square aluminum frame. Geo-reference of the point cloud was enabled using control points established by eight cross-shaped screws on the aluminum frame ( Figure.

Data processing
Structure-from-motion photogrammetry is revolutionizing the collection of detailed topographic data (Westoby et al., 2012;James et al., 2017). High resolution DEMs produced from photographs acquired with consumer cameras need careful handling (James and Robson, 2014). In this study, both manual and automatically derived photographs were imported into a software 150 program, Agisoft Photoscan Professional 1.4.0. This software allowed us to estimate camera intrinsic parameters, camera positions, and scene geometry. Agisoft Photoscan Professional is a commercial package which implements all stages of photogrammetric processing (James et al., 2017). It has previously been used to generate three-dimensional point clouds and digital elevation models of debris-covered glaciers (Miles et al., 2017;Quincey et al., 2017;Steiner et al., 2018), ice surfaces and braided meltwater rivers (Javernick et al., 2014;Smith et al., 2016). In our study, we found that after new snowfall, it was 155 difficult to match feature points in the photo sets. Three days of automatic data could not be processed. We estimated z0 data for the missing days based on data from snowfall days at the automatic site.

160
Methods for measuring roughness at plot-scale were first developed by soil scientists (Dong et al., 1992;Smith, 2014). Metrics such as the random roughness (RR) or root mean square height deviation (σ), the sum of the absolute slopes (ΣS), the microrelief index (MI), and the peak frequency (the number of elevation peaks per unit transect length) were used. Later these roughness indices were used to describe snow or ice surface roughness (Rees and Arnold, 2006;Fassnacht et al., 2009b;Irvine-165 Fynn et al., 2014).
Current photogrammetry methods produce high-resolution three-dimensional topographic data. Earlier two-dimensional profile-based methods for estimating surface roughness discard much of the potentially useful three-dimensional topographic data (Passalacqua et al., 2015). Smith et al. (2016) were able to use equation (1), developed by Lettau (1969), to make better use of the topographic data, using multiple point clouds and digital elevation models (DEM). In this method, z0 is quantified as: where: h * represents the effective obstacle height (m) and is calculated as the average vertical extent of micro-topographic variations; s is the silhouette area facing upwind (m 2 ); S is the unit ground area occupied by micro-topographic obstacles (m 2 ); and 0.5 is an averaged drag coefficient.
Based on the work of Lettau (1969), Munro (1989) simplified the equation (1) by assuming that h* can equal twice the standard deviation of elevations in the de-trended profile, with the profile's mean elevation set to 0 meter. The aerodynamic roughness 180 length for a given profile then becomes Where f is the number of up-crossings above the mean elevation in profile; X is the length (m) of profile, and σd is the standard derivation of elevations of profile.For manual photogrammetry, we put the aluminum frame horizontally over the ice surface, the plot is detrended by setting the control points at z axis of the same values. For automatic photogrammetry, the control field 185 of wooden frame was also laid horizontally over the ice surface that lowered as the ice melted and maintained a horizontal position between the control field and ice surface. A DEM based approach enables the roughness frontal area s to be calculated directly for each cardinal wind direction (Smith et al., 2016). The combined roughness frontal area was calculated across the plot, the ground area occupied by micro-topographic obstacles is 1m 2 . We used a DEM-based average ( ̅ 0_ ) of four cardinal wind directions to represent overall aerodynamic surface roughness. Based on the half-hour wind direction data at the 190 August-one ice cap, the daily upward wind direction DEM-based z0_DEM was also estimated at the automatic photogrammetry site. Considering that wind direction changed during the day, in this case we selected the prevailing wind direction to calculate frontal area s. The prevailing upwind direction DEM-based z0_DEM was applied to calculate turbulent heat flux. Using the Munro (1989) method, z0_Profile was calculated for every profile (n=1000) in both orthogonal directions for each plot at the automatic photogrammetry site.

Snow and ice surface energy balance calculation
The temporal variation of z0 at the automatic site was studied from energy balance perspective. The surface heat balance of a melting glacier is given by: Where, QM is the heat flux of melting; Qis is the incoming shortwave radiation; Qos is the outgoing shortwave radiation; QL is the net longwave radiation; QE is the latent heat flux; QH is the sensible heat flux; QP is the heat from rain; and QG is subsurface heat flux.

205
In a horizontally homogeneous and steady surface state, the surface heat fluxes QE and QH can be calculated using either the bulk aerodynamic approach or profile method, based on the Monin-Obukhov similarity theory (e.g., ; Arck and Scherer, 2002;Garratt, 1992;Oke, 1987). In this study, half-hour observations at 4 m level and daily upward wind direction DEM-based z0 were used to calculate QE and QH based on the bulk method. The heat from rain is given by Konya and Matsumoto (2010): Where, ρw is the density of water(1000 kg m -3 ); CW is the specific heat of water (4187.6 J kg -1 K -1 ); TW is the wet-bulb temperature (K); and Pr is the rainfall intensity (mm). The subsurface heat flux QG is estimated from the from the temperaturedepth profile and is given by where kT is the thermal conductivity, 0.4Wm -1 K -1 for old snow and 2.2W m -1 K -1 for pure ice (Oke, 1987).

215
In order to calculate Pr, we used the air temperatures recorded at the AWS. There is an elevation difference between the study site (4700 m) and the AWS (4790m); recorded air temperatures were corrected to account for the elevation difference, a lapse rate of -5.6 o C Km -1 was applied based on observation nearby (Chen et al., 2014) . The ice cap is flat and open terrain so in this case wind speed and relative humidity at the study sites were assumed to be close to those observed at the AWS.

Photogrammetry precision
We used seventeen plots to analyze the horizontal and vertical accuracy of our automatic photogrammetry, and thirty-one plots for our manual photogrammetry. Based on the Agisoft PhotoScan processing report, automatic photogrammetry average point 225 density of the final plot point clouds was over 1,000,000 points m -2 . DEMs of 1mm resolution were generated at plot scale.
The average geo-reference errors fluctuated at around 1 millimeter (see Tables 2 and 3). Total RMSE of the automatic control points was 3.0±2.1 mm, for check points 3.62±1.6 mm. Vertical error for control points was 3.58mm±3.01mm, and 4.83 ±2.9mm for check points (Tables 2 and 3 (Table 1).
Control points vertical accuracy of manual photogrammetry is about 1.65±1.3 mm. Total RMSE of manual photogrammetry check points is 0.99±0.3 mm, vertical accuracy is 0.66±0.3mm (see Tables 2 and3). Standard deviation for x, y and z axis were all within 5mm (Figure 4 b, 4d, 4f).  Note that the control and check point errors are larger for the automatic measurements than for the manual ones (See Figures   4). We believe that this is the case because, rather than using static f-stop and exposure times (as in automatic photogrammetry) 240 researchers engaged in manual photogrammetry could adjust exposure time based on ice surface conditions. This allowed production of better quality photos even on cloudy or foggy days. The difference of survey design also caused more precise results for manual than automatic photogrammetry. For the automatic measurements, the camera was moving linearly, and the density of tie-points was much higher in the foreground compared to the background. For the manual method, photos were taken by surrounding the target area. This type of surface provided a much more robust elevation model and points density.

Aerodynamic surface roughness as measured by automatic photogrammetry
Data for ice surface roughness was collected by the automatic photogrammetry camera site from July 12 to September 15, a period covered the whole melting season. Profile and DEM data show that z0 estimates vary by two orders of magnitude 255 over the study period ( Figure 5). The upwind DEM-based data showed a z0_DEM varying from 0.   was highest in the transition zone between snow cover and ice. This zone moved up the ice cap during the melting season. On July 12, ice surface roughness decreased from 3.2mm to 0.25mm as altitude 285 increased ( Figure. 6a, r= 0.8429, P=0.0006<0.01). Near the ice cap terminals of 4590m, the ice surface featured porous snow/ice and many cryoconite holes. As altitude increased, the number of cryoconite holes decreased and snow coverage increased. At 4700m the ice surface was predominantly snow covered, and only a few small patches were bare of snow. On July 25, ice surface roughness fluctuated between 0.27 to 0.65 mm at the ice cap terminals (4593m). At ~4700m, roughness increased to 1.85mm. Above that point, roughness gradually decreased to 0.25mm at the ice cap top, which was covered by 290 snow (Figure 6b).
On August 3, the August-one ice cap was predominantly bare ice; there was scattered snow crust at the ice cap top. The ice surface, (terminals to top) showed a heavy deposit of cryoconite. Potogrammetric data collected manually revealed that ice surface roughness increased with altitude ( Figure. 6c, r=0.7). From terminals to top, z0 varied from 0.06 mm to 2.2 mm. On August 29, the ice cap surface roughness showed no significant correlation with altitude ( Figure.

Z0 and weather
Figure 7 compared ̅ 0_ and corresponding meteorological conditions of precipitation, air temperature, downward solar 305 radiation, relative humidity and wind speed. Detailed analysis indicates snowfall was recorded from July 12 to 24. In general, snowfall reduced roughness if it resulted in a fully snow-covered surface. However, if a patchy, shallow snow cover was formed, it tended to increase z0 after a short drop. For example, on August 11 and 12, two successive sleety days created a patchy snow cover which soon increased z0. Between July 26 and August 31 there were sixteen rainfall events, which tended to lower ice surface z0.

310
Daily temperatures during the study period ranged from -6.5 °C to 7.1 °C (mean: 1.3, Figure. 7c). It was 1.2 °C on July 11. It increased to 3.6 °C on July 24 (the date when z0 was highest). It continued increasing until July 29, when it reached its highest Daily downward mean solar radiation fluctuated dramatically during the study period due to cloud and overcast ( Figure. 7d).
Incident solar radiation fluctuated between 129W m -2 and 753 W m -2 (mean: 469 W m -2 ). From July 29 to end of August, the weather was cloudy, warm, calm, and humid most of the time ( Figure. 7b, 7c, 7e 7f), and ̅ 0_ was relatively stable except when there was intermittent snowfall-induced fluctuation. After in September, the weather was again becoming cold and dry 320 and z0 was quite variable.

Ice-surface energy balance at automatic z0 observation study site
The following section analyzes the changes in surface energy balance at the automatic site. Meteorological  for -28% to 32% (mean: 15%) of the net energy flux. Latent heat was generally small throughout the study period. Daily mean of latent heat varied from -80.1 to 11.1 W m -2 (mean: -13.2 W m -2 ). It account for a mere 0.9% for the total incoming flux. It was negative from July 11 to 26 when the ice surface was snow covered. After July 26 the latent heat was mainly positive in From July 25 to August 5 rainfall energy varied from 0 to 11.7 W m -2 (mean: 0.3W m -2 ). Rainfall accounted for a mere 0.2% of total incoming flux. One event accounted for much of the total: on July 28 a 31mm rainfall event added a flux of 11.7 W m -2 , which resulted in visible smoothing of the ice surface (Figure 9). Compared to other energy components, QG was very small,

350
with a daily mean of -0.65 W m -2 and a maximum and minimum of -0.4 and -2.1 W m -2 , respectively.

355
Based on the previously listed measurements of energy fluxes we calculated the probable surface ablation at the automatic photogrammetry site. We took into account observed net radiation, bulk method calculated turbulent heat fluxes, heat from rainfall, and subsurface heat flux. There was good agreement between the model and observed results ( Figure 10).

390
Note that aerodynamic surface roughness on days when snow fell was strongly affected by the amount of the snowfall. If we exclude snowfall days and snow covered period, we see a significant exponential relationship between ice surface z0_DEM and LS (Figure 12a, r= -0.34). Scatter diagrams showed significant exponential relationship between ice surface z0_Profile and LS and net longwave radiation (Figure12b, r=-0.69). ̅ 0_ vs. LS also showed a significant exponential relationship (Figure 12c, r=-0.46). Scatter diagrams in Figure 12 also showed z0 did not keep decreasing when LS was above 0.2. z0_DEM, z0_Profile and 395 ̅ 0_ was around 0.56±0.21mm, 0.33±0.03 mm and 0.6±0.26 mm, respectively.
The z0 (z0_DEM, z0_Profile ̅ 0_ ) vs. LS graph indicates that when turbulence and rainfall heat increased, aerodynamic surface roughness decreased. As soon as LS is above 0.2, the ice surface will not keep smoothing and z0 sustained its lowest stage.
Time series correlation of all main energy items and z0_Profile were performed. Table 4 shows an example of the lagged correlations between z0_profile and five variables. The z0 and net shortwave radiation displayed a positive correlation with 0 to 400 1 days lag time. The z0 response to QE with a correlation of -0.6 showed a lag of 0 to 1 days. The z0_Profile also had a negative relationship with QL with no lag or 1 day lag time. The z0_Profile response to LS with a correlation of -0.58 was with a lag of 0 to 2 days. 0 to 2 days lag time gives an indication of the main energy items efforts limitations over ice surface z0. In other words, a sunny and cold day facilitates rough ice surfaces; warm and cloudy days tend to produce a smoother ice surface.
When net shortwave radiation is higher, and if latent and sensible heat were smaller, z0 would tend to be higher for the next 2 405 days. When net shortwave radiation is smaller, as on cloudy days, any snowfall or rainfall is usually associated with smaller z0 for the following 2 days. Under a negative QM, the surface z0 would be not affected by melting process.

420
Photogrammetric techniques such as Structure from Motion (SfM) (James and Robson, 2012) and Multi-view Stereo (MVS) represent a low-cost option for acquiring high-resolution topographic data. Such approaches require relatively little training and are extremely inexpensive (Westoby et al., 2012;Fonstad et al., 2013;Passalacqua et al., 2015). We used both automatic and manual photogrammetric methods to sample spatial and temporal z0 variation at the August-one ice cap. Adjustments to 425 exposure time based on ice surface conditions and survey design of the area surrounding the target made the manual photogrammetry is more precise than automatic photogrammetry (Tables 2 and 3). However, precision is not always the major concern. The glacier surface was a harsh, even punishing environment for the researchers doing manual photogrammetry. In addition, manual photogrammetry took much longer. Automatic methods reduced hours of field work, spared researchers, and produced nearly continuous data. Cloudy or frosty weather affected automatic photogrammetry exposures, and heavy snowfalls 430 resulted in a texture-less surface. Nevertheless, it is likely that photogrammetry techniques will continue to improve and that these drawbacks may be mitigated.

435
Previous studies of glacier surfaces roughness have rarely covered the whole glacier, from terminals to top, in one melting season (Föhn, 1973;Smeets et al., 1999;Denby and Smeets, 2000;Greuell and Smeets, 2001;Albert and Hawley, 2002;Brock et al., 2006;Smeets and Van den Broeke, 2008;Smith et al., 2016). This whole-glacier study allowed us to follow the movement of the transition zone, where snow was melting and exposing ice, from terminals to top. The transition zone moved up as the melting season proceeded, so roughening the surface of the glacier and raising z0.At the start of the melting season, snow cover 440 first disappeared, leaving an ice surface, at the terminals end of the August-one ice cap, that is, at the lower altitude. This newly exposed surface was rougher (z0 was higher) than on the upper part of glacier, which was still snow covered (see the black line Figure 6a for z0 distribution at different altitudes). As the snowline shifted to higher altitudes, ice surface increased, as did z0 (see the dashed black curve in Figure 6b). As the melting continued, the snow and ice transition belt reached the top of glacier (see the dotted curve in Figure 6c). When the ice cap was completely free of snow, z0 and elevation were no longer correlated 445 (see the dotted-dashed line in Figure 6d). In summary, maximum z0 was recorded at the cross-glacier transition zone between snow and ice. This zone shifted from lower altitude to higher altitude, from terminals to top, during the melting season. The spatial pattern of z0 distribution affected turbulent fluxes. The transition zone had maximum z0 and the zone also migrated across much of the glacier, highlighting the importance of transient surface characteristics.
Micro-topography, wind profile, and eddy covariance methods generate a wide range of z0 values for snow and ice surfaces 450 (Grainger and Lister, 1966;Munro, 1989;Bintanja and Broeke et al., 1995;Schneider, 1999;Hock and Holmgren, 2005;Brock et al., 2006;Andreas et al., 2010;Gromke et al., 2011)., In this study, z0_profile, z0_DEM, and ̅ 0_ showed similar variation pattern during the melting season. The difference of z0_profile, z0_DEM, and ̅ 0_ were within one order of magnitude. The latent and sensible heat calculated by z0_profile, z0_DEM, and ̅ 0_ were highly relevant among these methods. The automatic photogrammetry estimated z0 for snow-covered surfaces ranged from 0.1 to 0.55. New snowfall at snow surface in July formed 455 the lowest z0 values. Previous studies have shown that freshly fallen snow is subject to rapid destructive metamorphism (McClung and Schaerer, 2006), which can dramatically change the roughness of fresh snow surfaces (Fassnacht et al., 2009b).
Our study showed that z0 followed an increasing trend during melting season. Intermittent snowfall first decreased snow surface z0, which then began to increase as the snow surface deteriorated. In the data from Clifton et al. (2008), snow surface z0 was estimated at between 0.17 to 0.6 mm in a wind tunnel experiment. In an analysis of ultra-sonic anemometer recorder 460 data over snow-covered sea-ice, Andreas et al. (2010) found z0 values ranging from 10 -2 to 10 1 mm. In a wind-tunnel For patchy snow-covered ice surfaces, z0 varied from 0.5 to 2.6mm and ice surface z0 varied from 0.24 to 1.1mm. During the melting season, there were no blowing snow events and snow surface z0 was relatively smaller than in patchy snow-covered surface or ice surface. Ice surface z0 was generally larger than snow surface and smaller than patch snow-covered surface. Our

470
results match values reported in studies reporting results ranging from. 0.1mm to 6.9mm in Qilian mountain glaciers (Guo et al., 2018;Sun et al., 2018). Our results showed that z0 reached its maximum at the end of the summer melt, which matched wind profile measurements by Smeets and Broeke (2008).
The aerodynamic surface roughness is influenced by both boundary layer and the surface. In this study, the microtopographic estimated aerodynamic surface roughness only considers surface topography at plot scale, but its variability influenced by its 475 surrounding topography and boundary layer. Thus, the results of z0 estimated in this study still need validated by wind tower or eddy covariance observations. However, microtopographic roughness metrics are a very strong proxy for z0 (e.g. Nield et al, 2013), so we have much more confidence in the temporal and spatial variability presented by this work.
Surface geometry roughness develops due to local melt inhomogeneities in melting season. In early work, researchers argued that a variety of ablation forms, such as sun cups, penitents, cryoconite holes or dirt cones are formed by the sun (Matthes, 1934;Lliboutry, 1954;Mcintyre, 1984;Rhodes et al., 1987;Betterton, 2000). These ablation forms develop in regions with bright sunlight and cold, dry weather conditions are apparently required (Rhodes et al., 1987). These structures are observed 485 to decay if the weather is cloudy or very windy (Matthes, 1934;Lliboutry, 1954;Mcintyre, 1984).
The August-one ice cap dust concentrations are high in the melting season. Cryoconites are unevenly distributed over the ice surface leading to differential absorption of shortwave radiation at microscale. This process results in the roughening of the ice surface; a process that enhances turbulent heat exchange across the atmospheric boundary layer-ice interface. When the air temperature is above 0 o C, the ice surface keeps melting. The turbulent heat smooths the ice surface and increases the cryoconite 490 concentration over the ice surface and decreases ice surface albedo, enhancing shortwave radiation absorption (Figure 9). This roughening and smoothing process makes ice surface z0 to fluctuate at around 0.56 mm as long as the air temperature is above 0 o C. When temperature drops below 0 o C, bright sunlight and dry weather shutdown the ice surface smoothing process. The shortwave radiation induces even rougher ice and larger z0 until snow covers the ice surface. At the August-one ice cap, the turbulent heat contributes a small portion of incoming energy, but the smoothing ice surface process decreases ice surface 495 albedo and seems enhance ice surface shortwave radiation. The z0 fluctuation in the melt season is similar with cryoconite holes developing when the radiative flux is dominant and decaying when turbulent heat is dominant (McIntyre, 1984;Takeuchi et al., 2018). The glacier surface energy balance components vs. z0 analysis in this study confirms that main energy items of net shortwave radiation and turbulent heat flux affect the same day and following 2 days z0. This study found an exponential relationship between z0 and LS. The delicate role of z0 played in the ice surface balance is still not fully known. Further 500 comparative studies are needed to investigate the z0 variation through eddy covariance, profile method and DEM-based z0 estimation.

505
Manual and automatic measurements of snow and ice surface roughness at the August-one ice cap showed spatial and temporal variation in z0 over the melting season. Manual measurements, taken from terminals to glacier top, show that the nature of the surface cover features are correlated with z0 rank in this order: transition region > pure ice area or pure snow area. The transition region forms a zone of maximum z0, which shifts, over the melting season, from terminals to top. The observed z0 vs energy items analysis indicated that LS (turbulent heat index) was also an important determinant of ice aerodynamic surface roughness.

510
Aerodynamic surface roughness is a major parameter in calculations of glacier-surface turbulent heat fluxes. In previous studies investigators used a constant z0 value for the whole surface of the glacier. This study captures much smaller scale variation spatial and temporal glacier surface aerodynamic roughness through automatic and manual photogrammetric observations.
Such close observation of variation in z0 certainly enhanced the accuracy of the surface energy balance models developed in the course of this study.
of z0 was estimated than it would be for debris covered glaciers (Miles et al., 2017;Quincey et al., 2017). Uneven or heterogeneous ice surface such as sastrugis, crevasses, channels, and penitents could greatly affect ice surface aerodynamic surface roughness and it would be hard to estimate its z0 based on a profile method. SfM estimation of z0 might be a good 520 choice at macro-scale. In the accumulation season, more attention would be needed to be paid to spatial and temporal variations of z0 as z0 is a key parameter for sublimation calculation during this period. Studies have indicated that the Lettau (1969) approach calculated z0 dependent on plot scale and resolution. In this study, we only select 1×1 m scale at 1mm resolution to study its spatial and temporal variability. Further comparative studies of z0 are needed at different scales and resolutions.