The euphotic zone under Arctic Ocean sea ice: Vertical extents and seasonal trends

Eight Ice‐Tethered Profilers were deployed in the Arctic Ocean between 2011 and 2013 to measure vertical distributions of photosynthetically active radiation (PAR) and other bio‐optical properties in ice‐covered water columns, multiple times a day over periods of up to a year. With the radiometers used on these profilers, PAR could be measured to depths of only ∼20–40 m in the central Arctic in late summer under sea ice ∼1 m thick. At lower latitudes in the Beaufort Gyre, late summer PAR was measurable under ice to depths exceeding 125 m. The maximum depths of measurable PAR followed seasonal trends in insolation, with isolumes shoaling in the fall as solar elevation decreased and deepening in spring and early summer after insolation resumed and sea ice diminished. PAR intensities were often anomalously low above 20 m, likely due to a shading effect associated with local horizontal heterogeneity in light transmittance by the overlying sea ice. A model was developed to parameterize these complex vertical PAR distributions to improve estimates of the water column diffuse attenuation coefficient and other related parameters. Such a model is necessary to separate the effect of surface ice heterogeneity on under‐ice PAR profiles from that of the water column itself, so that euphotic zone depth in ice‐covered water columns can be computed using canonical metrics such as the 1% light level. Water column diffuse attenuation coefficients derived from such autonomously‐collected PAR profile data, using this model, agreed favorably with values determined manually in complementary studies.

The light field of an ice-covered ocean differs from that of an ice-free ocean in several fundamental ways. Most obviously, the presence of sea ice and surface snow strongly attenuates the transmission of incident sunlight into the water column below (Little et al. 1972;Grenfell and Maykut 1977). Sediment or algae in or on the ice layer will further attenuate this transmitted sunlight and may further alter its spectral distribution (Belzile et al. 2000). Optical heterogeneity in the ice layer itself and the varying presence of snow cover, between-floe leads, and surface melt ponds all introduce substantial horizontal variability in under-ice irradiance on scales of meters to tens of kilometers (Petrich et al. 2012;Katlein et al. 2016). These same factors also vary in time, generating temporal variability in under-ice light fields on scales of hours and longer. Sea ice is a highly scattering optical medium that affects the directionality of light fields under ice (Perovich 2003;Katlein et al. 2014), and sea ice can also alter the vertical distribution of light in the underlying water column, creating apparent subsurface maxima at locations where light attenuation by the ice layer above the sensor is stronger than it is in sea ice nearby (Frey et al. 2011). Even the diel character of water column light fields differs under sea ice, not due to the ice per se but because ice-covered oceans typically occur at high latitudes where the "normal" diel light cycles of lower latitudes occur for only a few weeks in spring and fall, between dark polar winter and the full-day insolation of polar summer.
These influences that sea ice has on light availability in the underlying water column will in turn affect spatial and temporal trends in primary production in pelagic ecosystems that experience ice cover either seasonally or perennially. These linkages between under-ice light fields and water column primary production are becoming increasingly relevant in the Arctic Ocean given the significant long-term decreases being observed there in both the areal coverage and thickness of sea ice (Tucker et al. 2001;Stroeve et al. 2007;Comiso et al. 2008;Kwok and Rothrock 2009;Laxon et al. 2013;Krishfield et al. 2014). Such changes in the sea ice layer act to both brighten and deepen the euphotic zone, providing more light to support primary production not only in waters now free of seasonal sea ice but also in regions that either seasonally or perennially still retain now-thinner sea ice cover (Perrette et al. 2011;Long et al. 2015). Such changes in Arctic sea ice could dramatically alter existing patterns of production and biogeochemistry in pelagic Arctic ecosystems in several ways: by advancing algal growth to earlier in the year (Leu et al. 2011), by decreasing the interval between the yearly ice algal bloom and the subsequent phytoplankton bloom in the underlying water column (Søreide et al. 2010), or by redistributing to water column phytoplankton much of the light utilization and production that currently occurs within the sea ice by ice algae. All three scenarios represent fundamental restructurings of primary production in the Arctic Ocean.
Technical and logistical challenges have long hindered any direct assessment of the timing of and seasonal trends in photosynthesis and production in ice-covered Arctic pelagic ecosystems. Physical-biological ecosystem models therefore play an especially important role in examining how such changes in sea ice may affect Arctic production and ecology (e.g., Arrigo et al. 2008;Zhang et al. 2010;Deal et al. 2011;Jin et al. 2012). The robustness of these models depends on having accurate representations of the vertical, spatial, seasonal, regional, and interannual variability of light fields under sea ice, but unfortunately direct observations of these modes of variability in under-ice light fields are broadly lacking for the Arctic. Autonomous systems are well suited for obtaining direct observations of irradiance in ice-covered water columns on the spatiotemporal scales that are needed to better constrain such models. Yet to date, key properties of ice-covered water columns such as the light field or euphotic zone depth have received little attention in long-term Arctic Ocean observing programs. This represents a significant gap in international monitoring efforts in these sensitive and rapidly changing ecosystems (National Research Council 2006;National Science Foundation 2007). Motivated by the ongoing success with using Ice-Tethered Profilers (ITPs) for long-term hydrographic assessments in ice-covered regions of the Arctic Ocean (Toole et al. 2006;Krishfield et al. 2008;Timmermans et al. 2010;Toole et al. 2011), and by progress with integrating bio-optical sensors into autonomous profilers for use in ice-free oceans (e.g., Boss et al. 2008;Bishop and Wood 2009), a prototype biooptical sensor suite was developed specifically for use with ITPs in ice-covered water columns in the Arctic . Eight ITPs were outfitted with these sensor suites and were deployed in the Arctic Ocean between 2011 and 2013 to obtain high-resolution profiles of optical and bio-optical properties, including PAR irradiance, under sea ice. These observational data, collected over seasonal to annual time scales with temporal resolutions of daily or better, facilitated an improved characterization of seasonal and sub-seasonal variability in the intensities and depths of the euphotic zone in ice-covered Arctic Ocean ecosystems.

