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.02^∘ W). Empirical irradiance attenuation coefficients at 350–750 nm are ∼ 0.9–8.0 m⁻¹ 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.

1 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 (Frey et al., 2011; Grenfell, 1979), biophotochemistry in ice and snowpack (France et al., 2011), and theoretical predictions of “Snowball Earth” palaeoclimates (Dadic et al., 2013; Warren et al., 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:

$$\begin{array}{}\text{(1)}& k_{\mathrm{att}}\left(\mathit{\lambda }\right)=k_{\mathrm{abs}}\left(\mathit{\lambda }\right)+k_{\mathrm{sca}}\left(\mathit{\lambda }\right),\end{array}$$

where k_abs (m⁻¹) is the absorption coefficient, k_sca (m⁻¹) is the scattering coefficient, and all are functions of wavelength λ. This study reports on the irradiance attenuation coefficient k_att of bare glacier ice in the Greenland Ice Sheet ablation zone, a critical parameter needed to calculate subsurface absorption and scattering of transmitted radiation that to our knowledge has received no direct field study.

Measurements of k_att in snowpack and sea ice indicate three main sources of variation with relevance to geophysical applications. First, the magnitude of k_att is primarily controlled by ice microstructure (e.g. the size, shape, orientation, and number of air bubbles, ice grains, and cracks), via its control on k_sca (Dadic et al., 2013; Libois et al., 2013; Light et al., 2004, 2008). For the range of air bubble sizes (∼ 10⁻³–10⁻⁴) and ice grain sizes (∼ 10⁻¹–10⁻³) observed in glacier ice, k_sca is effectively independent of wavelength in the visible and near-infrared spectrum (Bohren, 1983; Dadic et al., 2013; Perovich, 1996). Spectrally, k_att is low in the near-ultraviolet and blue-green spectrum (∼ 250–600 nm), where k_abs is extremely low (< 10⁻⁸), and progressively higher for wavelengths > 600 nm, where k_abs rapidly increases up to its maximum value (∼ 10⁻²) at the far end of the solar spectrum (Warren and Brandt, 2008). Vertically, k_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 upward) 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 centimetres of depth is typically sufficient to reach the asymptotic regime where monochromatic k_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 (Frey et al., 2011), but this source of variation is not examined here.

In addition to experimental values obtained from measurements of light transmission in ice or snow, k_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:

$$\begin{array}{}\text{(2)}& m\left(\mathit{\lambda }\right)=m_{\mathrm{re}}\left(\mathit{\lambda }\right)-i{m}_{\mathrm{im}}\left(\mathit{\lambda }\right),\end{array}$$

where m_re is the real part of the complex refractive index (∼ 1.31 in the visible), m_im is the imaginary part, and

$$\begin{array}{}\text{(3)}& k_{\mathrm{abs}}^{\mathrm{ice}}\left(\mathit{\lambda }\right)={\displaystyle \frac{\mathrm{4}\mathit{\pi }}{\mathit{\lambda }}}{m}_{\mathrm{im}}\left(\mathit{\lambda }\right)\end{array}$$

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-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 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 solar radiation in ice in both 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., 1999a, b; Liston and Winther, 2005). However, theoretical values for k_att are rarely validated experimentally, and to our knowledge no such experimental values exist for near-surface glacier ice.

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 k_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 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 (ICESat) transmitted 1064 nm laser pulses to measure the distance (range) between the satellite and ice sheet surfaces (Schutz et al., 2005). Photons with a wavelength of 1064 nm penetrate into 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 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, similarly 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 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 k_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 k_att obtained from the two-stream analytical solution (see Eq. 26 in 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 k_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 measurement uncertainty, and in Sect. 5 we discuss how our k_att values differ from 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.

Figure 1 The irradiance sensor is comprised of a remote cosine receptor and fibre-optic light guide fitted inside an insulated white PVC tube of 2 m in length. Holes are drilled level and horizontal into the ice; the tube is inserted, and drill shavings are packed around the hole to prevent stray reflections. The cosine receptor collects the downwelling light and guides it to the fibre-optic cable that transmits the light to an Ocean Optics^® JAZ spectrometer, and a computer running the Ocean Optics^® OceanView software records the spectra. Background downwelling surface spectra are recorded on a 2 m mast drilled into the ice approximately 3 m to the northwest of the in-ice measurement location (see photo background). This photograph was taken on 21 July 2018 at ∼ 13:22 local time (UTC−3).

2 Methods

2.1 Transmittance measurements

Ice transmittance was measured on 20 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 above sea level. (67.15^∘ N, 50.02^∘ W). Subsurface (in-ice) spectral irradiance was measured at a ∼ 0.35 nm spectral resolution in the wavelength range of 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 consisting of a 90^∘ collimating lens adapter and a remote cosine receptor (RCR) with a Spectralon^® 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^∘ away from the sun and at a 2.5 m horizontal distance from the measurement location.

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 and levelled horizontally with a digital spirit level, and the sequence was repeated to a 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 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 a 0.0228 s integration time per scan, yielding a 0.46 s total integration time per irradiance measurement. Irradiance measurements were recorded at a 1 Hz frequency for 30 s yielding 30 irradiance profiles at each depth, after which the tube was removed; the next hole was drilled, and the sequence was repeated, working from the top hole toward the bottom on 20 July and from the bottom hole toward the top on 21 July. The measurements were completed between 13:45 and 14:35 local time (UTC−3) on 20 July and between 13:09 and 14:00 on 21 July, at solar zenith angles of ∼ 48–51^∘. Solar noon at this time and location is ∼ 13:26.

Background upwelling and downwelling spectral irradiances were measured continuously at 2 m height above the ice surface ∼ 3 m away from the in-ice measurements with a dual-channel Ocean Optics JAZ spectrometer. These data were recorded at a 1 min frequency using a 30-scan average with a 0.011 s integration time. Light was guided to the spectrometer via two 3 m fibre-optic cables attached to two RCRs mounted in an 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.

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 and 21 July immediately 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 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 experimental setup as described above and with the in-ice spectrometer held level above the ice surface in an upward-looking orientation. These residual dark-current spectra are treated as systematic errors and are subtracted from all irradiance profiles as a pre-processing step prior to analysis (Fig. 2a).

Figure 2 (a) Residual dark-current spectra for the surface-based reference spectrometer and the in-ice spectrometer. Dark-current spectra are recorded prior to each absolute irradiance measurement as input to the OceanView software dark-current correction module. Shown here are residual dark-current spectra after automated software correction, which are treated as systematic errors and subtracted from irradiance profiles prior to fitting experimental k_att values. (b) Ideal angular response function (ideal cosine) and empirical angular response function provided by Ocean Optics from laboratory measurements on the same type of irradiance sensor used in this study. The dashed red line in (b) is used as an empirical probability density function for the angular response of the cosine receptor in our Monte Carlo simulations of detector interference.

Figure 3 Photographs of an ice core collected at the field site. (a) The upper few centimetres of ice is semi-granular, with ∼ 4 cm of unrecovered granular ice not shown. (b) The 122 cm ice core was broken into three segments corresponding to depths of 4–45, 45–74, and 74–122 cm below the ice surface (the far right of the image in b is at 74 cm). The density of these segments is 801, 884, and 888 kg m⁻³, respectively. Below ∼ 10 cm, the bubbly ice appears foliated, indicating variations in bubble density and size distribution that affect scattering. Black box in (b) is approximately the image area in (a).

2.2 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 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 10 in-ice irradiance profiles for that measurement, prior to the rapid increase in light intensity, and discarding the remainder.

2.3 Ice thickness and density

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

$$\begin{array}{}\text{(4)}& \mathrm{\Delta }z=\mathrm{\Delta }h{\displaystyle \frac{\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{ice}}}},\end{array}$$

