Spectral attenuation coefficients from measurements of light transmission in bare ice on the Greenland Ice Sheet

Light transmission into bare glacial ice affects surface energy balance, biophotochemistry, and light detection and ranging (lidar) laser elevation measurements but has not previously been reported for the Greenland Ice Sheet. We present measurements of spectral transmittance at 350– 900 nm in bare glacial ice collected at a field site in the western Greenland ablation zone (67.15 N, 50.02W). Empirical irradiance attenuation coefficients at 350–750 nm are ∼ 0.9– 8.0 m−1 for ice at 12–124 cm depth. The absorption minimum is at ∼ 390–397 nm, in agreement with snow transmission measurements in Antarctica and optical mapping of deep ice at the South Pole. From 350–530 nm, our empirical attenuation coefficients are nearly 1 order of magnitude larger than theoretical values for optically pure ice. The estimated absorption coefficient at 400 nm suggests the ice volume contained a light-absorbing particle concentration equivalent to ∼ 1–2 parts per billion (ppb) of black carbon, which is similar to pre-industrial values found in remote polar snow. The equivalent mineral dust concentration is ∼ 300–600 ppb, which is similar to values for Northern Hemisphere warm periods with low aeolian activity inferred from ice cores. For a layer of quasi-granular white ice (weathering crust) extending from the surface to ∼ 10 cm depth, attenuation coefficients are 1.5 to 4 times larger than for deeper bubbly ice. Owing to higher attenuation in this layer of near-surface granular ice, optical penetration depth at 532 nm is 14 cm (20 %) lower than asymptotic attenuation lengths for optically pure bubbly ice. In addition to the traditional concept of light scattering on air bubbles, our results imply that the granular near-surface ice microstructure of weathering crust is an important control on radiative transfer in bare ice on the Greenland Ice Sheet ablation zone, and we provide new values of flux attenuation, absorption, and scattering coefficients to support model development and validation.