Instrumentation and deployments
Long-term time series of under-ice irradiance profiles were collected in the Arctic Ocean using a variant of the ITP (Toole et al. 2006;Krishfield et al. 2008). In this variant, the standard conductivity-temperature-depth (CTD) sensor on the ITP's profiling vehicle was augmented with a custom sensor package (a "biosuite") consisting of a cosine radiometer for the photosynthetically active radiation (PAR) wavelengths between 400 nm and 700 nm (PAR-LOG-d, Satlantic) and a three-channel active sensor (FLBBCD, WETLabs) that measured the fluorescences of chlorophyll and dissolved organic matter and the optical backscattering at 700 nm. Both of these sensors were affixed to the top endcap of the profiler, oriented upward, with their optical faces at the same height as the topmost end of the CTD assembly (see Laney et al. 2014). This placement minimized potential artifacts in PAR that might arise from shading by the profiler body itself. A copper shutter that covered the faces of both sensors was opened only during profiling to provide physical protection for these optical sensors and to minimize potential biofouling in this already low-fouling environment. Each PAR datum in these profiles represents the average of 10 digitized acquisitions sampled in a burst, with bursts performed continuously at 4 Hz throughout the entire profile, providing an effective vertical resolution of these PAR data of 25-30 cm.
The depths over which these ITPs profiled were determined by the length of their tethers ($760 m) and by a physical safety stop on the upper end of each tether typically at $6 m. The profilers were programmed to finish every upward profile at 7 m but the safety stop further minimized chances that the profiler would damage itself or its upward facing sensors by impacting hardware affixed to the bottom of the 5 m urethane-reinforced section of the tether that extended through the ice. This safety stop also unavoidably limited the profiler from measuring irradiance at depths above 7 m, i.e., fully to the ice-ocean interface. All ITPs in this study were programmed to conduct a one-way vertical profile four times daily either to or from a pre-programmed depth of 200 m. For every third profiling cycle, this depth was increased to 750 m to provide periodic sampling deeper into the upper ocean. During November to February inclusive, this 6 h interval between one-way profiles was tripled in order to extend these profilers' operational lifetimes. By sampling in the upper 200 m primarily and by reducing sampling effort over wintertime, an ITP deployed in late summer can reserve adequate power to resume high-resolution profiling at 6 h intervals the following spring, providing it physically survives ice ridging, floe degradation, damage of the surface buoy by polar bears, and other modes of failure. Profiles were programmed to occur in synchrony with Universal Time and were not continually adjusted to coincide with local solar noon throughout each ITP's drift trajectory although this is technically feasible.
Eight such ITPs were deployed in the Arctic Ocean in late summer between 2011 and 2013 within different Arctic Ocean observing efforts (Table 1; Fig. 1). All ITPs were deployed with their surface buoys placed directly on the snow or ice surface and not on pallets as was the practice with earlier versions as described in Krishfield et al. (2008). Five of these eight systems were deployed in the Beaufort Gyre as part of the Beaufort Gyre Exploration Project (BGEP; www.whoi.edu/beaufortgyre). The first system (ITP52) was deployed in August 2011 slightly north of 788N after which it traveled a total of 824 km over 110 d, ceasing profiler transmissions in late November presumably when the ice floe on which it was deployed was reconfigured (as indicated by other instruments on the same floe ceasing to function at around the same time). The second system (ITP65) was deployed in August 2012 at nearly 818N and travelled 2474 km over the subsequent 307 d. The third system (ITP64) was also deployed in August 2012, further south in an area where ice floes were not suitable for safe through-ice deployment of an ITP. This third system was instead released directly into open water among broken pack ice where it subsequently survived fall freeze-up and operated for a total of 361 d along a drift track of 3147 km. It ceased receiving at all from its bio-optical sensor suite and thus did not contribute to subsequent analyses. The remaining three ITP systems were deployed in the central Arctic Ocean over the same 3 yr in two additional observational programs. ITP48 was deployed in September 2011 at $858N as part of the Hybrid Arctic/Antarctic Float Observation System program (HAFOS; Fahrbach and Boebel 2004). It operated for 437 d over a total distance of 2745 km before communications with the underwater unit ceased in mid-November 2012, presumably after grounding its tether in shallower waters north of Ellesmere Island. The PAR sensor on this ITP provided valid data for only the first 239 profiles. A second HAFOS system (ITP60) was deployed in September 2012 and the profiler communicated for 106 d over 1128 km before ceasing transmissions, but none of its profiles included acceptable PAR data. The third system (ITP72) was deployed in late August 2013 in collaboration with the Nansen and Amundsen Basins Observational System (NABOS; http://nabos.iarc.uaf.edu/) project and drifted 1150 km over 108 d before ceasing to transmit profiler data in mid-December 2013. ITP72 is an exception to the upper limit of 7 m set for profiling, as this deployment involved additional sensors placed on the tether at depths down to $12 m which set the upper limit of ITP72 profiles to $12 m.

Data set preparation
These PAR profiles were subjected to a quality control examination in order to prepare a working data set for use in subsequent analyses (Fig. 2). We first identified and removed samples measured during the 2-min sensor characterization cycles that were performed before and after each profile, when the optical sensors were powered while the antifouling shutter remained closed (see Laney et al. 2014). Profiles were then examined visually for spurious data and such profiles were omitted from the data set. This is a necessarily subjective step given the types of failures that can occur in long-term polar deployments given, for example, the lower reliability of subsea cable connectors in such cold waters. Such spurious data are best identified by visual inspection. Next, calibrations were applied and we removed any PAR measurements corresponding to intensities less than these sensors' threshold of detection. This threshold was determined empirically from measured PAR profiles and was found to be 0.0552 lE m 22 s 21 for five of the six ITPs that returned PAR data and slightly higher for the sixth (ITP52, 0.0561 lE m 22 s 21 ). These empirically determined thresholds are close to a theoretical threshold sensitivity of 0.05 lE m 22 s 21 that can be computed from the typical maximal measurable irradiance of these sensors (specified at 5000 lE m 22 s 21 ) and by assuming a five order of magnitude dynamic range.
The data set at this level of quality control was used for visualization and for general qualitative assessment of spatiotemporal characteristics and trends in under-ice PAR. For more quantitative subsequent analyses, individual profiles were additionally flagged but not removed from the data set when either of the following two conditions was met. First, a flag was given to any profile that occurred when local solar elevation at the time of sampling was expected to be zero or negative, as computed using a solar elevation model based on date, time of day, longitude, and latitude (Kirk 1994). This step identified those PAR profiles that occurred when local insolation was effectively zero, i.e., during winter or nighttime during the spring and fall. Second, a flag was given to any profile that did not contain PAR data from above the 15 m depth horizon and at least 10 vertical meters of continuous PAR samples. The 15 m minimum sampling horizon identified situations where the profiler vehicle may have not completed a full vertical excursion, as sometimes happens when the sea ice is drifting more than $30 cm s 21 relative to the underlying water column, or in open water situations where wave action is substantial. Profiles so flagged were excluded from the PAR profile modeling analysis described below to avoid artefactual results from incorporating profiles with inappropriate data (Table 2). This horizon sits below the topmost profiling depth of ITP72 ($12 m) and therefore allows us to retain this time series in our modeling analysis.

Spatiotemporal coverage and data quality
A total of 3269 PAR profiles passing quality control steps for visualization purposes were obtained from the six ITPs that transmitted data for at least one profile: 478 from the central Arctic and 2791 from the Beaufort Gyre (Fig. 3a, thin grey traces). Because ice-or instrument-related failures with ITPs are more likely to occur during the overwintering period, the majority of PAR profiles in these six time series represent the late-summer, fall, and early winter months following an ITP deployment. Sampling for PAR was not suspended during winter and thus many of these nominally "valid" PAR profiles were obtained when the sun was consistently below the horizon and were flagged accordingly. Of the 3269 nominally valid profiles only $48% (101 in the central Arctic Ocean and 1471 in the Beaufort Gyre) occurred at locations and times when the sun was above the horizon and under-ice PAR was likely nonzero (Fig. 3a, thick black traces). Of these six ITPs only ITP64 provided substantial PAR data throughout the subsequent summer following its deployment.

PAR: penetration depths and vertical distributions
The maximal depth to which PAR could be measured by these ITPs (i.e., the PAR penetration depth) is operationally defined, set by the minimum threshold of detection of the specific radiometers used on these profilers (see Methods, above). Therefore, the absolute values of these penetration depths have no ecological meaning, but relative differences in PAR penetration depths among and within these six time series can be used to quantify regional and temporal variability in the euphotic zones of these ice-covered Arctic Ocean water columns. For purposes of visualizing changes in the penetration depths among and within these time series, the one profile daily that occurred closest to solar noon was identified in each time series, serving as a daily representative subset of profiles. All other factors aside, sunlight would penetrate deepest into the ocean at the time of day when solar elevation is highest. However, because these ITPs were programmed to profile in synchrony with Universal Time and not with local solar noon, none of the four daily profiles ( Fig. 3b-i, paired panels, top: dotted traces) could be expected to occur exactly at solar noon (paired panels, top: topmost grey trace). Typically though, with every ITP deployed at least one of these four daily profiles occurred within 2 h of solar noon and provided an acceptable proxy for "daily" PAR profiles, for visualization.
These near-noon penetration depths ( Fig. 3b-i, paired panels, bottom) exhibited seasonal and regional variability broadly consistent with the expected seasonal and latitudinal forcing in solar insolation (Fig. 3, top paired panels: grey traces). Overall, PAR penetrated deeper into the water column at lower latitudes and during months of the year when solar elevations were highest. The depth horizons of minimum measurable PAR generally shoaled between September and October in these ITP time series, presumably reflecting decreases in solar elevation and thus in insolation as both fell to zero. The most northerly deployed ITP (ITP48 at  3b). By the beginning of October, PAR was wholly unmeasurable by ITP48 at its topmost sampled depth of $7 m. Two ITPs that were deployed somewhat further south (ITP72 and ITP65, both around 818N) were able to measure PAR down to $40 m in the water column in early September (Fig. 3c,e).
Within 2 weeks, the latter was unable to measure PAR above the threshold level of detection while the former remained able to detect measurable PAR until early October. Two ITPs that were deployed even further south (ITP52 at 788N and ITP69 at 758N) were able to measure PAR throughout September to depths well below 100 m in the Beaufort Gyre (Fig. 3d,f) but by late October neither system was able to Remaining panels (b-i) are paired, showing maximum and minimum daily local solar elevation (grey traces, top panel) and actual solar elevation at every profile (up to four times daily, black dots, top panel). Bottom panels in each pair show the daily maximum depths over which these above-zero PAR data were observed.
Laney et al. The euphotic zone under Arctic sea ice measure PAR in the water column at sampling depths any deeper than their $7 m upper stop. The remaining system (ITP64, released in August into open water in broken pack ice at roughly 798N) was able to measure PAR to only $50 m depth through the month of September and became unable to detect water column PAR at all by mid-October (Fig. 3h).
In several of these time series, profiles were observed which did not fully complete their expected vertical excursions (e.g., mid-September for ITP52 or intermittently throughout the summer for ITP64). These likely represent situations where the profiling vehicle's traction drive experienced slippage on the tether while profiling, as can occur when the sea ice is drifting strongly relative to the underlying water column.
Only three of the initial eight ITP systems returned PAR profiles in the springtime following their deployment the previous fall (Fig. 3a). The earliest springtime profile was transmitted by ITP69 in late February which indicated that PAR was measurable to depths of $30 m this early in the season (Fig. 3g). The single profile returned from ITP65 in mid-April of 2013 indicated PAR measurable to depths of only $11 m (profile not shown). The third system (ITP64) transmitted valid PAR profiles beginning in late April 2013, measurable in some cases to depths exceeding 50 m (Fig. 3i). Profiles with nonzero PAR intensities were observed intermittently between April and late May, and given the short time scales and large magnitudes of change in under-ice PAR seen during this period it is likely that these rapid changes in under-ice PAR reflect gross changes to the surface ice layer, such as the opening and later refreezing of leads in the overlying sea ice that is known to occur at this time of year (Assmy et al. 2017). Changes in the abundance of ice algae, which would also affect transmittance over this early part of the growth season, would presumably introduce less dramatic modulations in under-ice PAR that moreover would occur with somewhat longer time scales. Following this period of highly intermittent above-zero PAR data, measurements of PAR became considerably more consistent and penetrated progressively deeper into the upper ocean to remain measurable down to $70 m for the rest of the summer.
This variability in the penetration depths of PAR was accompanied by strong differences in its characteristic vertical distribution. PAR profiles did not always display the canonical exponential decrease with depth that would be expected in ice-free oceans but instead often exhibited anomalously low intensities within the upper 30 m (Fig. 4). These anomalous vertical features are most evident in late fall profiles from ITP52 and ITP69 in the Beaufort Gyre (Fig.  4a,c), with the former exhibiting broader apparent subsurface "maxima" compared to latter whose profiles appeared more sharply diminished at the topmost-sampled depths. In contrast, late-fall PAR profiles from ITP65 in the Beaufort Gyre exhibited only slight deviations from exponential in its topmost-sampled depths (Fig. 4b). With ITP64, which survived its open-water deployment in the Beaufort Gyre and the subsequent polar winter to begin consistently detecting PAR under ice by late May, its characteristic vertical distributions of PAR changed over the Arctic summer. Maximum PAR values occurred initially at $30 m and subsequently shoaled, becoming more pronounced by June and mid-August (Fig. 4e). Profiles returned by the most northerly ITPs deployed in the Transpolar Drift (ITP48 and ITP72) did not exhibit such strong near-surface anomalies in PAR and instead generally displayed a seemingly exponential attenuation of PAR with depth, similar to what would be expected in ice-free water columns (Fig. 4f,g). However, slight deviations from purely exponential can be seen in the topmost depths of some of those profiles.
This data set included PAR profiles from Arctic water columns that experienced only partial ice cover, i.e., in the ITP64 time series during its initial release into open water among broken pack ice. This near-surface anomaly in PAR was sometimes seen in the topmost-sampled meter, albeit intermittently, during this period in broken pack ice (Fig. 4d). It is important to note that for the remaining ITP64 open water profiles where no near-surface anomaly was visually apparent, as with all other ITP profile data which similarly appeared "exponential" in its vertical distribution, the absence of such a near-surface anomaly does not necessary indicate that none was present. If any deviation in PAR occurred higher in the water column than the ITPs' upper stop at $7 m, this would not be evident in a PAR profile.

Parameterizing PAR profiles made under sea ice
These near-surface anomalies in irradiance intensity ( Fig.  4) appear similar to what Frey et al. (2011) observed in profiles taken manually in the Chukchi Sea, through sea ice that was characterized by a substantial areal fraction of melt ponds. The canonical exponential model for parameterizing light profiles in ice-free water columns (e.g., Kirk 1994;Mobley 1994) fitted their observed under-ice irradiance profiles poorly, and so Frey et al. (2011) developed an alternate parameterization based on a conceptual model of irradiance profiles performed through a region of bare ice surrounded by melt ponds that occupy an areal fraction P of the ice cover: Here, I is the irradiance intensity transmitted through the ice layer and detectable immediately under the center bare ice region at the ice-ocean interface. P is defined as above, and an "enhancement factor" N quantifies the additional light transmittance through the melt ponds relative to transmittance through bare ice. The water column diffuse attenuation coefficient k is defined identically as in Eq. 1, and a Laney et al. The euphotic zone under Arctic sea ice fifth parameter h i arctan R=z ð Þ quantifies the depthdependent influence of shading by the center bare ice region of radius R.
Our initial approach was to parameterize our PAR profile time series using this Frey et al. (2011) model, by fitting Eq. 2 to our six time series of I PAR (z) data. Because the number of samples in each ITP-measured profile greatly exceeds the number of parameters in Eq. 2, the model would be overconstrained and in principle we could recover robust estimates of all five parameters of Eq. 2 by fitting it to observed profiles of I PAR (z). We observed that Eq. 2 could sometimes provide qualitatively acceptable fits for PAR profiles, given careful attention to setting proper parameter bounds (e.g., Fig. 5a), but we were not able to obtain consistently good fits using Eq. 2. Moreover, we identified several issues with this model that we believed made it unsuitable for further use for parameterizing these ITP-obtained PAR profiles. The first concern was mathematical in that Eq. 2 reduces to Eq. 1 both when N 5 1 or when P 5 0. This presents a challenge when using Eq. 2 to interpret field observations when N and P are neither known a priori nor measured directly, as is the case with our data set, even when the model is overconstrained by the data available for fitting. This mathematical issue becomes especially problematic when measured I(z) Fig. 4. Stagger plots of PAR profiles measured by the six ITP systems deployed in the Beaufort Gyre (a-e) and in the central Arctic (f,g). Panels for ITP64 (d,e) are split to omit the winter months when no insolation occurred. For clarity, these panels show only one "midday" profile every 2 d, not the entire complement of four daily profiles. These plots are intended to illustrate the characteristic vertical distribution of PAR and are not scaled identically in magnitude across all ITPs. Exact intensities for each PAR profile time series are represented by color, with Roman numerals in parenthesis indicating the corresponding color bar scale.
profiles only deviate slightly from exponential, because the source of this deviation cannot be attributed specifically to either N or to P using Eq. 2. This joint influence of N and P became strongly evident during our initial attempts to fit our profile time series using Eq. 2, where best-fit values of N and P were found to be very tightly correlated which In panels (c-i), fits of profiles that were best fit by the exponential model of Eq. 1 are shown with dashed grey lines; profiles that were best fit with Eq. 3 are indicated with solid grey lines.
indicated poor parameter independence and in turn, an inappropriate model to use when certain model parameters (here, N and P) are not known a priori.
A second concern with Eq. 2 was that PAR profiles returned by these ITPs do not include data from the topmost 7 m which is where several of its variables, particularly I, have their greatest effect. Therefore, even though it is possible in principle to solve this N-P independence problem by deriving a comparable equation which expresses the joint effect of N and P in a single new parameter, this second issue of not having observations above 7 m would still hinder the use of this hypothetical new model by not having sufficient data to constrain estimates of I and to a lesser extent the joint effect of N and P.
Our third concern with using the Frey et al. (2011) model to parameterize our data set was that doing so would implicitly assign a specific mechanism as being responsible for the vertical distributions we observed in our time series, when in fact we could not be assured that this is the case. Our six ITP time series were collected over seasonal time scales in broadly different regions of the Arctic Ocean, and we could not assume a priori that melt ponds per se were the sole or even most frequent mechanism causing these observed nearsurface anomalies in PAR. Other processes in the sea ice cover can also introduce unusual anomalies in vertical profiles of PAR (e.g., Katlein et al. 2016) and would presumably require a different model based on other mechanistic processes than melt ponding alone to explain the vertical distributions seen our entire ensemble of I PAR (z) profiles.
Given these concerns with adopting the Frey et al. (2011) melt pond model we instead developed a semi-empirical parameterization to explain these anomalies independent of any specific physical mechanism. Similarly to both Eqs. 1 and 2, this semi-empirical parameterization assumes a water column with vertically uniform PAR attenuation properties that are expressed by a single diffuse attenuation coefficient k D , defined identically to that presented in Eq. 1: I 0 can be roughly interpreted as the approximate irradiance that would be observed at the ice-ocean interface if the overlying ice cover were spatially homogeneous. The I 0 in Eq. 3 differs slightly from the nominal "surface" irradiance I 0 of Eq. 1 in that the former is by necessity referenced to the zero isobar that, strictly speaking, is located within the surface ice layer (sea ice floats isostatically on the ocean surface with freeboard above and keels below the zero isobar). I 0 could instead be referenced to the ice-ocean interface by introducing a vertical offset similar to the approach suggested by Katlein et al. (2016), but as ice thickness data were largely lacking throughout these six time series such offsetting was ignored for simplicity in our subsequent analyses. We introduced a second exponential coefficient k NS to describe the near-surface decrease of this predicted I 0 to its actual intensity at the ice-ocean interface (I NS ), due to whichever mechanisms in the ice are causing this nearsurface diminishment of PAR, unspecified in this model. A notable aspect of this model is that k NS could in principle be negative, representing a situation where for example light transmittance is greater in the immediate area around the vertical axis of a PAR profile, and lesser in areas further distant. This would create a vertical distribution in PAR that is not anomalously diminished in the upper water column but rather amplified, as might be predicted when a ship clears heavy ice cover from around itself and then performs a hydrocast within a cleared area to measure a PAR profile. However, to explain the characteristic shapes seen in these autonomously sampled ITP profiles k NS would always be positive, reducing I 0 in the topmost depths beyond what would be attributed to water column attenuation alone. This new semi-empirical parameterization provided seemingly acceptable fits to PAR profiles throughout these six time series (e.g., Fig. 5b), but not all PAR profiles in these time series exhibited noticeable near-surface anomalies in PAR and many did not visually appear to be different from a simple exponential (e.g., examples seen in Fig. 4). It is possible to force a fit of Eq. 3 to all PAR profiles in these time series using a bounded curve fitting approach (Coleman and Li 1994;Coleman and Li 1996), but we instead chose to use a model fitting and selection approach to identify those profiles for which Eq. 1 might instead provide a better fit in a statistical sense. This would in principle distinguish times and locations where putative heterogeneity in surface shading by sea ice was significant in these time series, from times and locations where it was not. To accomplish this, we examined every measured profile of I PAR (z) to determine if it was better fit by the new candidate model of Eq. 3 or by the simpler exponential model of Eq. 1 instead. We did not include the parameterization expressed by the Frey et al. (2011) "melt pond" model (Eq. 2) in this comparison for the reasons described above. The lsqcurvefit routine in Matlab (R2014b, The Mathworks, Natick, Massachusetts) was used to Table 3. Parameters for the semi-empirical model for underice irradiance profiles (Eq. 3), showing units and bounds used when fitting this model to these six measured I PAR (z) profile time series. The maximum bound on I NS for any given profile was set at the maximum PAR intensity measured during that individual profile.

Model parameter Units
Bounds during fitting exercise Laney et al. The euphotic zone under Arctic sea ice compute the best-fit parameters for Eqs. 1, 3 for every PAR profile, chosen because it allows boundary conditions to be placed on any or all model parameters (Table 3). The adjusted coefficient of determination ( R 2 ) was then computed for each pair of fits and the model with the higher R 2 was deemed to be a better representation of that particular I PAR (z) profile in a statistical sense. The adjusted coefficient accounts for any spurious advantage that Eq. 3 would have over Eq. 1 by virtue of having more explanatory variables (Zar 1999). For the purposes of this model fitting and selection exercise, profiles that were flagged as described in Methods (those without data above the 15 m horizon, or otherwise with a vertical span less than 10 m) were excluded in order to avoid potentially spurious parameter estimates or inferences that might arise from retaining inadequately characteristic profiles in this analysis. By using 15 m as the upper horizon, we are able to retain profiles from ITP72, roughly one third of our usable PAR profiles from the Central Arctic. This model selection exercise revealed systematic differences in whether Eq. 3 or Eq. 1 provided a better fit to any given profile in these six time series (Fig. 5c-i), in turn indicating where and when the under-ice light field exhibited statistically meaningful near-surface anomalies in I PAR (z). For profiles collected in the central Arctic Ocean, the exponential model of Eq. 1 provided a statistically better characterization of I PAR (z) for 21% and 100% of profiles from ITP48 and ITP72, respectively (Table 4). In contrast, PAR profiles collected by the three Beaufort Gyre ITPs that were deployed initially in sea ice (ITP52, ITP65, and ITP69) were almost uniformly better fit by the semi-empirical model of Eq. 3 for 97%, 100%, and 100% of all profiles, respectively. For the initial period of ITP64's deployment after its release in the Beaufort Gyre in open water among broken pack ice, Eq. 3 provided a statistically better fit for 54% of its measured profiles. In 2013 the following spring, after insolation resumed, 94% of ITP64's PAR profiles were better fit by Eq. 3, indicating consistent near-surface anomalies.
The choice of using only profiles with data above a 15 m minimum upper depth horizon and having at least 10 m of continuous vertical data is necessarily subjective. To examine the sensitivity of this fitting exercise to these restrictions on profile length and depth, we performed an identical analysis in which we set the upper horizon 5 m higher in the water column, and required profiles to have at least 15 m of continuous vertical data. These more conservative restrictions both lengthen the profiles used to assess the fit of Eq. 3 vs. Eq. 1 and shift these profiles to higher in the water column, where the near-surface shading effect would be more pronounced. This new set of profiles however excluded all from both ITP48 and ITP72, our two Central Arctic profilers whose PAR profiles did not penetrate deep enough into the water column meet these criteria. Results using this set of profiles indicated that Eq. 3 best fit this higher-quality data set in even more cases than originally: now 100% for ITP65 and ITP69 as before but increasing ITP52 to 100% and ITP64 (2013) to 99%. Only ITP64's 2012 period remained relatively unchanged, now fit by Eq. 3 in 55% instead of 54% of cases with the prior shorter and deeper set of profiles.
Even though we do not believe the Frey et al. (2011) model is an appropriate parameterization to apply to these ITPenabled PAR profiles for the reasons detailed above, for completeness we examined how well Eq. 2 fit these PAR profile data in a statistical sense using our original restrictions on profile upper depth (profiles with samples above 15 m) and length (profiles with minimum 10 m vertical extent). In a three-way comparison between the three models of Eqs. 1-3, we observed no clear preference for Eq. 2 over Eq. 3 among these six time series (Table 4). This suggested that even if we were to discount the above arguments regarding the appropriateness of Eq. 2, roughly half of the observed I PAR (z) profiles are still not fit best by the melt pond-based parameterization of Frey et al. (2011) and are better represented by the semi-empirical model. Table 4. For these six PAR profile time series, the number of profiles considered appropriate for the curve fitting and selection exercise, and the number and percent fit best by Eq. 1 (purely exponential model) and Eq. 3 (semi-empirical parameterization), as well as by a three-way comparison between Eqs. 1-3 (final three columns). Parameter estimates: magnitudes and trends This model fitting and selection process helped to quantify the seasonal and regional differences in the vertical character of I PAR (z) that are apparent visually in Fig. 4. The parameters derived from this model fitting exercise further quantify these differences in terms of specific aspects of these under-ice light fields such as the diffuse attenuation coefficient k D . During the months of August to November when these six time series overlapped the most, the greatest difference in estimates of k D was seen between the central Arctic Ocean and the Beaufort Gyre (Fig. 6a). Estimates of k D from the central Arctic were higher, from $0.12 m 21 to $0.17 m 21 for ITP48 and from $0.11 m 21 to $0.13 m 21 for ITP72, compared to the much lower values of $0.03 m 21 and $0.07 m 21 that characterized water columns in the Beaufort Gyre. These observed trends in central Arctic  Table 3). (e) Mean-square error for the fits generating these best-fit parameter estimates. (f) Depth of the 1% light level z eu . Closed symbols indicate profiles better fit by Eq. 1 than Eq. 3. estimates of k D were not strongly biased by the model used to determine them (compare filled symbols for Eq. 1 with open symbols for Eq. 3). The lowest estimates of k D in these time series were within ITP64's time series and tended to occur early in the growth season, between late April and July ($0.02 m 21 ). Generally, over the subsequent summer period estimates of k D derived from the ITP64 time series rose to $0.05 m 21 and remained steady.
The estimates of I 0 that were derived from these six time series encompassed almost four orders of magnitude (Fig.  6b). Qualitatively, trends generally followed the expected seasonal forcing by insolation and sea ice cover. The sole data available from early in the year (from ITP64 in the Beaufort Gyre) indicate that when PAR is present very early in the year in late April, estimates of I 0 can be as high as were seen much later in the summer. This is consistent with the rapid opening and closing of leads as suggested earlier.
Estimates of I 0 extending into May arise primarily from nonexponential I PAR (z) distributions (open triangles), persisting well after the daily light cycles had transitioned into the full 24 h summer insolation. Some April profiles in contrast were better fit by Eq. 1, indicative of intermittent, more exponential-like distributions of PAR (closed triangles). By late May I 0 rose rapidly, likely reflecting abrupt increases in sea ice transmittance associated with, for example, the loss of surface snow and/or thinning of the sea ice itself. Moreover, loss of snow and thinning of ice also often coincides with the decline in algal communities resident in the sea ice layer (Cota and Smith 1991), and their release into the water column may further increase light transmittance by the surface ice layer. PAR profiles from ITP64 were better fit by the semi-empirical model of Eq. 3 (open triangles) for the majority of the summer and fall. Decreasing trends in I 0 seen during August to November in the remaining five ITPs presumably mirror the seasonal decrease in solar elevation and thus insolation.
The two empirically defined parameters of Eq. 3 (I NS and k NS ) were obtained only from those PAR profiles where the semi-empirical model of Eq. 1 was statistically more appropriate (i.e., open symbols). Thus, time series of these parameters contain fewer points than were generated for k D and I 0 .
Estimates of I NS quantify the approximate PAR irradiance at the top of the water column, and best-fits of I NS encompassed three orders of magnitude albeit often converging at the lower bound that was set for this parameter (0.01 lE m 22 s 21 ; Fig. 6c). Ignoring the values at this lower bound, the remaining estimates of I NS exhibited a seasonal trend that roughly mirrored that in I 0 , rising from very low levels in early June and later decreasing over the September to October time frame. A notable difference in the time series of I 0 and I NS was apparent in early summer, where the initial increase in I NS lagged that of I 0 and exhibited a more gradual rate of change than was seen earlier in I 0 at the end of May.
The fourth parameter of Eq. 3 (k NS ) in a sense quantifies the sharpness of the near-surface shading in PAR. As defined in Eq. 3, this parameter could in principle be either positive or negative, and for completeness Eq. 3 was first fit to all six PAR profile time series with no restriction on the sign of k NS to identify any PAR profiles with negative best-fit k NS that might indicate greater light transmittance by sea ice in the area immediately local to the ITP installation. In no instance was the best fit value of k NS negative when allowed to be and therefore k NS was assumed to always act to reduce nearsurface irradiance in these six time series of ITP-obtained PAR profiles. Estimates of k NS spanned one and a half orders of magnitude and best fit values always remained above the lower-set parameter bound of zero (Fig. 6d). For the springsummer months in the ITP64 time series, estimated magnitudes of k NS were generally very low, indicating relatively weak near-surface shadings in PAR. Estimates of k NS increased toward the end of summer and into fall, indicating sharper near-surface shadings and quantifying the changes over time that are visually apparent in the vertical distributions of PAR observed by ITP64 (Fig. 4e). For the remaining three Beaufort Gyre ITPs, estimates of k NS generally increased over the late-summer and fall period, suggesting progressively stronger near-surface anomalies in under-ice PAR. For the two ITPs in the central Arctic, the subset of PAR profiles that exhibited near-surface anomalies coincided with best fit values of k NS that were generally higher than those observed in the Beaufort Gyre, indicating sharper near-surface anomalies in PAR. However, comparable k NS magnitudes during this time of the year can also be seen in the Beaufort Gyre at the very end of the ITP64 time series, before transmissions from this profiler ceased.
The mean squared error (MSE) of fitting Eq. 3 to these daily PAR profiles was also examined to see if these observed trends in estimates of k NS or I NS were potentially affected by biases in the quality of fit (Fig. 6d). With ITP64, observed MSE stayed at a relatively constant level for the summer, indicating relatively good quality of fit. Notably, the poorest fits for ITP64 PAR profiles were associated with those best fit by the simple exponential model of Eq. 1 (filled triangles), not those best fit by the semi-empirical model of Eq. 3 (open triangles). In the remaining five PAR profile time series MSEs, all decreased considerably and progressively between early August and November, indicating increasingly better fits to the semi-empirical model of Eq. 3 over this period. Visually, within the range of MSEs deemed acceptable for this analysis, even fits with the "poorest" quality (e.g., examples shown in Fig. 5) capture the first-order character of these complex under-ice vertical distributions of PAR.