where Δh is in situ ice thickness between detector positions, ρ is in situ ice density, ρ_ice is solid ice density (917 kg m⁻³), and Δz is solid ice thickness between detector positions. Two separate observers made 10 independent measurements of Δh. In addition, one observer made 41 replicate measurements of an ablation stake using the same metre stick, yielding a mean difference and standard error in Δh. The ice density ρ was measured on a 1.2 m ice core extracted at the measurement location with a Kovacs Mark IV corer (https://kovacsicedrillingequipment.com/, last access: 7 April 2021) (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.

2.4 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):

$$\begin{array}{}\text{(5)}& I_z\left(\mathit{\lambda }\right)=I_{\mathrm{0}}\left(\mathit{\lambda }\right)\exp \left[-k_{\mathrm{att}}\left(\mathit{\lambda }\right)(z-z_{\mathrm{0}})\right],\end{array}$$

where k_att(λ) is the asymptotic flux attenuation coefficient, I_z(λ) is in-ice spectral irradiance at depth z, I₀(λ) is background downwelling spectral irradiance, z₀ is the ice surface, and

$$\begin{array}{}\text{(6)}& T_z\left(\mathit{\lambda }\right)=I_z\left(\mathit{\lambda }\right)/I_{\mathrm{0}}\left(\mathit{\lambda }\right)\end{array}$$

is spectral transmittance. The optical depth τ_z(λ) is a dimensionless path length that scales the physical thickness of a layer by its attenuation rate:

$$\begin{array}{}\text{(7)}& {\mathit{\tau }}_z\left(\mathit{\lambda }\right)=-\ln T_z\left(\mathit{\lambda }\right)=k_{\mathrm{att}}\left(\mathit{\lambda }\right)\left(z-z_{\mathrm{0}}\right).\end{array}$$

Estimates of k_att(λ) for each spectral band are obtained by solving a linear equation of the form

$$\begin{array}{}\text{(8)}& {\mathit{\tau }}_z\left(\mathit{\lambda }\right)={\mathit{\tau }}_{\mathrm{0}}\left(\mathit{\lambda }\right)+k_{\mathrm{att}}\left(\mathit{\lambda }\right)(\mathrm{\Delta }z+{\mathit{\epsilon }}_{\mathrm{\Delta }z})+{\mathit{\epsilon }}_{\mathrm{\Delta }\mathit{\tau }},\end{array}$$