Introduction
Understanding the transmission, absorption, and scattering of light in ice is important for snow and ice energy balance modelling (Brandt and Warren, 1993), lidar remote sensing of snow surface elevation and grain size (Deems et al., 2013;Yang et al., 2017), primary productivity beneath sea ice Grenfell, 1979), bio-photochemical cycling in ice and snow (France et al., 2011), and theoretical predictions of "Snowball Earth" paleoclimates (Dadic et al., 2013;Warren et al., 35 2002). Each of these applications requires knowledge of the vertical distribution of light attenuation in ice, which for a medium (such as glacier ice) that both absorbs and scatters light is specified by the spectral attenuation coefficient: where abs (m -1 ) is the absorption coefficient, sca (m -1 ) is the scattering coefficient, and all are functions of wavelength, λ.
This study reports on the irradiance attenuation coefficient att of bare glacier ice in the Greenland Ice Sheet ablation zone, a 40 critical parameter needed to calculate subsurface absorption and scattering of transmitted radiation that to our knowledge has received no direct field study.
Measurements of att in snowpack and sea ice indicate three main sources of variation with relevance to geophysical applications. First, the magnitude of att is primarily controlled by ice microstructure (e.g., the size, shape, orientation, and 45 number of air bubbles, ice grains, and cracks), via its control on sca (Dadic et al., 2013;Libois et al., 2013;Light et al., 2004Light et al., , 2008. For the range of air bubble sizes (~10 -6 ) and ice grain sizes (~10 -1 -10 -3 ) observed in glacier ice, sca is effectively independent of wavelength in the visible and near-infrared spectrum (Bohren, 1983;Dadic et al., 2013;Perovich, 1996).
Spectrally, att is low in the near-ultraviolet and blue-green (~250-600 nm) where abs is extremely low (<10 -8 ), and progressively higher for wavelengths >600 nm, where abs rapidly increases up to its maximum value (~10 -2 ) at the far end of 50 the solar spectrum (Warren and Brandt, 2008). Vertically, att is at a maximum at the incident boundary (the snow or ice surface) where a portion of upwelling radiation (i.e., transmitted flux reflected upwards) escapes the ice volume before re-reflection downward. Within this near-surface optical boundary layer (Bohren and Barkstrom, 1974), attenuation rates rapidly decrease with depth to an asymptotic value as multiple scattering establishes an isotropic (diffuse) radiation field (Briegleb and Light, 2007;Warren, 1982). For fine-grained dry snow, a few cm depth is typically sufficient to reach the 55 asymptotic regime where att is constant (Brandt and Warren, 1993). For sea ice the depth required is typically larger and can exceed >20 cm depending on near-surface ice microstructure and the vertical location of the refractive boundary if present (Grenfell, 1991;Grenfell and Maykut, 1977). Attenuation coefficients are also influenced by the horizontal distribution of ice type and surface cover  but this source of variation is not examined here.

60
In addition to experimental values obtained from measurements of light transmission in ice or snow, att is obtained analytically from optical theory (Bohren, 1987;Warren et al., 2006). Light attenuation in pure ice is specified analytically by the complex refractive index: where re is the real part of the complex refractive index (~1.31 in the visible), im is the imaginary part, and: 65 abs ice (λ) = 4 im (λ) (3) is the absorption coefficient of pure ice (Warren et al., 2006;Warren and Brandt, 2008).
Light attenuation in glacier ice differs from pure ice owing to compositional and structural factors that control scattering and absorption, such as the size, geometry, and vertical distribution of embedded light absorbing particles (LAPs) and light 70 scattering air bubbles and ice grains of size larger than wavelength (Askebjer et al., 1997;Picard et al., 2016;Price and Bergström, 1997b;Warren et al., 2006). Analytical methods typically approximate ice and snowpack as homogeneous plane-parallel slabs of spheres having the same volume-to-surface-area ratio (i.e., optically-equivalent grain size) as the collection of non-spherical ice grains and air bubbles in realistic ice (Brandt and Warren, 1993;Grenfell and Warren, 1999;Wiscombe and Warren, 1980). Mie theory is used to calculate the single-scattering properties and two-stream radiative transfer 75 approximations are used to calculate multiple scattering and bulk absorption in the ice volume (Bohren and Barkstrom, 1974;Mullen and Warren, 1988;Wiscombe and Warren, 1980). The single-scattering properties can also be derived from the ratio of surface area to mass (i.e., specific surface area) with or without the assumption of spherical scattering geometry (Kokhanovsky and Zege, 2004;Malinka, 2014), as applied to the highly scattering granular surface layer on sea ice (Malinka et al., 2016). Models of the prior form have been used to calculate subsurface meltwater production caused by penetration of 80 solar radiation in ice both in Greenland (van den Broeke et al., 2008;Kuipers Munneke et al., 2009) and Antarctica (Brandt and Warren, 1993;Hoffman et al., 2014;Liston et al., 1999aListon et al., , 1999bListon and Winther, 2005). However, theoretical values for att are rarely validated experimentally, and to our knowledge no such experimental values exist for near-surface glacier ice.

85
In addition to ice surface energy balance, understanding light attenuation in ice is important for interpreting interactions between visible-wavelength light sources and ice surfaces, for example laser altimetry measurements of ice surface elevation (Deems et al., 2013;Gardner et al., 2015;Greeley et al., 2017). The reciprocal of att is the attenuation length, or the average distance travelled by a photon before attenuation by absorption or scattering (Ackermann et al., 2006). In the context of altimetry, the attenuation length is sometimes referred to as the penetration depth, or the average depth to which the 90 electromagnetic signal penetrates before it is backscattered to the atmosphere (Ridley and Partington, 1988;Rignot et al., 2001;Zebker and Weber Hoen, 2000). The laser altimeter onboard Ice, Cloud, and Land Elevation Satellite-1 (ICESat-1) transmitted 1064 nm laser pulses to measure the distance (range) between the satellite and ice sheet surfaces (Schutz et al., 2005). Photons with wavelength 1064 nm penetrate snowpack no more than a few centimetres (Brandt and Warren, 1993;Järvinen and Leppäranta, 2013). This length scale is smaller than typical laser altimetry surface elevation errors due to ice and snow surface 95 roughness and geolocation uncertainty (Deems et al., 2013). In contrast, the laser altimeter onboard ICESat-2 transmits 532 nm laser pulses (Markus et al., 2017). Ice is ~10 times more transparent at 532 nm than at 1064 nm (Warren and Brandt, 2008), and photons at 532 nm may penetrate many tens of centimetres into glacier ice. These subsurface scattered photons may introduce a range bias in ICESat-2 surface elevation retrievals over glacier ice, similar to radar penetration into snow (Brunt et al., 2016;Gardner et al., 2015;Greeley et al., 2017;Smith et al., 2018). To our knowledge no in situ observations of 532 100 nm optical penetration depth for bare glacier ice exist, precluding field validation of penetration depth obtained from theoretical radiative transfer models.
The purpose of this investigation is to provide experimental values for att obtained from measurements of solar flux attenuation in bare ice in the Greenland Ice Sheet ablation zone, and to compare them with theoretical values for att obtained 105 from the two-stream analytical solution (c.f. Eq. 26 Bohren, 1987;Schuster, 1905). We benchmark our field estimates against the two-stream solution because of its wide use in surface energy balance models applied to snow and ice. In Sect. 2 we describe the field measurements and the optical theory used to interpret the solar flux attenuation. In Sect. 3 we report values for att obtained from our measurements, compare them with values obtained from two-stream theory, and propose a simple empirical model that accounts for enhanced near-surface attenuation. In Sect. 4 we discuss how our att values differ from 110 prior experimental values acquired in sea ice, snowpack, and deep South Pole glacial ice, and the implication of these differences for modelling radiative transfer in bare glacier ice. To demonstrate the broader implications of our study, we suggest how our findings can be used to improve models for subsurface heating of ablating glacier ice.

Transmittance measurements 115
Ice transmittance was measured on 20 July and 21 July 2018 in the Kangerlussuaq sector of the western Greenland Ice Sheet.
The study site is located ~1 km from the ice sheet margin at 840 m a.s.l. (67.15 o N, 50.02 o W). Subsurface (in-ice) spectral irradiance was measured at ~0.35 nm spectral resolution in the wavelength range 350-900 nm with an Ocean Optics® JAZ spectrometer. Light was guided from the ice interior to the spectrometer with a 3 mm diameter Kevlar-sheathed fibre optic cable fitted inside a 2 m long insulated white PVC tube (Fig. 1). The fibre was attached at one end to an irradiance sensor 120 consisting of a 90 o collimating lens adapter and a remote cosine receptor (RCR) with a Spectralon TM diffusing element. The RCR lens barrel was wrapped in white PTFE tape and set 2 mm out from the PVC tube exterior to act as a contact horizon between its diffusing element and the ice. The system was operated from a battery-powered computer running the Ocean Optics® OceanView software placed on a tripod platform oriented 180 o away from the sun and 2.5 m horizontal distance from the measurement location. 125 To access the interior of the ice, holes were drilled horizontally into a 2-m high sidewall of a natural ice feature with a battery powered hand drill fitted with a 3 cm diameter Kovacs auger bit. To drill these holes, the auger was advanced into the sidewall approximately 20 cm, levelled horizontally with a digital spirit level, and the sequence was repeated to 2 m horizontal depth.
The PVC tube-fibre optic assembly was then inserted into the hole, RCR facing upward, and a 2 m long ruler was shimmed 130 under the bottom of the PVC tube to ensure the RCR barrel preserved contact with the overlying ice thus minimizing stray light contamination into the RCR field of view. Ice shavings were packed around the drill hole to prevent light reflection into the hole. Spectral irradiance was measured using a 20-scan average with 0.0228 s integration time per scan, yielding 0.46 s total integration time per irradiance measurement. Irradiance measurements were recorded at 1 Hz frequency for thirty seconds yielding thirty irradiance profiles at each depth, after which the tube was removed, the next hole was drilled, and the sequence 135 was repeated, working from the top hole toward the bottom on 20 July, and from the bottom hole toward the top on 21 July.
Background upwelling and downwelling spectral irradiances were measured continuously at 2 m height above the ice surface 140 ~3 m away from the in-ice measurements with a dual-channel Ocean Optics JAZ spectrometer. These data were recorded at 1 min frequency using a 30-scan average with 0.011 s integration time. Light was guided to the spectrometer via two 3 m fibre optic cables attached to two RCRs mounted in upward-looking and downward-looking orientation on a 2 m long horizontally levelled boom attached to a vertical mast frozen into the ice. The horizontal boom became unstable on 21 July and the upward-looking RCR was moved to the vertical mast; the downward-looking RCR was decommissioned. 145 The surface-based spectrometer was calibrated for absolute irradiance in a controlled setting prior to the field experiment using an Ocean Optics HL-3P radiometrically calibrated halogen light source. During the field experiment, the in-ice spectrometer was cross-calibrated to the surface spectrometer by holding it level above the ice surface in an upward-looking orientation ~3 m away from the surface spectrometer. Cross-calibration irradiance profiles were collected on 20 July and 21 July immediately 150 prior to subsequent in-ice measurements. All in-ice irradiances are cross-calibrated to the surface spectrometer as a pre-processing step prior to further analysis.
Dark current spectra were recorded prior to each irradiance measurement as input to the OceanView automated dark current correction module. To measure dark current, the RCR lens barrel was capped with a custom-fit opaque metal cap provided by 155 Ocean Optics. OceanView adjusts these spectra in real-time for changes in integration time and for charge leakage if detected, corrects the nonlinear analogue-to-digital response of the linear silicon charge coupled device, and applies a boxcar smoothing over adjacent pixels to further reduce noise. Following these automated corrections, the opaque cap was left in place and residual dark current (noise) was recorded with the reference spectrometer in its identical setup during the experiment and with the in-ice spectrometer held level above the ice surface in an upward-looking orientation. These residual dark current spectra 160 are treated as systematic errors and are subtracted from all irradiance profiles as a pre-processing step prior to analysis (Fig.   2a).

Weather conditions
The 20 July experiment was conducted under low, thick cloud cover with light rain and no direct sun, ideal conditions for estimating the attenuation of diffuse light in ice. The 21 July experiment was conducted under higher, thinner cloud cover with 165 no rain and very brief periods of intermittent direct sun (see Fig. 1). The effect of intermittent direct sun was easily identified in the in-ice irradiance measurements as a rapid increase in light intensity, which only occurred during the third measurement on 21 July. This was mitigated by averaging over the first ten in-ice irradiance pro files for that measurement, prior to the rapid increase in light intensity, and discarding the remainder.

Ice thickness and density 170
The ice thickness between detector positions was measured to the nearest millimetre with a metre stick and converted to units of solid ice thickness with the relation: where Δh is in-situ ice thickness between detector positions, is in-situ ice density, ice is solid ice density (917 kg m -3 ), and Δz is solid ice thickness between detector positions. Two separate observers made ten independent measurements of Δh. In 175 addition, one observer made 41 replicate measurements of an ablation stake using the same metre stick, yielding a mean difference and standard error on Δh. The ice density was measured on a 1.2 m ice core extracted at the measurement location with a Kovacs Mark IV corer (www.kovacsicedrillingequipment.com) (Fig. 3). The ice core was split along natural breaks into three segments that were measured to the nearest millimetre with a calliper and weighed to the nearest gram on an Acculab digital scale. 180

Experimental asymptotic flux attenuation coefficients and ice surface albedo
Spectral asymptotic flux attenuation coefficients are estimated by fitting a Bouguer-law exponential decay model to the in-ice irradiance profiles (Grenfell and Maykut, 1977): where att (λ) is the asymptotic flux attenuation coefficient, z (λ) is in-ice spectral irradiance at depth z, 0 (λ) is background 185 downwelling spectral irradiance, z 0 is the ice surface, and: is spectral transmittance. The optical depth z (λ) is a dimensionless path length that scales the physical thickness of a layer by its attenuation rate: Estimates of att (λ) for each spectral band are obtained by solving a linear equation of the form: where 0 is a parameter (y-intercept), Δz = z − z 0 is ice thickness, Δz is an error term that represents ice thickness measurement uncertainty, and Δ is an error term that represents optical path measurement uncertainty. Equation 8 is solved by Maximum Likelihood Estimation (MLE), which gives an unbiased estimate of the slope when measurement errors are 195 present in both the independent and dependent variables (see Sect. 2.9) (York et al., 2004).
The attenuation length att (λ) is the inverse of att (λ) and is analogous to the photon mean free path or transport length (Ackermann et al., 2006). It is equivalent to the path length in ice required to attenuate irradiance to 37% (1/ ) of its incident intensity, i.e., the path length at which = 1/ and = 1: 200 . (9) The ice surface spectral albedo is the ratio of the upwelling spectral irradiance to the downwelling spectral irradiance: and the broadband albedo is: (11) 205

Asymptotic flux attenuation coefficients
Theoretical att (λ) values are calculated using the asymptotic solution to the delta-Eddington two-stream radiative transfer approximation (Joseph et al., 1976;Schuster, 1905): where ext (λ) is the extinction efficiency, eff is the effective scattering particle radius (m), (λ) is the average cosine of the 210 scattering angle, also referred to as the asymmetry parameter, and (λ) is the single-scattering albedo: where att (λ) and sca (λ) are the single-scattering attenuation coefficient (m -1 ) and scattering coefficient (m -1 ), respectively.
Equation (12) describes light attenuation by multiple scattering and absorption in a homogeneous plane-parallel slab of absorbing spheres far from any boundaries (Bohren, 1987). 215 To estimate eff , Eq. (12) is inverted and solved by iteration for the value of eff that minimizes the difference between measured and calculated att at λ = 600 nm. This method assumes that all absorption at 600 nm is due to ice (Warren et al., 2006). If absorption was influenced by LAPs eff would be over-estimated. Values for ext (λ), (λ), and (λ) are obtained from Mie scattering algorithms (Mätzler, 2002) provided as MATLAB® code with input (λ) from Warren and Brandt 220 (2008). The Mie solutions at each wavelength are integrated over a Gaussian size distribution (N=1000) of scattering radii ( r = eff , r = 0.15 eff ) to eliminate ripples associated with Bessel function solutions to the Mie equations (Gardner and Sharp, 2010). The optimal eff value is 3.1 mm, which corresponds to a specific surface area of 1.05 m 2 kg -1 . These values are used in all subsequent calculations. Warren et al. (2006) developed a method to estimate abs for pure ice (i.e. abs ice ) from measurements of flux attenuation in snow in Antarctica. The method relies on three assumptions: 1) the value of abs ice at the reference wavelength (λ 0 = 600 nm)