Euphotic zone trends and bio-optical influences
For ice-covered water columns that exhibit such grosslynon-exponential vertical distributions of PAR like those seen here in the Arctic Ocean, computing euphotic zone depths , are needed to isolate the influence of diffuse attenuation by the water column itself from that of spatial heterogeneity in shading by the overlying sea ice. Once a water column attenuation coefficient can be estimated, Eq. 1 can then be used to compute the depth of the euphotic zone in terms of any given percent light level. For this study, we used the 1% light level to demarcate the bottom of the euphotic zone. The horizons of z eu(1%) we computed were on average shallower in the central Arctic compared to the Beaufort Gyre: $40-50 m vs. $70-110 m, respectively (Fig. 6f). The euphotic zone appeared slightly deeper on average early in the year in ITP64's time series, and progressively deeper excursions of z eu(1%) were seen the three other Beaufort Gyre ITPs starting in the late September. This phenomenon of z eu(1%) being deeper in the spring and fall is consistent with light penetration being stronger in these shoulder seasons when the attenuation of this light by phytoplankton in the water column is presumably less, due to lower phytoplankton biomass during these periods. Such influences can in principle be explored with these data sets by examining the profiles of chlorophyll, colored dissolved organic matter, and particulate backscattering that were concurrently measured by these ITPs. These three biooptical properties should to some degree influence water column light attenuation and so we examined their potential control on k D using a statistical approach. We first computed the depth-integrated values of these bio-optical properties over the top 100 m (i.e., from $7 m at the top of each ITP profile to 100 m depth) to serve as explanatory variables in a multiple linear regression analysis to examine how changes in water column chlorophyll concentration, optical backscattering, and colored dissolved organic matter (CDOM) may affect k D . This depth range parallels the vertical domain over which our model selection and fitting exercises were conducted (Fig. 5) and encompasses the depths of the vertical excursion in the chlorophyll maximum that occurs in the Beaufort Gyre over the course of a year ). Ideally, we could then use these depth-integrated time series of chlorophyll, CDOM, and scattering as explanatory variables in a basic multiple linear regression analysis where our model-estimated k D was the response variable, for each of the six time series, However, for several of these time series there were too few profiles having sufficient quality to obtain estimates of k D , to expect much confidence in results from a standard multiple regression analyses even despite the seemingly favorable p values sometimes computed in these regressions (Table 5) Only the two time series with the largest number of k D estimates (ITP52 and the latter 2013 period of ITP64's deployment) resulted in regressions significant at the p < 0.01 level, although ITP64's initial 2012 time series was only marginally above this threshold. These three nominally significant regressions all represented time series from the Beaufort Gyre (ITP52 and ITP64), but beyond that there was no uniform pattern among these three in the particular bio-optical variable(s) that had regression coefficients that were individually significant.
To explore in more detail the relative influence of these three bio-optical variables, we performed an additional stepwise regression (Johnson 1998) using the Matlab stepwisefit function (Matlab Statistics toolbox v9.1, The Mathworks, Natick, Massachusetts). This stepwise approach again identified the same three time series (ITP52 and both subsets of ITP64) as being significant at the p < 0.05 level (Table 5). Table 5. Results of a multiple linear regression analyses where the estimates of k D recovered from the model fitting and selection analysis (Table 4) are regressed against concurrently collected data of chlorophyll (chl), optical backscattering (b b ), and colored dissolved organic matter (CDOM) concentrations, used as explanatory variables. Italics indicate overall and coefficient-specific p values less than 0.01 in both the basic and stepwise multiple linear regressions.