where τ₀ is a parameter (y intercept), $\mathrm{\Delta }z=z-z_{\mathrm{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 present in both the independent and the dependent variables (see Sect. 2.9) (York et al., 2004).

The attenuation length l_att(λ) is the inverse of k_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 % $(\mathrm{1}/e)$ of its incident intensity, i.e. the path length at which $T=\mathrm{1}/e$ and τ=1:

$$\begin{array}{}\text{(9)}& l_{\mathrm{att}}\left(\mathit{\lambda }\right)={\displaystyle \frac{\mathrm{1}}{k_{\mathrm{att}}\left(\mathit{\lambda }\right)}}.\end{array}$$

The ice surface spectral albedo is the ratio of the upwelling spectral irradiance to the downwelling spectral irradiance:

$$\begin{array}{}\text{(10)}& \mathit{\alpha }\left(\mathit{\lambda }\right)={\displaystyle \frac{I\uparrow \left(\mathit{\lambda }\right)}{I\downarrow \left(\mathit{\lambda }\right)}},\end{array}$$

and the broadband albedo is

$$\begin{array}{}\text{(11)}& \mathit{\alpha }=\underset{{\mathit{\lambda }}_{\mathrm{1}}}{\overset{{\mathit{\lambda }}_{\mathrm{2}}}{\int }}\mathit{\alpha }\left(\mathit{\lambda }\right)I_{\mathrm{0}}\left(\mathit{\lambda }\right)\text{d}\mathit{\lambda }/\underset{{\mathit{\lambda }}_{\mathrm{1}}}{\overset{{\mathit{\lambda }}_{\mathrm{2}}}{\int }}{I}_{\mathrm{0}}\left(\mathit{\lambda }\right)\text{d}\mathit{\lambda }.\end{array}$$

2.5 Asymptotic flux attenuation coefficients

Theoretical k_att(λ) values are calculated using the asymptotic solution to the δ-Eddington two-stream radiative transfer approximation (Joseph et al., 1976; Schuster, 1905):

$$\begin{array}{}\text{(12)}& k_{\mathrm{att}}\left(\mathit{\lambda }\right)={\displaystyle \frac{\mathrm{3}}{\mathrm{4}}}{\displaystyle \frac{Q_{\mathrm{ext}}\left(\mathit{\lambda }\right)}{r_{\mathrm{eff}}}}\sqrt{\mathrm{3}\left(\mathrm{1}-\stackrel{\mathrm{‾}}{\mathit{\omega }}\left(\mathit{\lambda }\right)\right)\left(\mathrm{1}-g\left(\mathit{\lambda }\right)\stackrel{\mathrm{‾}}{\mathit{\omega }}\left(\mathit{\lambda }\right)\right)},\end{array}$$

where Q_ext(λ) is the extinction efficiency; r_eff is the effective scattering particle radius (m); g(λ) is the average cosine of the scattering angle, also referred to as the asymmetry parameter; and $\stackrel{\mathrm{‾}}{\mathit{\omega }}\left(\mathit{\lambda }\right)$ is the single-scattering albedo,

$$\begin{array}{}\text{(13)}& \stackrel{\mathrm{‾}}{\mathit{\omega }}\left(\mathit{\lambda }\right)={\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{sca}}\left(\mathit{\lambda }\right)}{{\mathit{\sigma }}_{\mathrm{att}}\left(\mathit{\lambda }\right)}},\end{array}$$

where σ_att(λ) and σ_sca(λ) are the single-scattering attenuation coefficient (m⁻¹) and scattering coefficient (m⁻¹), 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).

To estimate r_eff, Eq. (12) is inverted and solved by iteration for the value of r_eff that minimizes the difference between measured and calculated k_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, r_eff would be overestimated. Values for Q_ext(λ), g(λ), and $\stackrel{\mathrm{‾}}{\mathit{\omega }}\left(\mathit{\lambda }\right)$ are obtained from Mie scattering algorithms provided as MATLAB^® code (Mätzler, 2002) with input m(λ) from Warren and Brandt (2008). The Mie solutions at each wavelength are integrated over a Gaussian size distribution (N=1000) of scattering radii $N\left({\mathit{\mu }}_r=r_{\mathrm{eff}};{\mathit{\sigma }}_r=\mathrm{0.15}{r}_{\mathrm{eff}}\right)$ to eliminate ripples associated with Bessel function solutions to the Mie equations (Gardner and Sharp, 2010). The optimal r_eff values are ∼ 9.3 and ∼ 10.6 mm with corresponding specific surface areas ∼ 0.35 and ∼ 0.31 m² kg⁻¹ for 20 and 21 July experimental values, respectively.

2.6 Flux absorption coefficients

Warren et al. (2006) developed a method to estimate k_abs for pure ice (i.e. $k_{\mathrm{abs}}^{\mathrm{ice}}$) from measurements of flux attenuation in snow in Antarctica. The method relies on three assumptions: (1) the value of $k_{\mathrm{abs}}^{\mathrm{ice}}$ at the reference wavelength (λ₀=600 nm) is known accurately, (2) the value of k_att at λ₀ is not affected by LAPs in the measured snow or ice, and (3) $\stackrel{\mathrm{‾}}{\mathit{\omega }}\left(\mathit{\lambda }\right)$ varies so little as to be effectively independent of wavelength in the spectral range considered (here the near-UV and visible spectra). Warren et al. (2006) verified the validity of these assumptions for the spectral range of 350–600 nm and obtained the following relation (Eq. 15 of that paper) between flux attenuation and flux absorption:

$$\begin{array}{}\text{(14)}& {\left[{\displaystyle \frac{k_{\mathrm{att}}\left(\mathit{\lambda }\right)}{k_{\mathrm{att}}\left({\mathit{\lambda }}_{\mathrm{0}}\right)}}\right]}^{\mathrm{2}}\approx \left[{\displaystyle \frac{k_{\mathrm{abs}}\left(\mathit{\lambda }\right)}{k_{\mathrm{abs}}\left({\mathit{\lambda }}_{\mathrm{0}}\right)}}\right].\end{array}$$

Equation (14) assumes that k_abs is not affected by LAPs at the reference wavelength (600 nm), but the relation can be used to estimate k_abs at any other wavelength, including those where absorption is affected by LAPs. At those wavelengths, Eq. (14) will predict values for k_abs higher than $k_{\mathrm{abs}}^{\mathrm{ice}}$ if LAPs are present in the measured snow or ice volume, due to the influence of LAPs on k_att.