Flux absorption coefficients 225
is known accurately, 2) the value of att at λ 0 is not affected by LAPs in the measured snow or ice, and 3) (λ) varies so little as to be effectively independent of wavelength in the spectral range considered (here the near-UV and visible). Warren et al. (2006) verified the validity of these assumptions for the spectral range 350-600 nm and obtained the following relation (Eq. 230 15 of that paper) between flux attenuation and flux absorption: Equation 14 assumes that abs is not affected by LAPs at the reference wavelength (600 nm) but the relation can be used to estimate abs at any other wavelength, including those where absorption is affected by LAPs. At those wavelengths, Eq. (14) will predict values for abs higher than abs ice if LAPs are present in the measured snow or ice volume, due to the influence of 235 LAPs on att .
The inferred abs values can be related to a mass absorption cross-section (MAC) (Doherty et al., 2010): where is the spectral MAC (m 2 kg -1 ) and is the mass-mixing ratio of LAPs in the ice volume (g LAPs g -1 ice). We exploit 240 this to interpret differences between our theoretical and experimental values of att on the basis of differences between abs ice (Warren et al., 2006) and the abs values that we obtain for glacier ice from Eq. (14). To provide context, we use representative values of for black carbon BC and insoluble mineral dust (hereafter 'dust') dust to estimate corresponding equivalent mass mixing ratios eq BC and eq dust (Di Mauro et al., 2017;Doherty et al., 2010). The "equivalent" mass mixing ratio is the mass-mixing ratio of each LAP species required to explain the difference between abs ice and our inferred abs values at a 245 reference wavelength, and follows a similar approach used to infer LAP absorption in snowpack (Tuzet et al., 2019). For BC , we use 6 m 2 g -1 as a representative MAC at 550 nm and an absorption Ångstrom exponent range 0.8-1.9 to scale this value to 400 nm (Doherty et al., 2010). For dust , we use 0.013 m 2 g -1 at 550 nm (Di Mauro et al., 2017) and an absorption Ångstrom exponent range 2-3 (Doherty et al., 2010). We note that these descriptive estimates provide context for discussion; actual LAP species concentrations were not measured. 250