Multiple linear regression
Stepwise regression Again however we observed no uniform pattern among these time series as to whether variability in chlorophyll, CDOM, or particulate backscatter were most correlated with estimated changes in k D that we derived during our model fitting and selection exercise.

Assessment of potential artifacts
Because this study relied heavily on autonomouslysampled observations, it was critical to consider ways in which potential sources of artifact could affect these PAR profiles and the inferences we drew from them. We attributed the observed near-surface anomalies in PAR to a known shading phenomenon that can occur when surface sea ice immediately above the axis of an irradiance profile is optically less transmittant than in surrounding areas (Frey et al. 2011). However, the surface buoys of these ITPs will also contribute some degree of analogous shading, because they sit on the sea ice directly above the profiling vehicle. One key assumption we make is that any shading effect arising from the surface buoys themselves contributes negligibly to the near-surface phenomena we observed in PAR. We believe this assumption is reasonably supported by several lines of argument. First, if buoy shading had any substantial influence on the under-ice light field, we should expect it to also create a strong shading effect in the absence of sea ice, such as during the initial open-water deployment period of ITP64. Over this period, ITP64's surface buoy was presumably floating free in the ocean before winter freeze-up, yet at no time during this initial open-water period were strong nearsurface anomalies in PAR evident (Fig. 4d). Only small anomalies in PAR were seen at the very tops of some profiles and moreover not consistently in every profile as might be expected for any direct shading artifact that might arise from the buoy itself that is always present.
We can extend this analysis and estimate the expected magnitude of the buoy's direct shading effect using a geometric approach accounting for the dimensions of the profiler, tether, and buoy and allowing for the cosine characteristics of the PAR sensor (Fig. 7, left). Two types of surface buoys were used for these six ITPs: an inverted truncated conical buoy (left: "A") for all systems other than for ITP52 which used an earlier, shorter cylindrical design ("B"). The base diameter of both buoy designs is comparable ($0.66 m) and so both buoy types can be considered roughly equivalent in the amount of direct shading each introduces. The profiler experiences the maximal amount of direct buoy shading at the top of every profile ($7 m), and when floating free among broken pack ice as in ITP64's initial deployment period, this buoy shading appears to the profiler as an opaque disc (Fig. 7, left: "direct shading cone," short dashed lines). The percent of the light field sensed by the PAR sensor at 7 m can be computed using the following relationship: where we here integrate the downwelling irradiance only from 08 azimuth to 458 instead of the full downwelling hemisphere, to provide a conservative over-estimate of any buoy-driven shading artifact. Given these stipulations, the profiler at $7 m depth will experience only $0.5% of shading of incident PAR due to the surface buoy sitting directly overhead (Fig. 7a, solid line). This is considerably less shading than is typically observed in near-surface PAR when it deviates from an expected exponential distribution (e.g., Fig.  4) and thus buoy-driven shading would appear to contribute negligibly to the near-surface anomalies observed in our I PAR (z) time series. This minimal contribution of direct buoy self-shading is too small even to account for the minimal near-surface anomalies seen intermittently during ITP64's open-water phase among broken pack ice. One hypothesis is that those intermittent small near-surface anomalies observed during ITP64's open water period do in fact reflect a degree of shading by surface sea ice, in situations when the buoy was driven by the wind against, or became trapped between, ice floes in the area.
To compare more quantitatively this estimated degree of direct buoy shading with the putative shading seen in our time series we used an approximative approach, independent of our semi-empirical model, to estimate what I PAR (z) might look like in the absence of local surface shading. We created families of estimated "no-shading" PAR profiles by extrapolating I PAR (z) profile data from depths below a certain horizon, starting at 20 m and progressing down to 50 m in 10 m intervals, back up to the surface (Fig. 7b, dashed lines). The percent difference can be then computed between these estimated no-shading PAR profiles and profiles where we did observe near-surface anomalies (Fig. 7c). This approach estimated the amount of typical shading from all sources to be typically 50% or higher at the 7 m depth horizon, which is two orders of magnitude greater than the estimated $0.5% direct shading at this depth that the buoy would contribute. This analysis more quantitatively suggests that direct buoy shading is unlikely to contribute meaningfully to the nearsurface shading anomalies we see in these PAR profiles.
In situations when sea ice is present, the buoy no longer appears as a sharply defined opaque disc from the perspective of the profiler but rather appears as a larger diffuse disc because the high degree of scattering by the sea ice introduces an enlarging effect on the apparent shadow (Fig. 7, left: "scattered shading cone," long dashed lines). However, and most importantly, to the first order this enlarging effect introduces no net change in the amount of light shaded by a buoy. The buoy's shaded area appears larger to the profiler but also brighter in intensity than what would observe in the direct shading cone in the absence of sea ice. This phenomenon is very difficult to model robustly but the areal increase and magnitude of brightening in the buoy's shadow can be roughly estimated using radiative transfer models. A recent Monte Carlo simulation study by Petrich et al. (2012) examined horizontal distributions of transmitted light under an ideal, semi-infinite obstruction sitting on sea ice of 1 m thickness, which is analogous to the scenario with these ITPs. These Monte Carlo simulations indicate that the shadow cast by a semi-infinite obstruction on 1 m sea ice will be $90% diminished at a distance 1 m from the edge of the obstruction. Inversely, the bottom of the ice layer 1 m horizontally underneath the obstruction will appear 10% as bright as the far field irradiance well away from the obstruction. Our ITP buoys are finite obstructions (i.e., cylindrical, not infinitely deep walls) and these above modeling results therefore likely overestimate both the diameter and the magnitude of any shadowing expansion that actually occurs with these ITP buoys in 1 m thick sea ice. Regardless, when these simulation results are used to compute the percent of PAR shaded by the buoy when sitting on sea ice $1 m thick, we can approximate the diffusely shaded region as having a diameter of 2.66 m, i.e., the buoy physical diameter plus a 1 m diffuse shadow radially around it (depicted in Fig. 7, left; greyed region under buoy "A"). If we overestimate and assume that this shadow is equally as dark as the smaller direct cone (i.e., ignore the brightening effect), then at the 7 m safety stop the profiler would experience only $8% shading of incident PAR (Fig. 7a dashed line). This 8%, which represents a substantial overestimation, is still too little to explain the typically 50% shading effect seen in many PAR profiles near the surface. In the more realistic scenario where the obstruction is finite with only a $0.66 m diameter and when the brightening effect is not ignored, the actual direct shading by the buoy can be expected to be much less. Therefore, these in-ice estimates of the buoy's direct shading in-ice (dashed) conditions. (b) Example of initial PAR profile collected by ITP52 (dark trace) and corresponding chlorophyll concentration (grey symbols, plotted at x10 concentration). Additional dashed lines indicate purely exponential PAR extinction curves extrapolated to the surface using PAR data at depths below 20 m, 30 m, 40 m, and 50 m as noted. (c) The percent difference between the extrapolated exponential PAR profiles in (b) and the actual measured PAR for this profile.
effect are still too small to be considered a meaningful contribution to the putative near-surface shading we observe throughout these PAR profile time series.
The above arguments address artifacts arising from direct shading by the buoy itself, but the presence of the buoy on the sea ice may also introduce indirect shading artifacts on these under-ice PAR profiles. For example, when surface snow is seasonally present the buoys may act to dam drifting snow which would make the area above the profiler considerably less transmittant than surrounding areas nearby. Similarly, during ice melt in spring and summer the buoy itself may insulate the ice immediately below it from surface melting, so that adjacent ice might become relatively thinner and relatively more transmittant. Both of these phenomena would affect vertical profiles of irradiance in a manner similar to how the Frey et al. (2011) melt pond mechanism introduces nearsurface anomalies in under-ice irradiance profiles. Alternatively, the surface buoys might enhance scouring of snow around themselves in windy conditions and therefore increase light transmittance immediately around themselves, introducing near-surface enhancement in profiles of PAR. These various sources of indirect artifact would not contribute to any putative near-surface shading in I PAR (z) seen initially at the beginning of any given time series, immediately after these ITPs are deployed, but they could in principle contribute to the longerterm variability we observed and in particular our modelderived estimates of the sea ice parameters I NS and k NS (Fig. 6).
The potential contribution of such indirect buoy shading artifacts on I PAR (z) is very difficult to predict given the potential interactions between various such artifacts and the absence of companion data regarding the local sea ice and snow cover throughout these time series. Moreover, the inherent complexity of how phenomena such as snow drifts or scouring would affect the transmittance of incident sunlight around the buoy is challenging to model. We can estimate the spatial scales of these indirect phenomena to be one to several meters around the buoy itself, which is considerably less than the characteristic areal scales of $1000 m 2 or so that radiative transfer models suggest are necessary to generate the types of near-surface anomalies we observed in our PAR profiles (Katlein et al. 2015), at depths exceeding 30-40 m. For this reason, we also consider that indirect buoy artifacts have a similarly negligible shading contribution on these I PAR (z) profiles, although we recognize that such artifacts cannot be ruled out absolutely. When using autonomous systems the potential for artifacts to arise are numerous, especially when examining an environmental property as highly sensitive to local conditions as irradiance.