The inferred k_abs values can be related to a mass absorption cross section (MAC) (Doherty et al., 2010):

$$\begin{array}{}\text{(15)}& k_{\mathrm{abs}}\left(\mathit{\lambda }\right)=k_{\mathrm{abs}}^{\mathrm{ice}}\left(\mathit{\lambda }\right)+\mathit{\beta }c{\mathit{\rho }}_i,\end{array}$$

where β is the spectral MAC (m² kg⁻¹) and c is the mass mixing ratio of LAPs in the ice volume (g LAPs g⁻¹ ice). We exploit this to interpret differences between our theoretical and experimental values of k_att on the basis of differences between $k_{\mathrm{abs}}^{\mathrm{ice}}$ (Warren et al., 2006) and the k_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 c_eqBC and c_eqdust (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 $k_{\mathrm{abs}}^{\mathrm{ice}}$ and our inferred k_abs values at a reference wavelength, and it follows a similar approach to that used to infer LAP absorption in snowpack (Tuzet et al., 2019). For β_BC, we use 6 m² g⁻¹ as a representative MAC at 550 nm and an absorption Ångström exponent range of 0.8–1.9 to scale this value to 400 nm (Doherty et al., 2010). For β_dust, we use 0.013 m² g⁻¹ at 550 nm (Di Mauro et al., 2017) and an absorption Ångström exponent range of 2–3 (Doherty et al., 2010). We note that these descriptive estimates provide context for discussion; actual LAP species concentrations were not measured.

2.7 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 k_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., 2013; Light et al., 2008; Mullen and Warren, 1988). To account for non-diffuse near-surface attenuation, we define a piecewise optical depth:

$$\begin{array}{}\text{(16)}& \mathit{\tau }\left(\mathit{\lambda }\right)=\underset{\mathrm{0}}{\overset{z^{\prime }}{\int }}{k}^{\prime }\left(\mathit{\lambda }\right)\text{d}z+\underset{z^{\prime }}{\overset{z}{\int }}{k}_{\mathrm{att}}\left(\mathit{\lambda }\right)\text{d}z,\end{array}$$

where k^′ 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 k^′ with a centred finite-difference form of Eq. (7):

$$\begin{array}{}\text{(17)}& k^{\prime }\left(\mathit{\lambda }\right)=-{\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }{z}^{\prime }}}\ln \left[{\displaystyle \frac{I_{z^{\prime }}\left(\mathit{\lambda }\right)}{I_{\mathrm{0}}\left(\mathit{\lambda }\right)}}\right].\end{array}$$

Here, Δz^′ is 12 cm and $I_{z^{\prime }}$ 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 d_λ to include the effect of near-surface attenuation. The attenuation length is the depth at which τ=1. Setting τ(λ) to 1 in Eq. (16) and solving for z yields

$$\begin{array}{}\text{(18)}& z={\displaystyle \frac{\mathrm{1}-\mathrm{\Delta }{z}^{\prime }\left[k^{\prime }\left(\mathit{\lambda }\right)-k_{\mathrm{att}}\left(\mathit{\lambda }\right)\right]}{k_{\mathrm{att}}\left(\mathit{\lambda }\right)}}=d_{\mathit{\lambda }}.\end{array}$$

Equation (16) gives estimates of spectral transmittance that account for non-diffuse near-surface attenuation but relies on knowledge of k^′, which is sensitive to the spectral composition and directional distribution of I₀ 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:

$$\begin{array}{}\text{(19)}& {\mathit{\chi }}_{\mathrm{0}}=\underset{{\mathit{\lambda }}_{\mathrm{1}}}{\overset{{\mathit{\lambda }}_{\mathrm{2}}}{\int }}{I}_{\mathrm{0}}\left(\mathit{\lambda }\right)\exp \left[-{k^{\prime }}_{\mathit{\lambda }}\mathrm{\Delta }{z}^{\prime }\right]\text{d}\mathit{\lambda }/\underset{{\mathit{\lambda }}_{\mathrm{1}}}{\overset{{\mathit{\lambda }}_{\mathrm{2}}}{\int }}{I}_{\mathrm{0}}\left(\mathit{\lambda }\right)\text{d}\mathit{\lambda }.\end{array}$$

The χ₀ parameter is analogous to the i₀ 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:

$$\begin{array}{}\text{(20)}& i_{\mathrm{0}}=\underset{{\mathit{\lambda }}_{\mathrm{1}}}{\overset{{\mathit{\lambda }}_{\mathrm{2}}}{\int }}\left[\mathrm{1}-{\mathit{\alpha }}_{\mathit{\lambda }}\right]I_{\mathrm{0}}\left(\mathit{\lambda }\right)\exp \left[-{k^{\prime }}_{\mathit{\lambda }}\mathrm{\Delta }{z}^{\prime }\right]\mathrm{d}\mathit{\lambda }/\underset{{\mathit{\lambda }}_{\mathrm{1}}}{\overset{{\mathit{\lambda }}_{\mathrm{2}}}{\int }}\left[\mathrm{1}-{\mathit{\alpha }}_{\mathit{\lambda }}\right]I_{\mathrm{0}}\left(\mathit{\lambda }\right)\mathrm{d}\mathit{\lambda }.\end{array}$$