Near surface effects
Equations 7 and 12 are applicable at distances far enough from the incident boundary (here the ice surface) that the radiation field is diffuse and att is constant with depth. Near the ice surface the radiation field is converted via multiple scattering from direct to diffuse flux, and attenuation may be enhanced by direct reflection, enhanced scattering and/or absorption by the granular near-surface ice microstructure, or specular reflection at the ice surface, depending on its roughness (Dadic et al., 255 2013;Light et al., 2008;Mullen and Warren, 1988) . To account for non-diffuse near-surface attenuation, we define a piecewise optical depth: where ′ is an effective attenuation coefficient for the near-surface non-diffuse layer and z ′ is a depth chosen to partition this layer from the interior diffuse region. We estimate ′ with a centred finite-difference form of Eq. (7): 260 Here, Δz ′ = 12 cm and z ′ is the 12 cm in-ice irradiance measured on 20 July. Accordingly, the asymptotic attenuation length (Eq. 9) is distinguished from an effective penetration depth λ to include the effect of near-surface attenuation. The attenuation length is the depth at which = 1. Setting (λ) = 1 in Eq. (16) and solving for z yields: Equation 16 gives estimates of spectral transmittance that account for non-diffuse near-surface attenuation but relies on knowledge of ′ , which is sensitive to the spectral composition and directional distribution of 0 and the structure and composition of the near-surface ice (Grenfell and Maykut, 1977;Light et al., 2008). To generalize the magnitude of near-surface attenuation, we calculate the fraction of downwelling spectral irradiance that transmits the non-diffuse layer weighted by the downwelling spectral irradiance: 270 The 0 parameter is analogous to the 0 parameter introduced by Grenfell and Maykut (1977) to partition the fraction of solar irradiance absorbed in the upper 10 cm of sea ice, which they termed the "Surface Scattering Layer" (SSL), and the ice interior, in which radiation is exponentially absorbed at a constant rate: The 0 parameter has been widely adopted in energy balance models of glaciers and sea ice to compute subsurface flux divergence (heating rates) when radiation penetration is considered important (Bintanja and Van Den Broeke, 1995;Hoffman et al., 2014;Holland et al., 2012). For example, the sea ice component of the Community Earth System Model (CESM) uses 0 = 70% for the visible (200-700 nm) and 0 = 0% for the infrared (700-5000 nm) (Briegleb and Light, 2007). The important distinction is that 0 partitions the absorbed flux whereas 0 partitions the downward flux (Brandt and Warren, 280 1993). For both 0 and 0 , we set Δz ′ to 10 cm for consistency with prior work (Grenfell and Maykut, 1977;Light et al., 2008;Maykut and Untersteiner, 1971).