Irradiance distributions under Arctic sea ice
This study examined six long-term observational time series of the vertical distribution of PAR in ice-covered regions of the Arctic Ocean. The transmittance of light through Arctic sea ice is well studied both in general (e.g., Grenfell and Maykut 1977;Perovich et al. 1998;Light et al. 2008;Nicolaus et al. 2010a) and in terms of using its spectral characteristics to assess resident ice algal communities (e.g., Maykut and Grenfell 1975;Mundy et al. 2007;Campbell et al. 2014), but in comparison few studies have examined the light that is transmitted through Arctic sea ice for its own importance in establishing the euphotic zone in the underlying water column. This study generates important information on the structure and variability of the euphotic zone under sea ice in the Arctic Ocean, which in turn provides new insight into how these may modulate photosynthesis, production, and phytoplankton ecology in Arctic pelagic ecosystems regionally over seasonal scales. To the best of our knowledge, this work represents the first-ever observational study to examine light penetration into underice water columns in the Arctic on sub-daily to seasonal scales, as well as the first detailed examination of its vertical distribution under sea ice in different regions of the Arctic Ocean.
Given the substantial challenges with measuring basic ecological processes in ice-covered water columns, biogeochemical models are especially valuable for assessing climate-driven responses in Arctic marine ecosystems (e.g., Zhang et al. 2010;Jin et al. 2012;Long et al. 2015;Zhang et al. 2015). Uncertainties in these models can be reduced by having better observational constraints on key ecosystem variables such as under-ice irradiance intensities and PAR penetration depth. Polar-suitable autonomous profiling systems such as the ITP and others (e.g., Kikuchi et al. 2007;Prisenberg et al. 2007) can enable the long-term collection of such core ecosystem properties. Yet autonomous approaches remain largely underutilized for long-term ocean ecosystem assessment in the Arctic despite a well-recognized need for better observational capabilities in this rapidly changing region (National Research Council 2006;National Science Foundation 2007). This study serves as an important demonstration that such core ecosystem properties as euphotic zone depth, and other properties as well such as phytoplankton biomass (through chlorophyll or scattering), can be effectively measured in ice-covered regions of the Arctic Ocean autonomously and over the long temporal scales needed to better constrain regional biogeochemical models involving Arctic under-ice light fields.
Beyond the collection of these novel time series themselves, this study provides valuable insight into how polarspecific phenomena can influence the distribution of light in marine water columns. The putative sea ice shading effect seen here is a unique feature of under-ice light fields which can be explained by the overlying sea ice layer being horizontally heterogeneous in its optical transmittance, admitting less light into the water column immediately above the vertical profiling axis than in other areas nearby (Frey et al. 2011;Katlein et al. 2016). A key finding of our study was the frequency and magnitude with which such distinctive nearsurface irradiance shading effects were present in these icecovered water columns (e.g., Fig. 5). This shading effect was observed in the Beaufort Gyre persistently throughout the entire insolated season but much less frequently in the central Arctic. As the optical heterogeneity needed in the surface ice layer to cause such near-surface anomalies in I PAR (z) likely has areal scales horizontally on the order of $1000 m 2 or so (Katlein et al. 2015), these observations should not be considered representative of all ice-covered euphotic zones everywhere in the Arctic. Profiles made in a region with sea ice cover that is more optically homogeneous, characterized by considerable snow cover or thicker ice, for example, would presumably not display such near-surface shadings in PAR or might display them more weakly. Such regions might instead be expected to exhibit vertical distributions of PAR under sea ice that are more closely exponential, due to higher spatially uniformity in the light transmittance characteristics of sea ice on these local scales (Katlein et al. 2016). This potentially explains the observations from ITP48 and ITP72 in the central Arctic where an absence of apparent near-surface anomalies in PAR was noted in times of the year when factors that might introduce substantial optical heterogeneity in the surface ice layer, such as melt ponds, would not be expected.
The primary focus of this study was to examine the intensity and vertical distribution of PAR in ice-covered water columns over seasonal to annual scales. The prevalence of the observed near-surface shading phenomenon has important ramifications for how under-ice PAR should be measured and modeled, but this phenomenon will also presumably affect water column photosynthesis and primary production to some unknown extent. Estimating this influence is not trivial. These near-surface shading features represent the integrated result of sea ice transmittance characteristics within a larger region around the ITP surface buoy, and to create these PAR profile shapes this region presumably contains within it smaller areas that are relatively more and less transmittant. Phytoplankton in the water column immediately below such a sea ice layer will experience intermittently dimmer and brighter light levels as they advect underneath these smaller areas of more and less optically transmittant sea ice. Phytoplankton at depths generally below any apparent subsurface maxima in PAR will presumably experience weaker fluctuations, given that this effect increasingly averages out spatially with increasing depth. For phytoplankton at those depths and below, these ITP-observed PAR levels might reasonably serve as acceptable spatial averages and be suitable for use in broader scale biogeochemical models. For the phytoplankton in the upper water column however, which experience the highest shading-driven fluctuations in PAR, this issue of light fluctuations introduces considerable complexity in terms of modeling primary production. This fluctuating light effect due to water column moving underneath a spatially variable surface ice layer is analogous to the long studied yet still unresolved phenomenon of phytoplankton in ice-free waters experiencing light fluctuations due to the passage of intermittent clouds. Passing clouds can introduce similarly intermittent and rapid perturbations in the light levels phytoplankton experience (Stramska and Dickey 1998), but to date there are no good physiological models to quantify how such fluctuations affect photosynthesis, nor is there even agreement as to whether such fluctuations in light increase or retard light harvesting (e.g., Walsh and Legendre 1983;Qu eguiner and Legendre 1986;Pearcy et al. 1994). For phytoplankton cells in the top optical depths underneath a sea ice layer with light transmittance properties that are spatially heterogeneous, a spatially averaged light intensity might not be appropriate for use in biogeochemical models, but at present there is no adequate alternative. The water column immediately under sea ice provides a new system to study for future efforts seeking to understand how dynamical aspects of photosynthesis in fluctuating light environments affect photosynthesis and primary production, complementing other previous examined mechanisms for introducing fluctuations such as vertical displacement and intermittent cloud cover.