The i₀ 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 i₀=70 % for the visible (200–700 nm) and i₀=0 % for the infrared (700–5000 nm) (Briegleb and Light, 2007). The important distinction is that i₀ partitions the absorbed flux whereas χ₀ partitions the downward flux (Brandt and Warren, 1993). For both χ₀ and i₀, we set Δz^′ to 10 cm for consistency with prior work (Grenfell and Maykut, 1977; Light et al., 2008; Maykut and Untersteiner, 1971).

2.8 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 k_att values, following methods developed to simulate light propagation in biological tissue, ocean waters, and sea 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 Q_ext(λ), g(λ), and $\stackrel{\mathrm{‾}}{\mathit{\omega }}\left(\mathit{\lambda }\right)$ inferred from our optical measurements (Sect. 2.5). A complete technical description is given in the Supplement, where model accuracy is verified by comparison with benchmark solutions to the radiative transfer equation (van de Hulst, 1980).

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 simulations 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 detector rod. The detector rod albedo ${\stackrel{\mathrm{‾}}{\mathit{\omega }}}_{\mathrm{rod}}\approx \mathrm{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 model is integrated over 10 000 interactions at nine wavelengths in 50 nm increments from 350 to 750 nm, allowing us to fit the wavelength dependence of the estimated systematic uncertainty in simulated k_att values.

Figure 4 (a) Field spectra of surface downwelling (z=z0) and in-ice irradiance at four depths below the ice surface collected on 20 July 2018 (Layer A) between 13:45 and 14:35 local time in the western Greenland ablation zone (67.15^∘ N, 50.02^∘ W). Raw data were recorded at a 1 Hz frequency for 30 s, yielding 30 irradiance profiles at each depth. Shown here are 30 s averages at a ∼ 0.35 nm spectral resolution for each depth (black dots) and 1 nm interpolated values smoothed with a 3 nm centred moving mean filter for clarity (continuous lines). (b) Transmittance at each depth, with 30 s averages (black dots), 1 nm interpolated values (continuous lines), and shaded bounds (±2σ) representing propagated measurement uncertainty deduced from the standard deviations of the 1 Hz raw data (N=30 for each value). Results for the 21 July 2018 (Layer B) experiment (not shown) are qualitatively similar.

2.9 Uncertainty propagation

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 s_i and sample size N≥30, the standard error is $s_i/\sqrt{N}$. For N<30, standard errors are scaled by a critical t value drawn from 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 I(λ) is

$$\begin{array}{}\text{(21)}& {\mathit{\sigma }}_{\text{I}}=\sqrt{{\left({\mathit{\sigma }}_{\text{I}}^{\ast }\right)}^{\mathrm{2}}+{\left({\mathit{\sigma }}_{\text{D}}\right)}^{\mathrm{2}}},\end{array}$$

where ${\mathit{\sigma }}_{\text{I}}^{\ast }$ is the standard deviation of the high-frequency irradiance spectra before dark-noise correction and σ_D 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 T(λ) is

$$\begin{array}{}\text{(22)}& {\mathit{\sigma }}_T=T\sqrt{{\left({\displaystyle \frac{{\mathit{\sigma }}_{I_z}^{\mathrm{cal}}}{I_z^{\mathrm{cal}}}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{{\mathit{\sigma }}_{I_{\mathrm{0}}}}{I_{\mathrm{0}}}}\right)}^{\mathrm{2}}},\end{array}$$

where ${\mathit{\sigma }}_{I_z}^{\mathrm{cal}}$ and ${\mathit{\sigma }}_{I_{\mathrm{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

$$\begin{array}{}\text{(23)}& {\mathit{\sigma }}_{\mathit{\tau }}={\displaystyle \frac{{\mathit{\sigma }}_T}{T}},\end{array}$$

and the combined uncertainty for k_att is

$$\begin{array}{}\text{(24)}& {\mathit{\sigma }}_k=\sqrt{({\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\tau }}{)}^{\mathrm{2}}+({\mathit{\sigma }}_{\mathrm{\Delta }z}{)}^{\mathrm{2}}}.\end{array}$$

Equation (24) gives a first-order description of σ_k due to statistical propagation of measurement uncertainties, neglecting higher-order interaction terms. A description of the statistical uncertainty in fitted k_att(λ) values is given by the MLE of the regression slope of Eq. 8, which can be expressed in terms of an error model as

$${\widehat{\mathit{\tau }}}_z\left(\mathit{\lambda }\right)={\mathit{\tau }}_{\mathrm{0}}\left(\mathit{\lambda }\right)+k_{\mathrm{att}}\left(\mathit{\lambda }\right)\left(\mathrm{\Delta }\widehat{z}+{\mathit{\epsilon }}_{\mathrm{\Delta }\mathrm{z}}\right)+{\mathit{\epsilon }}_{\mathrm{\Delta }\mathit{\tau }},$$

where $\widehat{\mathit{\tau }}$ and $\mathrm{\Delta }\widehat{z}$ are the true but unobserved (due to measurement error) optical depth and ice thickness and ${\mathit{\epsilon }}_{\mathrm{\Delta }z}\sim N\left(\mathrm{0},{\mathit{\sigma }}_{\mathrm{\Delta }z}^{\mathrm{2}}\right)$ and ${\mathit{\epsilon }}_{\mathrm{\Delta }\mathit{\tau }}\sim N\left(\mathrm{0},{\mathit{\sigma }}_{\mathrm{\Delta }\mathit{\tau }}^{\mathrm{2}}\right)$ 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 k_att(λ) are adjusted for N−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 k_att(λ) values.

Figure 5 (a) Sample linear regressions between measured transmittance (indicated by rectangles with width and height proportional to measurement uncertainty in both variables) and depth in the range of 12–77 cm (Layer A) at five representative wavelengths spanning the measured spectral range. The slope of each line is the attenuation coefficient k_att. Shaded bounds are ±95 % confidence intervals from a two-sided t distribution. (b) Red box inset in (a) shows the y-axis intercept of each regression is less than 100 %, indicating the magnitude of deviation from Bouguer's law near the surface. (c) Spectral k_att (blue dots with shaded uncertainty; left axis) and spectral albedo (red dots with shaded uncertainty; right axis). Beyond ∼ 700 nm, in-ice transmitted irradiance is too low to reliably estimate k_att (see Figs. 4a and 2c), as indicated by the increased scatter in k_att values.

3 Results

3.1 Spectral transmittance

Four in-ice irradiance spectra were collected at 12, 36, 58, and 77 cm depth below the ice surface on 20 July (hereafter referred to as Layer A) (Fig. 4a) and at 53, 67, 82, and 124 cm on 21 July (hereafter referred to as Layer B). At all depths, spectral transmittance T is maximal at 430 nm and maintains relatively stable and high values at wavelengths up to ∼ 500 nm in the visible, beyond which T decreases into the red end of the visible spectrum following the well-known exponential increase in ice absorptivity (Fig. 4c). Maximum T values vary from 78 % at 12 cm to 45 % at 77 cm. For wavelengths > 500 nm, T 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 deviation of the 1 Hz raw data is < 1 W m⁻² nm⁻¹ at all wavelengths, consistent with field observations of thick cloud cover and diffuse light conditions described in Sect. 2.2. Instrumental noise and high-frequency measurement variations propagate as ±1.6 % uncertainty on T 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.

Figure 6 Visible and near-UV irradiance attenuation coefficient k_att spectra from measurements of light transmission in bare glacier ice collected on 20 July 2018 at 12–77 cm depth below the ice sheet surface. Field estimates are compared with asymptotic two-stream theory for optically clean bubbly ice (continuous dotted line) and with values at nine wavelengths from four simulations with a three-dimensional Monte Carlo radiative transfer model (solid symbols). Monte Carlo values for clean bubbly ice are shown for 350, 400, and 450 nm to demonstrate detector interference at these wavelengths; values at wavelengths > 550 nm converge with field spectra and are omitted for clarity. Two measures of uncertainty are shown: (1) statistical linear model uncertainty ε_LM (shaded uncertainty bounds; ±1 standard error in the linear regression) and (2) ε_LM combined with systematic uncertainty ε_MC due to detector interference estimated with Monte Carlo simulation (error bars; μ±ε). The same comparison for the 21 July experiment (inset; Layer B) suggests detector interference is within statistical uncertainty at wavelengths > 400 nm.

3.2 Experimental flux attenuation coefficients and albedo

Fitted k_att values for Layer A range from 0.98 ± 0.17 to 7.86 ± 0.43 m⁻¹ for wavelengths between 350–750 nm (Fig. 5a), with uncertainty bounds that represent combined statistical and systematic uncertainty (see Sect. 4 for a discussion of systematic error). These values correspond to attenuation lengths of 1.02 ± 0.18 to 0.13 ± 0.007 m, respectively. Layer B k_att values are ∼ 12 % lower than Layer A values at 350–500 nm and within 1 % at 650 nm (see inset Fig. 6). For Layer A, the minimum in k_att is at 390 nm, blue-shifted relative to the maximum in T at 430 nm. For Layer B, the minimum is at 397 nm. The coefficient of determination (r²) 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 and 77 cm depth were too low to reliably fit k_att values (see Figs. 4b and 5c). For Layer B values, low light levels prevented fits beyond ∼ 650 nm.

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 k_att and the maximum in T. 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).

3.3 Theoretical flux attenuation coefficients

Asymptotic k_att values predicted by two-stream theory for optically clean bubbly ice are nearly 1 order of magnitude lower 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 βcρ_i in Eq. 15) varies from 0.009–0.015 m⁻¹ at 350–530 nm. The equivalent black carbon concentration c_eqBC inferred at 400 nm is 1–2 ng g⁻¹ for both Layer A and Layer B, where the range covers uncertainty in both the absorption spectra and the absorption Ångström exponent (Doherty et al., 2010). The equivalent mineral dust concentration c_eqdust is ∼ 344–620 ng g⁻¹ for Layer A and 303–545 ng g⁻¹ for Layer B. Monte Carlo simulations without detector interference replicate both asymptotic theory for clean bubbly ice (i.e. when forced with $k_{\mathrm{abs}}^{\mathrm{ice}}$) and field estimates when forced with k_abs values inferred from our optical measurements (solid blue squares, Fig. 6). Monte Carlo simulations of detector interference are discussed further in Sect. 4.

Figure 7 (a) Effective attenuation coefficient k^′ for the near-surface 0–12 cm region compared to k_att values estimated for the interior 12–77 cm (Layer A) and 53–124 cm (Layer B) depth regions. (b) Effective k^′ values are ∼ 1.6 times larger than Layer A values at wavelengths beyond about 600 nm but are ∼ 2–4 times larger between 400–600 nm. The spectral dependence suggests a higher influence of absorptive impurities on attenuation enhancement near the ice surface than in the ice interior. The shaded bounds on k^′ represent propagated ±1.2 cm vertical measurement uncertainty.

3.4 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 k^′ values (Eq. 17) for the quasi-granular 0–12 cm layer are ∼ 1.5 times higher than k_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 I_z is overestimated by 10 %–60 % if Bouguer's law is applied to the incident downwelling irradiance I₀ using asymptotic k_att values, with median overestimation of 23 % (Fig. 8a). In contrast, the piecewise optical depth (Eq. 16) predicts I_z to within 12 % of measured values for all wavelengths between 350–700 nm with a median error of 3 %. Integrated over these wavelengths, χ₀ is 0.68 and i₀ is 0.66, suggesting 66 % of the total incoming irradiance within this spectral region is absorbed at depths below 10 cm. If k_att is used rather than k^′ to calculate χ₀ and i₀, the respective values are 0.81 and 0.79.

Stated in terms of penetration depth, d_eff varies from 12–84 cm between 350–700 nm. These values are 13 %–44 % lower than attenuation lengths l_att inferred from empirical asymptotic k_att values. Specifically at 532 nm, d_eff is 52 cm, or 10 cm lower than the 62 cm empirical l_att value and 14 cm lower than the 66 cm theoretical l_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, 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), 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; Light et al., 2008).