Monte Carlo simulations of detector interference
We developed a Monte Carlo radiative transfer model to estimate the effect of detector interference on measured irradiances and fitted att values, following methods developed to simulate light propagation in biological tissue, ocean waters, and sea 285 ice (Leathers et al., 2004;Light et al., 2003;Wang et al., 1995). Photon scattering is specified by a Henyey-Greenstein scattering phase function with single-scattering properties ext (λ), (λ), and (λ) inferred from our optical measurements (Sect. 2.5). A complete technical description is given in Appendix 2, where model accuracy is verified by comparison with benchmark solutions to the radiative transfer equation (van de Hulst, 1980).

290
In the Monte Carlo simulations, photons are launched from an irradiance sensor on a detector rod with dimensions identical to those reported in this study. In the ideal (baseline) simulation, photons originate from an isotropic point source and propagate through ice until they transmit the surface or are terminated by absorption. Detector interference is investigated by repeating the Monte Carlo with an ideal cosine source function describing the angular response to radiance of the RCR, and with a non-ideal (empirical) angular response function (Fig. 2), with and without scattering and absorption interference by the PVC 295 detector rod. The detector rod albedo rod ≈ 0.4 is calculated from the absorption spectra of polyvinyl chloride (Zhang et al., 2020); scattering by the rod is assumed isotropic. The Monte Carlo is integrated over 10 000 interactions at nine wavelengths in 50 nm increments from 350 nm to 750 nm, allowing us to fit the wavelength dependence of the estimated systematic uncertainty in simulated att values.

Uncertainty propagation 300
Unless stated otherwise, all statistical uncertainties reported in this paper are standard errors that correspond to 68% confidence intervals around the mean (Taylor and Kuyatt, 1994). For an individual measurement with standard deviation and sample size ≥ 30 the standard error is /√ . For < 30, standard errors are scaled by a critical t-value drawn from the Student's t-distribution. Standard errors for combined quantities are propagated in quadrature and hereafter referred to as combined uncertainty. The combined uncertainty for spectral irradiance (λ) is: 305 where * is the standard deviation of the high-frequency irradiance spectra before dark-noise correction and is the standard deviation of the high-frequency dark-noise spectra. An analogous procedure is used to estimate the combined uncertainty for calibrated irradiance. The combined uncertainty for spectral transmittance (λ) is: where z cal and 0 are the combined uncertainties for calibrated in-ice irradiance and dark-noise corrected surface downwelling irradiance, respectively. The combined uncertainty for optical depth is: and the combined uncertainty for att is: Equation 24 gives a first order description of σ due to statistical propagation of measurement uncertainties, neglecting higher-order interaction terms. A description of the statistical uncertainty in fitted att (λ) values is given by the MLE estimate of the regression slope of Eq. 8, which can be expressed in terms of an error model as: where ̂ and Δẑ are the true but unobserved (due to measurement error) optical depth and ice thickness and Δz~( 0, Δz 2 ) 320 and Δ~( 0, Δ 2 ) are normally distributed error terms. Unlike Ordinary Least Squares, MLE gives an unbiased estimate of the slope and standard error of a linear functional relationship between two variables measured with error (York et al., 2004).
The method has been used in similar studies to infer optical coefficients (Zieger et al., 2011). The MLE standard errors for att (λ) are adjusted for − 2 degrees of freedom with a two-sided t-statistic (Cantrell, 2008) and combined in quadrature with systematic uncertainty estimated from Monte Carlo simulation to estimate total combined uncertainty for reported att (λ) 325 values.

Spectral transmittance
Four in-ice irradiance spectra were collected at 12 cm, 36 cm, 58 cm, and 77 cm depth below the ice surface on 20 July (hereafter referred to as Layer A) (Fig. 4a), and at 53 cm, 67 cm, 82 cm, and 124 cm on 21 July (hereafter referred to as Layer 330 B). At all depths, spectral transmittance is maximum at 430 nm and maintains relatively stable and high values up to ~500 nm in the visible, beyond which decreases into the red end of the visible spectrum following the well-known exponential increase in ice absorptivity (Fig. 4c). Maximum values vary from 78% at 12 cm to 45% at 77 cm. For wavelengths >500 nm, rapidly decreases both with wavelength and with depth; beyond ~800 nm nearly all incident light is attenuated within 36 cm of the ice surface, although substantial attenuation is apparent in the 12-36 cm depth region (Fig. 4b). The standard 335 deviation of the 1 Hz raw data is <1 W m -2 nm -1 at all wavelengths, consistent with field observations of thick cloud cover and diffuse light conditions described in Sect. 2.6. Instrumental noise and high-frequency measurement variations propagate as ± 1.6% uncertainty on for wavelengths between 400-600 nm, ± 1-8% for wavelengths between 350-400 nm, where instrumental noise is higher, and ±1-12% uncertainty for wavelengths between 600-750 nm, where noise is higher and light levels are low. 340 A, the minimum in att is at 390 nm, blue-shifted relative to the maximum in at 430 nm. For Layer B, the minimum is at 397 nm. The coefficient of determination (r 2 ) ranged from 0.96-1.0 (p<0.01), with a median value of 0.98, suggesting the data are described appropriately by the Bouguer-law exponential decay model up to ~700 nm, beyond which measured values of in-ice irradiance at 58 cm and 77 cm depth were too low to reliably fit att values (see Fig. 4b and Fig. 5c). For Layer B values, low light levels prevented fits beyond ~650 nm. 350

Experimental flux attenuation coefficients and albedo
Albedo spectra correspond closely to patterns in transmittance and attenuation (Fig. 5c). The near-UV and blue wavelengths that efficiently transmit ice mostly re-emerge as reflected light, owing to the extremely low values of ice absorptivity in the wavelength range 350-500 nm where albedo is maximum (Gardner and Sharp, 2010;He and Flanner, 2020;Warren et al., 2006). The maximum measured albedo value is 0.81±0.004 at 452 nm, further red shifted from the minimum in att and the 355 maximum in . All three quantities have low variability near the minimum; albedo is 0.79 at 390 nm. The broadband albedo (Eq. 11) for the 350-900 nm wavelength range is 0.70±0.006, which is high but not atypical for melting white ice under overcast skies (Bøggild et al., 2010).