Implications for measuring and modeling PAR profiles under sea ice
Our main finding that vertical distributions of PAR under Arctic sea ice so frequently deviate from the canonical exponential has an important corollary: that knowing the actual vertical distribution of PAR is critical for recovering water column properties such as k D . This has ramifications for how autonomous approaches can and should be used to measure under-ice light fields in Arctic Ocean observational programs, if one goal is to examine euphotic zone depths over seasonal scales. In ice-free water columns, k D can be derived from as few as two separate radiometers placed at different depths (e.g., Zheng et al. 2002) but such an approach becomes problematic if the light profile does not decrease monotonically and exponentially with depth, as is often seen here under sea ice. Two sensors at separate fixed depths would be insufficient to resolve even the basic vertical structure of the under-ice PAR profiles seen in these ITP time series, and any estimates of k D from such data would likely include considerable bias and/or error. From a strict mathematical perspective, four sensors at different depths would be the minimal number needed to constrain the fourparameter model we present in Eq. 3, although the specific choice of depths would also strongly affect the robustness of k D estimates that are derived from four discrete-depth PAR measurements. The observed seasonal evolution of vertical distributions of PAR (as seen in ITP64's time series, Fig. 4g) suggests that in fact there may be no optimal fixed depths for placing a small number of discrete radiometers that will Laney et al.
The euphotic zone under Arctic sea ice ensure robust assessment of k D over periods of many months. Innovations that use optical multiplexers to measure PAR at a greater number of fixed depths (e.g., 10 or more as described by Wang et al. 2014) may improve discrete sampling of PAR when its vertical distributions and temporal variability are as complex and temporally variable as is seen here. A chain of discrete radiometers could also enable PAR sampling at a larger number of fixed depths, with costs proportional to the number of depth sampled. The I PAR (z) profiles in these ITP-enabled PAR time series had considerable spatial resolution in the vertical ($25-30 cm), but even with such highly resolved profiles an appropriate model is still is necessary to estimate common euphotic zone metrics such as the depth of the 1% light level. Computing the 1% light level depth requires robust estimates of either I 0 or k D , neither of which can be reliably determined from our observed under-ice PAR profile time series using traditional models such as Lambert-Beer (e.g., Eq. 1). It is tempting simply to ignore in these PAR profiles those samples that occur above the apparent subsurface maximum, but this surface shading effect influences PAR at depths well below that of any apparent maximum (see Fig.  7b) and so it is mathematically inappropriate to compute k D from PAR profiles truncated in such a way. Other ad hoc, model-independent approaches to eliminate those samples affected by surface shading similarly risk introducing new artifacts into estimates of k D . This concern motivated our development of a more rigorous, model-based approach to parameterize these observed near-surface anomalies in PAR, even if our eventual solution used an empirical representation of the near-surface shading effect instead of having a specific mechanistic basis to explain the phenomenon we observed.
This four-parameter semi-empirical model (Eq. 3) provided a quantitative means to separate the effects of water column attenuation on I PAR (z) from any shading effect introduced by local spatial heterogeneity in the transmittance of the overlying sea ice. It represents an expansion of the Lambert-Beer model that has long been used to model water column light attenuation (Preisendorfer 1976), modified with an empirical description of the near-surface shading phenomenon arising in the surface ice layer. This empirical description provides acceptable fits to observed I PAR (z) even when the corresponding MSE is comparatively high within these data and when the nominal quality of fit in turn was comparatively low (Fig. 6e), suggesting that it reasonably represents the PAR profiles we observed in these ITP time series. We chose to forego examining this nearsurface shading effect in the context of any specific physical mechanism such as the melt-pond dynamics modeled by Frey et al. (2011) because there was no reason to expect that any single mechanism would be responsible for the entire range of I PAR (z) distributions that we observed in these ITP time series. In fact, when the Frey et al. (2011) melt pond model and this semi-empirical model of Eq. 3 were allowed to compete on equal terms to best fit these time series-in spite of our concerns with using that model in such a manner-the melt-pond model provided a better statistical representation of these observed I PAR (z) distributions in only $50% of all profiles compared to Eq. 3 (Table 4). Were more complete PAR profiles available, i.e., profiles that extend fully to the ice-ocean interface, it is possible that the Frey et al. (2011) model would demonstrate a much higher statistical advantage against Eq. 3, although the concern of parameter interdependence between N and P would still remain. Future efforts to improve the parameterization of under-ice irradiance distributions will likely require the continued assessment of multiple models to explain observations, in order to avoid introducing biases arising from adopting a single mechanistic model a priori without firm proof of its appropriateness.
One benefit of developing and using Eq. 3 for parameterizing under-ice I PAR (z) is that it preserves existing definitions of k D and I 0 , and so estimates of these parameters derived using this model can be compared directly with other estimates derived using Eq. 1 in cases where it was appropriate. We exploited this in our own analyses here, to compare estimates of k D and I 0 generated by both Eqs. 1, 3 in our model fitting and selection exercise (Fig. 6). It can also be exploited to compare our estimates of k D and I 0 directly with the same parameters generated by other workers who measured PAR profiles under Arctic sea ice. Of these two parameters, I 0 is more difficult to compare between studies. Its absolute magnitude in ice-covered water columns is a function of solar elevation and a number of atmospheric and sea ice properties that make any comparison among published values possible in only a very qualitative manner. In contrast, k D can typically be readily compared between studies, being a quasiinherent optical property whose absolute magnitudes are determined primarily by the water column itself. A recent summertime study in the Eurasian Basin above 858N used a hand-lowered sensor to collect PAR profiles through sea ice cover (Lund-Hansen et al. 2015), and PAR profiles from this region of the central Arctic Ocean were observed to be predominantly exponential in shape, similar to what was seen here with the two ITPs deployed at those highest latitudes (ITP48 and ITP72). Estimated magnitudes of k D from that study, properly derived from a Lambert-Beer model, were 0.172 m 21 (SD 5 0.0734, n 5 25), comparable to the $0.15 m 21 values we derived here from autonomously sampled profiles and using our model fitting and selection approach (Fig. 6b).
For an analogous comparison of k D values in the Beaufort Gyre, independent estimates can be derived from hydrocasts performed during the BGEP cruises on which these ITPs were deployed, as a limited number of those casts included a PAR radiometer. Estimates of k D derived from these hydrocasts were $0.072 m 21 (SD 5 0.019, n 5 6), comparable in magnitude to those estimated from these Beaufort Gyre ITP time series and derived primarily from Eq. 3 (Fig. 6a). Interestingly, these ship-based PAR profiles illustrated the opposite of the near-surface shading effect that is predicted by Eq. 3, and showed instead a near-surface enhancement of I PAR (z) in the topmost depths (Fig. 8). Such a near-surface enhancement can be expected for situations such as this one where an icebreaker clears ice from around itself in order to conduct hydrocasts. By clearing surface ice in its immediate area, the icebreaker in effect increases local light transmittance in the area above the vertical PAR profiling axis, and areas more distant from the ship remain more ice covered and therefore less optically transmittant. These ship observations provide the only examples in this study where PAR profiles are consistent with a near surface enhancement in PAR and Eq. 3 having an k NS less than zero. Although we did not observe negative estimates for k NS with the ITP-derived PAR profile time series, it might be unreasonable to expect to see negative values in those situations given that ITPs are typically sited on the thickest ice available on a floe in order to maximize their lifetimes when the ice eventually ridges or melts out in the future. Therefore, ITP deployment biases alone may be responsible for our seeing only positive estimates of k NS in these time series. Regardless, these icebreaker PAR profiles illustrate an interesting polar counterpoint to the better known concern with "ship shading" when making PAR profiles in water columns free of sea ice. In ice-covered water columns, although the ship itself is a source of shading, the net effect of the presence of the ship can be to increase the overall illumination of the water column below by admitting more light into the ocean than the ship itself blocks, i.e., introducing a "ship unshading" effect.
The general agreement between estimates of k D from our curve fitting and model selection approach, and independent estimates derived from more traditional methods, suggests that k D can be reasonably derived from autonomously collected, vertically complex PAR profiles, even those exhibiting the considerable spatiotemporal variability as was seen in the ITP-enabled time series examined here. These results should be viewed as preliminary and further refinement will require a larger and more comprehensive data set of in situ observations of I PAR (z) in under-ice Arctic water columns. In developing Eq. 3, one assumption we made was to characterize the water column using a single diffuse attenuation coefficient (Stavn 1988;Gordon 1989). We recognize that distributions of optical constituents under ice can exhibit considerable vertical structure, for example, as was seen in chlorophyll distributions during these ITP deployments . This assumption of a single k D is also shared by the melt pond model Frey et al. (2011) and can be considered a reasonable first-order approximation when developing new parameterizations. Future models for underice I PAR (z) might be improved by incorporating any concurrently measured bio-optical variables that should affect k D such as chlorophyll, non-algal particles, or CDOM, to account for vertical heterogeneity in water column attenuation. However, it should be noted that in these ITP time series the actual chlorophyll concentrations that were seen under sea ice were typically quite low, less than 0.5 lg L 21 in most cases. Moreover, given the characteristics of these particular fluorometers, these values may represent at least a twofold overestimate in chlorophyll when applying the calibrated values provided by the manufacturer (Boss and Ha€ entjens 2016), which is the nominal calibration approach used for these ITP deployments. It is certain that a single k D is an inappropriate assumption from a strictly theoretical perspective, but this data set is not adequate for a more detailed assessment of whether or not multiple attenuation coefficients will meaningfully improve these autonomously obtained estimates of euphotic zone depth.

Implications for observational programs assessing variability in under-ice PAR
For any initial assessment of euphotic zone depths under Arctic sea ice over broad spatiotemporal scales, known concerns with using a percentage of surface light level to demarcate the euphotic zone (Banse 2004) can be overlooked to the first order. Percent light levels are readily determined using the approaches outlined above, but other suitable metrics for euphotic zone illumination, such as the depth of a casts in the Beaufort Gyre during an ITP deployment cruise. Lines indicate best-fits of Eq. 1 over the corresponding range of PAR. A deviation from exponential is apparent in all four profiles near the surface, indicating an enhancement in PAR that presumably reflects a "ship unshading" effect of an icebreaker clearing sea ice cover immediately adjacent to the ship. Deeper in the water column the PAR profile becomes roughly exponential, while the deviation at depth reflects rolloff in sensitivity of the sensor as PAR values become unmeasurable. For comparison, the vertical dashed line indicates the empirically estimated threshold of sensitivity for the PAR radiometers used with the ITPs.
given isolume (where daily integrated photon flux is constant, e.g., Letelier et al. 2004), can also be generated from these ITP-enabled time series. Day lengths vary substantially over the course of the year in polar regions but these ITPs profile multiple times per day and thus also provide valuable diel-scale measurements of under-ice PAR that would be necessary for computing such alternate metrics for euphotic zone illumination.
Only one of the PAR profile time series in this ITP study spanned an entire annual period (ITP64), in part reflecting the challenges with autonomous systems surviving over harsh Arctic winter ice conditions (Toole et al. 2011). These observations from ITP64 indicated that strong perturbations in under-ice PAR can occur very early in the growth season (Fig. 4g), which we interpreted here as representing the opening and later refreezing of leads in the overlying sea ice (Assmy et al. 2017). These springtime observations also indicated that this near-surface shading phenomena in PAR can arise very early in the growth season and that their characteristic vertical shapes can vary over time well into the summer. This long-term temporal variability presumably reflects ongoing changes in sea ice properties, although we recognize the potential for indirect buoy-driven shading artifacts due to drifting snow and melting ice (see Results). Weekly-to seasonal-scale changes in sea ice driving this temporal evolution in the shape of I PAR (z) may involve changes in specific properties of a particular physical mechanism responsible for these near-surface shadings, for example, changes in the number and areal extent of melt ponds in the framework of the Frey et al. (2011) conceptual model during portions of the year when melt ponds are present. Alternatively, this temporal variability in I PAR (z) might reflect changes in the relative contribution of disparate mechanisms that introduce spatial heterogeneity in the light transmittance of the sea ice layer, such as a transition from patchy snow cover to patchy melt pond cover. An improved understanding of the optical properties of the sea ice layer and their influence on light profiles in the underlying water column might provide a means to infer such temporal trajectories in the character of the under-ice light field from measured proxies of certain sea ice parameters, instead of needing to rely solely on under-ice PAR measurements. Unfortunately, in this study, no other time series data of relevant sea ice properties are available to complement these ITP PAR profiles, and so we have little insight into the physical mechanisms responsible for these observed temporal changes in the shape of I PAR (z). Five of these eight ITPs were installed on the same floes with Ice Mass Balance (IMB) buoys, four of the "standard" IMB design (Richter-Menge et al. 2006) and one Seasonal IMB (SIMB; Polashenski et al. 2011). Ideally these IMB-ITP pairs would have provided companion time series data on sea ice thickness and snow cover, but usable co-temporal data were available in only one of these five instances and only for a period of a few months. These time series were insufficient to support any quantitative assessment of physical drivers of the trends observed in the shape of ITP64's PAR profiles.
Future studies to elucidate the specific physical mechanisms that affect the magnitude of I PAR (z) and its vertical character under sea ice would presumably require comparable time series of ice layer properties such as thickness and snow cover. Other relevant properties to examine include the absolute incident insolation above the ice layer and the presence and areal extent of melt ponds and leads. Measuring such variables in the Arctic environment over the necessary time scales involves considerable technical challenges, even with a parameter as nominally simple to measure as the incident insolation. State-of-the art approaches for measuring insolation above sea ice (e.g., Nicolaus et al. 2010b;Hudson et al. 2012) might not be robust enough to provide acceptable data over periods longer than a few months. Melt ponds and leads are likely to be important factors in driving variability in I PAR (z) in ice-covered regions of the Arctic Ocean, and future studies might consider implementing above-ice camera systems (e.g., Sankelo et al. 2010) to provide imagery for quantifying the melt pond fraction P. Direct measurements of P would provide valuable insight into the timing and influence of melt ponds on I PAR (z) and would also help to address the parameter independence issue that was encountered here when applying the Frey et al. (2011) melt pond model to these autonomously-observed I PAR (z) time series where both P and N were unconstrained.

Future directions
Although the optical and bio-optical sensor suite used in this study performed for periods up to a year in Arctic under-ice environments, as with any prototype system several technical issues were identified that merit further attention. These sensor suites used commercial, off-the-shelf sensors and incorporated rudimentary strategies to monitor sensor behavior and failure (see Laney et al. 2014). This limited self-diagnostic capability of these commercial sensors and the integrated sensor suite makes it difficult to ascribe specific modes of failure to the sensors themselves, to their supervisory systems, to the connecting cables, to mechanical or mounting issues, to the performance of the shutter actuator, or to other issues not envisioned pre-deployment. In future efforts to obtain similar long-term, unsupervised observational data in these ice-covered ecosystems, more robust and comprehensive self-diagnostic capabilities should be strongly considered and if unavailable commercially, developed de novo.
In this study, we were also concerned with the accuracy of these PAR measurements, being aware that sensor drift and fouling affect both the relative and absolute accuracy of these PAR measurements in long-term observational studies. For a subset of these time series, we do have initial, independent ship-based profiles of chlorophyll and less frequently PAR, obtained by hydrocasts, but these profiles were conducted within tens of kilometers of the buoy deployment site and often many hours earlier or later than the first programmed profile of any nearby ITP. Any correlation between the ITP and ship PAR profiles, for use in establishing initial in situ cross-calibrations, is therefore limited. The magnitudes of sensor drift and fouling are similarly difficult to confirm throughout each deployment given the numerous and often unexamined factors that control the absolute intensity of light under ice. These PAR sensors are calibrated with two parameters accounting for offset and gain, and we are in principle able to track any change in offset above the noise floor. It is not possible however to track changes in gain over time with these radiometers. As for fouling, our sensor cluster design includes a protective copper shutter over the optical faces of both the PAR sensor and fluorometer which is opened only during active profiling and remains closed at all other times . We believe this shutter is the best-possible solution for mitigating fouling over these long-term deployments, but it does not provide a robust solution for assessing sensor drift, either inherent drift due to component ageing or apparent drift due to longerterm mechanical changes to these radiometers in properties such as its cosine diffuser or optical windows.
This 4-yr Arctic bio-optical profiling project with ITPs paralleled initial efforts to obtain bio-optical profiles autonomously in the open ocean using profiling floats (e.g., Boss et al. 2008). In the open ocean, an important objective is to obtain highest-quality profile data of bio-optical properties for eventual calibration-validation applications with ocean color remote sensing (Claustre et al. 2010;IOCCG 2011). With virtually all of our observations being made in icecovered water columns unobservable from space, we did not adopt many of the more sophisticated quality control approaches that have been developed for processing openocean optical and bio-optical profile data for remote sensing purposes (e.g., Xing et al. 2011;Xing et al. 2012;Organelli et al. 2016). Nevertheless, many of the procedures developed in those studies are appropriate for application with underice data and should be considered in future Arctic efforts when under-ice optical measurements become more routine. Other aspects of these open-ocean procedures, such as accounting for surface wave focusing, are relatively unimportant in the under-ice environment and can be ignored. Additionally, there are features unique to the under-ice environment that are not incorporated in existing protocols for processing optical and bio-optical profile data in ice-free water columns, such as profiler slippage on the tether, which merit attention in any future effort developing quality control procedures for ITP optical data. Routine assessment of euphotic zone variability is still in its infancy in the icecovered Arctic Ocean, but the findings of this study suggest that autonomous systems can play a meaningful role in bettering our understanding of the under-ice light field in icecovered Arctic marine ecosystems.