Figure 8 Measured in-ice irradiance compared with (a) Bouguer's law (Eq. 7) with no modification and (b) the piecewise Bouguer's law (Eq. 16). The error structure (c, d) provides insight into the near-surface attenuation processes: relative errors (%) are positive (model underpredicts attenuation) at all wavelengths but are highest in the near-UV, lowest in the blue, and increase monotonically into the red end of the visible spectrum. (d) Errors are small and generally decrease monotonically with increasing wavelength. Taken together, near-surface attenuation enhancement is ∼ 10 %–60 %.

4 Uncertainty analysis

The effect of random and systematic uncertainties on our optical measurements and fitted k_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.

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⁻³ 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 value and with ρ values estimated for each Δh value from linear and cubic interpolation of the vertical density profile. The maximum Δz difference was 0.9 cm. The k_att values differed by < 1 %, and r² values were nearly identical. We use the depth-weighted average ρ values to calculate Δz, which is 835 and 855 kg m⁻³ for the measurements collected on 20 and 21 July, respectively.

Figure 9 Uncertainty budget for reported asymptotic attenuation coefficient k_att values for the 12–77 cm depth region (Layer A). Systematic uncertainties examined include spectrometer sensitivity to dark noise, ice thickness (detector position) measurement uncertainty, the non-ideal angular response of the irradiance sensor, and scattering and absorption interference by the polyvinyl chloride detector rod (estimated with Monte Carlo simulation). These systematics are combined with statistical uncertainty represented by 1 standard error in the maximum likelihood estimation (MLE) linear regression fits. Values for Layer B (53–124 cm) (not shown) are qualitatively similar but lower, with total uncertainty < 14 % in the region of maximum uncertainty.