Theoretical flux attenuation coefficients
Asymptotic att values predicted by two-stream theory for optically-clean bubbly ice are nearly one order of magnitude lower 360 than field estimates for wavelengths <500 nm, where very small concentrations of LAPs in the measured ice volume dominate absorption (compare dotted grey line to solid blue line, Fig. 6) (Warren et al., 2006). In contrast, field-estimates and two-stream theory converge at wavelengths >540 nm where absorption is dominated by grain-size effects (He et al., 2017;Libois et al., 2013). The magnitude of inferred absorption enhancement in the visible due to LAPs (the quantity i in Eq. 15) varies from 0.009-0.015 m -1 at 350-530 nm. The equivalent black carbon concentration eq BC inferred at 400 nm is 1-2 ng g -1 for both 365 Layer A and Layer B, where the range covers uncertainty in both the absorption spectra and the absorption Ångstrom exponent (Doherty et al., 2010). The equivalent mineral dust concentration eq dust is ~344-620 ng g -1 for Layer A and 303-545 ng g -1 for Layer B. Monte Carlo simulations without detector interference replicate both asymptotic theory for clean bubbly ice (i.e., when forced with abs ice ) and field estimates when forced with abs values inferred from our optical measurements (solid blues squares, Fig. 6). Monte Carlo simulations of detector interference are discussed further in Sect. 4. 370

Near-surface attenuation and effective penetration depth
Near the ice surface irradiance is not attenuated exponentially and Bouguer's law does not hold, as indicated by the y-intercepts of the straight lines in Fig. 5b at values <100%. Effective ′ values (Eq. 17) for the quasi-granular 0-12 cm layer are ~1.5 times higher than att values for interior bubbly ice at 12-77 cm depth for wavelengths >570 nm and are up to 4 times higher between 400-570 nm (Fig. 7). Owing to higher near-surface attenuation, transmitted irradiance z is 375 overestimated by 10-60% if Bouguer's law is applied to the incident downwelling irradiance 0 using asymptotic att values, with median over-estimation 23% (Fig. 8a). In contrast, the piecewise optical depth (Eq. 16) predicts z to within 12% of measured values for all wavelengths between 350-700 nm with median error 3%. Integrated over these wavelengths, χ 0 is 0.68 and 0 is 0.66, suggesting 66% of the total incoming irradiance within this spectral region is absorbed at depths below 10 cm. If att is used rather than ′ to calculate χ 0 and 0 , the respective values are 0.81 and 0.79. 380 Stated in terms of penetration depth, eff varies from 12-84 cm between 350-700 nm. These values are 13-44% lower than attenuation lengths att inferred from empirical asymptotic att values. Specifically at 532 nm, eff is 52 cm, or 10 cm lower than the 62 cm empirical att value, and 14 cm lower than the 66 cm theoretical att value for optically pure bubbly ice. These results point to the potential for reduced optical penetration due to enhanced scattering and absorption on or near the ice surface, 385 as well as within the ice volume where small LAP concentrations reduce optical backscattering due to enhanced absorption.
For smooth ice surfaces, attenuation may be enhanced by refraction at the ice-air interface (Mullen and Warren, 1988). If present, a refractive boundary would enhance near-surface attenuation via external specular reflection, and possibly via enhanced near-surface absorption of the internally reflected downward flux. Following Briegleb and Light (2007), Eq. 20-24, 390 we calculate the external diffuse specular reflectivity for a flat ice surface to be 0.063, meaning specular reflection could enhance attenuation by up to 6.3%. This value is smaller than the 18-44% near-surface attenuation implied by the y-intercepts in Fig. 5b, suggesting specular reflection alone cannot explain the discrepancy. Instead, we suggest that enhanced scattering by the granular near-surface ice microstructure, together with absorptive impurities, enhance near-surface light attenuation at our field site, consistent with observations of the granular and porous surface layer on sea ice (Grenfell and Maykut, 1977;395 Light et al., 2008).

Uncertainty analysis
The effect of random and systematic uncertainties on our optical measurements and fitted att values is evaluated with Monte Carlo simulation and statistical analysis. We considered systematic uncertainties in detector positions, spectrometer sensitivity to dark current, the non-ideal angular response of the irradiance sensor, and attenuation interference by the PVC detector rod. 400 The detector positions are known to within 0.9±0.4 cm from independent measurements of the vertical ice thickness Δh. The in-situ ice density varied from 801-888 kg m -3 between 4-124 cm where irradiances were measured. The variation in was examined by repeating the analysis with Δz values computed with a single depth-weighted average applied to each Δh, and with values estimated for each Δh from linear and cubic interpolation of the vertical density profile. The maximum Δz 405 difference was 0.9 cm. The att values differed by <1%, and r 2 values were nearly identical. We use the depth-weighted average values to calculate Δz, which are 835 kg m -3 and 855 kg m -3 for the measurements collected on 20 July and 21 July, respectively.
Detector position uncertainty was further assessed by fitting att values with an ensemble of 10 000 Δz values perturbed with 410 random errors drawn from a normal distribution ( = 0.9 cm, = 0.4 cm). At all wavelengths, the chance of obtaining a fitted att value >2% from the mean value was <5%. We take 2% as a conservative estimate of systematic uncertainty due to ice thickness measurement bias.
As described in Sect. 3, all irradiance spectra are corrected for residual dark noise. The noise may have varied during the 415 experiment, and dark noise measurements with the in-ice spectrometer were made on the surface, rather than within the ice.
To assess possible bias, we fit att values with and without residual dark noise correction. The mean difference was -0.01±0.13% averaged over the 350-700 nm wavelength range. For a few discrete wavelengths between 350-400 nm and 700-750 nm, differences approached 2%. These wavelengths are those with the highest dark noise in the reference spectrometer (Fig. 2). At wavelengths between 400-700 nm, differences were <0.5%. 420 Monte Carlo simulations indicate a possible +1-14% systematic bias due to detector interference for Layer A values, and +2-8% for Layer B values ( Fig. 9; also see purple stars minus solid squares, Fig. 6). The high end of this range applies to the wavelength region of minimum absorption ~350-450 nm. The simulated bias is within statistical uncertainty at wavelengths >450 nm for Layer A and at wavelengths >400 nm for Layer B (Fig. 9). The non-ideal cosine response of the RCR and the 425 presence of the detector rod both tend to increase att values relative to the ideal case, as expected given the low albedo of the detector rod. However, detector interference is masked somewhat by the presence of LAPs, as indicated by the larger simulated interference for bubbly ice without LAPs (see dotted grey line and associated Monte Carlo values, Fig. 6). Overall, the combined statistical and systematic uncertainty for the 350-450 nm region is <20% for Layer A values and <14% for Layer B values, and as low as ~5% for wavelengths >450 nm (Fig. 9). 430

Comparison with attenuation spectra for sea ice, snowpack, and deep glacial ice
We report spectral measurements of near-UV and visible light transmission in bare ablating glacier ice. These measurements are used to estimate irradiance attenuation coefficients att for the spectral range 350-750 nm. Prior studies quantified att for sea ice and snowpack (c.f. Fisher et al., 2005;Frey et al., 2011;Gerland et al., 2000;Grenfell and Maykut, 1977;Järvinen and 435 Leppäranta, 2013;King and Simpson, 2001;Light et al., 2008;Meirold-Mautner and Lehning, 2004;Pegau and Zaneveld, 2000;Picard et al., 2016;Tuzet et al., 2019;Warren et al., 2006). Scattering and absorption coefficients were quantified for compressed South Pole glacial ice at 800-2350 m depth by the AMANDA (Antarctic Muon and Neutrino Detector Array) experiment (Ackermann et al., 2006;Askebjer et al., 1995Askebjer et al., , 1997. For South Pole ice at 800-1000 m depth, visible and near-UV light scatters on air bubbles, below which bubbles transition under pressure to clathrates and light scatters on dust grains (Price 440 and Bergström, 1997b). In the bubbly ice regime studied by AMANDA, sca values at 532 nm are ~1-3 m -1 , comparable to the 1.6 m -1 value quantified in this study. Light scattering in the dusty-ice regime (>1000 m depth) is not comparable to this study; absorption by dust is discussed below. Fig. 10 compares our att spectra for glacier ice to seven previously published spectra for snowpack and sea ice. In general, 445 glacier ice is the most transparent structure examined with the exception of multi-year and first-year interior sea ice at wavelengths >540 nm . Light attenuation in sea ice is controlled by its unique vertical composition including brine inclusions, air pockets, solid salts, sea ice algae, dissolved organic matter, water saturation, and radiative interactions between the ice and underlying ocean (Perovich, 1996). The latter factor, together with differences in optically-equivalent grain size, may explain the low attenuation at longer wavelengths for sea ice shown here. Relative to 450 snowpack in Greenland and Antarctica (Järvinen and Leppäranta, 2013;Meirold-Mautner and Lehning, 2004;Warren et al., 2006), attenuation by glacial ice has similar spectral structure but is lower at all wavelengths, reflecting the higher specific surface area of fine-grained polar snow. Attenuation within the surface scattering layer (SSL) of sea ice is intermediate, with spectral structure similar to snowpack and glacial ice. Attenuation at 5 cm depth in snow near Summit, Greenland is highest of all, possibly due to direct light scattering in the near-surface optical boundary layer. The comparison demonstrates that att 455 values vary by nearly two orders of magnitude at visible wavelengths due to differences in ice structure and composition.
At visible wavelengths between 350-530 nm our field estimates of att are up to one order of magnitude larger than those obtained from two-stream theory for optically pure bubbly ice, consistent with selective absorption by mineral dust, black carbon, and microorganisms found on glaciers and ice sheet surfaces (Bøggild et al., 2010;Ryan et al., 2018;Stibal et al., 460 2017;Takeuchi, 2002;Yallop et al., 2012). For context, the absorptivity we document at 400 nm for Layer B can be explained by 1.2-1.8 ng g -1 (ppb) equivalent black carbon concentration. Values for Layer A are 1.4-2.0 ppb. Both estimates are relative to pure ice absorptivity values reported by Warren et al. (2006). These values are within the range 2±2 ppb reported for clean snow near Dye-2 on the interior Greenland Ice Sheet considered representative of pre-industrial fallout rates (Doherty et al., 2010). The equivalent mineral dust concentration is ~344-620 ppb for Layer A and 303-545 ppb for Layer B. 465 Relative to South Pole ice, our absorptivity values broadly agree with AMANDA values within two depth regions corresponding to peaks in atmospheric dust concentration during the Last Glacial Maximum (LGM) and Marine Isotope Stage 4 (MIS-4) glacial periods ~23 000 and ~66 000 years before present (Fig. 11). For these periods in Earth's history, Southern Hemisphere dust concentrations inferred from the Vostok and Dome-C ice cores are ~300-1500 ppb (Muhs, 2013;Petit et al., 470 1999). Hemispherical dust fluxes are generally synchronous at these timescales; similar peaks at LGM and MIS-4 are observed in Greenland ice cores (Ruth et al., 2003). However, Northern Hemisphere dust concentrations are several times higher (Muhs, 2013;Ruth et al., 2003), meaning correlation with South Pole absorptivity does not map age at our site. Rather, our optical measurements are consistent with the relatively low dust concentrations during Northern Hemisphere warm periods. For the western Greenland ablation zone, alternating bands of visibly dark and bright outcropping ice are associated with periods of 475 higher and lower aeolian activity during both the early Holocene (post-LGM) and late Pleistocene, with a characteristic band of older brighter interglacial ice ~0.7-1 km from the margin where our field site is located (Bøggild et al., 2010;Petrenko et al., 2006;Reeh et al., 2002;Wientjes et al., 2012). Taken together, this suggests the optical properties documented here are representative of Pleistocene interglacial ice with relatively low volumetric LAP concentration and smaller crystal diameters than Holocene ice associated with the 'dark zone' further inland (Gow et al., 1997;Petrenko et al., 2006;Wientjes et al., 2011). 480 Regarding pure ice absorptivity, our abs values support the lower bound pure ice estimate from Warren et al. (2006) (Fig.   12). The steeply sloping high values in the near-UV in the laboratory measurements (Grenfell and Perovich, 1981;Perovich and Govoni, 1991) are now understood as signatures of Rayleigh scattering on nanoscale defects in the laboratory-grown ice (Price and Bergström, 1997a). The South Pole values at 1755 m depth and 830 m depth are contaminated by trace dust 485 deposited during the late Pleistocene and early Holocene, respectively (Ackermann et al., 2006) (Fig. 11). The lowest values reported by Warren et al. (2006) (hereafter W06) were obtained by applying Eq. 7 to measurements of transmitted radiance in a single snow layer at ~90-135 cm depth near Dome C in Antarctica contaminated by ~0. Our absorption minimum is at 390 nm for Layer A values and 397 nm for Layer B values, in agreement with W06 and AMANDA. The wavelength shift in the P16 absorption minimum (430 nm) is apparent in all attenuation coefficient spectra shown in Fig. 10 that use the surface as a reference horizon, but is absent in those that use an interior layer as a reference 500 horizon. P16 used the latter method and excluded radiance measurements within 8 cm of the surface based on Monte Carlo simulations of detector interference and visual inspection of homogeneous attenuation rates, but their shifted minimum may indicate that radiance profiles in the near-UV and blue were impacted disproportionately by detector rod interference. The same can be said of our values, and may explain in part the spectral structure in our near-surface (0-12 cm) effective attenuation coefficient profile between 350-400 nm (dotted line Fig. 10). Similar spectral structure is apparent in diffuse attenuation 505 coefficients obtained in snowpack in the French Alps using the same method as P16 (c.f. Fig. 3b, Tuzet et al., 2019).
Differences aside, our inferred absorption spectrum provides new insight into the magnitude of this fundamental but uncertain optical property, and supports a conclusion that the minimum is likely <10 -2 m -1 and possibly lower (Ackermann et al., 2006;Picard et al., 2016;Warren et al., 2006).

510
In addition to the traditional concept of surface melt, visible light transmission provides an energy source for subsurface heating and internal melting of near-surface glacier ice (Cooper et al., 2018;Hoffman et al., 2014;Liston and Winther, 2005;Schuster, 2001). Prior estimates of subsurface meltwater production in bare ice used two-stream theory forced with values of abs for pure ice to calculate att and the absorbed solar flux as a function of depth below the ice surface in both Greenland and Antarctica (van den Broeke et al., 2008;Hoffman et al., 2014;Kuipers Munneke et al., 2009;Liston and Winther, 2005). The 515 influence of LAPs on subsurface meltwater production has not been quantified to our knowledge and is beyond our scope, but our results suggest LAPs enhance subsurface energy absorption in ablating glacier ice, consistent with enhanced surface melt rates caused by LAPs distributed on bare ice surfaces and within snowpack (Bøggild et al., 1996;Goelles et al., 2015;Goelles and Bøggild, 2017;Tuzet et al., 2019). From a practical perspective, this suggests that abs values for contaminated ice given here and snowpack given elsewhere (Picard et al., 2016) could provide realistic input for radiative transfer models absent 520 explicit knowledge of realistic LAP concentrations. In contrast, simulations that use the canonical pure ice absorptivity values compiled in Warren and Brandt (2008) will likely underestimate light attenuation and misrepresent the distribution of subsurface absorbed flux unless LAP concentrations are otherwise accounted for.

Conclusion
We report in-situ spectral measurements of near-UV and visible light transmission in near-surface bare glacial ice, collected 525 at a site in the western Greenland ablation zone during July 2018. In general, our empirical irradiance attenuation coefficients are nearly one order of magnitude larger in the range 350-530 nm than predicted by asymptotic two-stream theory using canonical values for the absorption coefficient of pure ice (Warren and Brandt, 2008). The absorption minimum is 0.013-0.014 ± 0.003 m -1 at 390-397 nm implying absorption length scales of 69-77 m. The volumetric scattering coefficient is 1.6 ± 0.2 m -1 at 532 nm, with an asymptotic attenuation length scale 0.62 ± 0.08 m. In addition to light scattering on air bubbles, 530 we find that light attenuation is enhanced by a layer of quasi-granular white ice that extends from the surface to ~10 cm depth at our field site. The effective penetration depth, which accounts for reduced optical transmission through this granular layer relative to interior bubbly ice, is 0.52 ± 0.07 m at 532 nm. Our co-located measurements of transmittance and albedo suggest that about 34% of cloudy sky downwelling solar irradiance at 350-700 nm is absorbed within 10 cm of the ice surface at this time and location, consistent with observations of the semi-granular surface layer on sea ice. The estimated absorption spectrum 535 suggests equivalent black carbon and mineral dust concentrations consistent with pre-industrial periods in Earth's history with low Northern Hemisphere aeolian activity, and therefore may provide a reasonable lower bound on volumetric absorption enhancement due to impurities embedded in outcropping glacial ice in the western Greenland ablation zone.

Data availability
The data is hosted by the Pangaea open access data repository (https://doi.pangaea.de/10.1594/PANGAEA.913508) and is 540 provided as a supplement to this manuscript in Appendix 2.     combined with systematic uncertainty due to detector interference estimated with Monte Carlo (error bars; ± ). The same comparison for the 21 July experiment (inset) suggests detector interference is lower and within statistical uncertainty at wavelengths >400 nm.    . Differences in attenuation magnitude at each wavelength are mostly controlled by structural differences that control scattering, whereas spectral differences are mostly controlled by differences in type and concentration of absorbing impurities.