Detector position uncertainty was further assessed by fitting k_att values with an ensemble of 10 000 Δz values perturbed with random errors drawn from a normal distribution $N(\mathit{\mu }=\mathrm{0.9}\mathrm{cm};\mathit{\sigma }=\mathrm{0.4}\mathrm{cm})$. At all wavelengths, the chance of obtaining a fitted k_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 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 k_att values with and without residual dark-noise correction. The mean difference was $-\mathrm{0.01}\pm \mathrm{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 %.

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 of ∼ 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 presence of the detector rod both tend to increase k_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).

Figure 10 Attenuation coefficient spectra for seven distinct ice structures (from lower left clockwise): interior of clean, dry snowpack ($r_{\mathrm{eff}}\approx {\mathrm{10}}^{-\mathrm{5}}$ m) near Dome C in Antarctica for two depth regions (90–135 cm and 40–90 cm) (Warren et al., 2006); near-surface (0–10 cm) dry snowpack ($r_{\mathrm{eff}}\approx {\mathrm{10}}^{-\mathrm{3}}$ m) in Dronning Maud Land, Antarctica (Järvinen and Leppäranta, 2013); near-surface (0–5 cm) and interior (0–18 cm) dry snowpack ($r_{\mathrm{eff}}\approx {\mathrm{10}}^{-\mathrm{4}}$ m) near Summit, Greenland (Meirold-Mautner and Lehning, 2004); surface scattering layer (SSL; 0–5 cm) of multi-year sea ice in the Arctic Ocean (Grenfell and Maykut, 1977); interior of ablating glacier ice in Greenland (this study) (12–77 cm in solid line; 0–12 cm in dotted line); interior of first-year sea ice in the coastal zone near Barrow, Alaska (Grenfell et al., 2006); and interior of multi-year sea ice in the Arctic Ocean (Grenfell et al., 2006). 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.

5 Discussion

5.1 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 k_att for the spectral range 350–750 nm. Prior studies quantified k_att for sea ice and snowpack (e.g. Fisher et al., 2005; Frey et al., 2011; Gerland et al., 2000; Grenfell and Maykut, 1977; Järvinen and 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., 1995, 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 non-scattering clathrates and light scatters on dust grains (Price and Bergström, 1997b). In the bubbly-ice regime studied by AMANDA, k_sca values at 532 nm are ∼ 1–3 m⁻¹, comparable to the 1.6 m⁻¹ 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.

Figure 10 compares our k_att spectra for glacier ice to seven previously published spectra for snowpack and sea ice. In general, glacier ice is the most transparent structure examined with the exception of multi-year and first-year interior sea ice at wavelengths > 540 nm (Grenfell et al., 2006). 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 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 the highest of all, possibly due to direct light scattering in the near-surface optical boundary layer. The comparison demonstrates that k_att values vary by nearly 2 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 k_att are up to 1 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., 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⁻¹ (parts per billion, 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 of 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.

Figure 11 Depth profiles of South Pole ice absorption coefficient at four wavelengths obtained from Eq. (26) and Table 3 of Ackermann et al. (2006). The 1 % K⁻¹ temperature dependence of pure-ice absorptivity (Woschnagg and Price, 2001) is removed from South Pole values for comparison with this study's lower (Layer B) and upper (Layer A) absorption coefficient estimate at each wavelength (dashed lines). Values reported in this study are consistent with South Pole values at depths corresponding to the Last Glacial Maximum and Marine Isotope Stage 4 when atmospheric-dust concentrations peaked in both hemispheres. Note that South Pole age and dust concentration do not map to ice near the Greenland Ice Sheet margin. Rather, Southern Hemisphere dust concentrations during the LGM and MIS-4 are consistent with Northern Hemisphere dust concentrations during warm interglacial periods and/or periods with low aeolian activity (Muhs, 2013; Reeh et al., 2002).

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., 1999). Hemispherical dust fluxes are generally synchronous at these timescales; similar peaks at the 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 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).

Regarding pure-ice absorptivity, our k_abs values provide support for 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 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.6 ppb of black carbon (75^∘ S, 123^∘ E; 3230 m above sea level). Picard et al. (2016) (hereafter P16) repeated the W06 experiment on 56 transmitted radiance profiles collected in snow near Dome C with variable impurity content. The values shown in Fig. 12 are their best estimate of pure-ice absorptivity from radiance profiles collected in snow with low impurity content (see “clean” subset, Fig. 17 of that paper). P16 were unable to reconcile their values with W06, after considering published values for impurity loadings in the vicinity of Dome C, suggesting the W06 values were unreasonably low. Regardless of that discrepancy, our values were undoubtedly influenced by LAPs but are lower than the P16 values. Treating our Monte Carlo estimate of detector interference as a known systematic error would bring our values closer to the AMANDA and W06 values and further from the P16 values (see dotted line, Fig. 12).

Figure 12 Estimates of ice absorption coefficient k_abs, obtained from five distinct sources: laboratory-grown pure ice (Grenfell and Perovich, 1981; Perovich and Govoni, 1991), as compiled in Warren (1984); snow in Antarctica, contaminated by trace concentrations of light-absorbing particles (LAPs) (Picard et al., 2016); glacial ice in Greenland with unknown concentration of LAPs (this study); compressed glacial ice at 1755 m depth and 830 m depth in the Antarctic Ice Sheet contaminated by dust deposited during the Late Pleistocene and Early Holocene, respectively (Ackermann et al., 2006); and snow in Antarctica with the effect of LAPs removed (pure-ice estimate) (Warren et al., 2006). Values from this study with the detector interference subtracted are shown as dotted blue line.

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 coincident with our maximum in transmittance. A similar shift is apparent in all attenuation coefficient spectra shown in Fig. 10 that report the surface as a reference horizon but is absent in those that report an interior reference horizon. P16 used an interior reference horizon, excluding radiance measurements within 8 cm of the surface based on Monte Carlo simulations of detector interference and visual inspection of homogeneous attenuation zones, but their shifted minimum may indicate that radiance profiles in the near-UV and blue spectra were impacted disproportionately by detector rod interference and/or other near-surface effects. The same effects may explain the spectral structure in our near-surface (0–12 cm) effective attenuation coefficient profile between 350–400 nm and its broad minimum between 430–490 nm (dotted line, Fig. 10). Similar spectral structure is apparent in diffuse attenuation coefficients obtained in snowpack in the French Alps using the P16 method (cf. Fig. 3b in 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⁻² m⁻¹ and possibly lower (Ackermann et al., 2006; Picard et al., 2016; Warren et al., 2006).

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 $k_{\mathrm{abs}}^{\mathrm{ice}}$ to calculate k_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 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 k_abs values for contaminated ice given here and snowpack elsewhere (Picard et al., 2016) could provide realistic input for radiative transfer models absent explicit knowledge of realistic LAP concentrations. In contrast, simulations that use the canonical $k_{\mathrm{abs}}^{\mathrm{ice}}$ 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.

6 Conclusion

We report the first in situ spectral measurements of near-UV and visible light transmission in near-surface bare glacial ice, collected at a field site in the western Greenland ablation zone on 20–21 July 2018. In general, our empirical irradiance attenuation coefficients are nearly 1 order of magnitude larger in the range of 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⁻¹ at 390–397 nm, implying absorption length scales of 69–77 m. The volumetric scattering coefficient is 1.6 ± 0.2 m⁻¹ at 532 nm, with an asymptotic attenuation length scale of 0.62 ± 0.08 m. In addition to light scattering on air bubbles, we find that light attenuation is enhanced by a layer of quasi-granular white ice (weathering crust) 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 deeper 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 was absorbed within this upper 10 cm surface layer at this time and location, consistent with observations of the semi-granular surface layer on sea ice. The estimated absorption spectrum suggests equivalent black carbon and mineral dust concentrations consistent with pre-industrial and warm interglacial periods with low Northern Hemisphere aeolian activity, and therefore it may provide a reasonable lower bound on volumetric absorption enhancement due to impurities embedded in outcropping glacial ice in the western Greenland ablation zone.

Appendix A

Table A1 Estimates of attenuation coefficient and absorption coefficient obtained from solar irradiance transmission measurements in glacier ice (Layer B; 53–124 cm below the ice sheet surface), 1 standard error in the linear regression, and coefficient of determination.