Model for deriving benthic irradiance in the Great Barrier Reef from MODIS satellite imagery

: We demonstrate a simple, spectrally resolved ocean color remote sensing model to estimate benthic photosynthetically active radiation ( bPAR ) for the waters of the Great Barrier Reef (GBR), Australia. For coastal marine environments and coral reefs, the underwater light field is critical to ecosystem health, but data on bPAR rarely exist at ecologically relevant spatio-temporal scales. The bPAR model presented here is based on Lambert-Beer’s Law and uses: (i) sea surface values of the downwelling solar irradiance, E s ( λ ); (ii) high-resolution seafloor bathymetry data; and (iii) spectral estimates of the diffuse attenuation coefficient, K d ( λ ), calculated from GBR-specific spectral inherent optical properties (IOPs). We first derive estimates of instantaneous bPAR . Assuming clear skies, these instantaneous values were then used to obtain daily integrated benthic PAR values. Matchup comparisons between concurrent satellite-derived bPAR and in situ values recorded at four optically varying test sites indicated strong agreement, small bias, and low mean absolute error. Overall, the matchup results suggest that our benthic irradiance model was robust to spatial variation in optical properties, typical of complex shallow coastal waters such as the GBR. We demonstrated the bPAR model for a small test region in the central GBR, with the results revealing strong patterns of temporal variability. The model will provide baseline datasets to assess changes in bPAR and its external drivers and may form the basis for a future GBR water-quality index. This model may also be applicable to other coastal waters for which spectral IOP and high-resolution bathymetry data exist


Introduction
Benthic irradiance is defined as photosynthetically active radiation (PAR, or benthic PAR) at 400 -700 nm wavelengths that reaches the seafloor.Benthic PAR provides the primary energy resource that drives benthic photosynthesis, thus essentially defines primary production in the benthic coastal ocean [1].Previous studies have emphasized the importance of light availability for benthic habitats.For example, light limitation determines photosynthesis, growth, and distribution of seagrasses [2][3][4][5], depth limits for coral reef development [6][7][8][9], and vulnerability of coral reefs to ocean acidification as a result of reduced calcification [10,11].On the other hand, high light levels in combination with heat stress from elevated seawater temperature can also lead to coral bleaching [12,13].Such juxtaposition thus highlights the need for an improved understanding of the underwater light field and benthic light availability, to better assess and manage marine ecosystem health [3,8].Efficient management of these benthic ecosystems would greatly benefit from a synopticscale understanding of spatial and temporal patterns of benthic light availability.Unfortunately, benthic light availability remains poorly understood [1,14], because datasets that could provide the much-needed information are either insufficient or lacking, including in one of the largest coral reef ecosystem on Earth, the Great Barrier Reef (GBR).As such, our main objective is to demonstrate the feasibility of estimating benthic light values within the GBR region using a remote sensing model to bridge this knowledge gap.
The GBR Marine Park (GBRMP) is a world heritage listed site adjacent to northeastern Australian coast spanning from 10°S to 24°S latitude.It is located on a continental shelf with an area of ~344,000 km 2 .While direct observations of the underwater light field have historically been measured by ships of opportunity and in situ moorings in the GBR, these existing data are spatially and temporally sparse.Although PAR in surface waters is frequently reported for marine ecosystems [15,16], benthic PAR is usually available only at a few point locations over relatively short periods [17][18][19].To understand the responses of whole ecosystems at regional scales, we need to estimate environmental conditions, including benthic PAR, at relevant scales (in this case, the scale of the whole GBR).Ocean color remote sensing presents an option for observing benthic PAR in the GBR at ecologically relevant spatio-temporal scales, which can provide near-daily images of the entire GBR at a nominal pixel resolution of 1 km 2 .While a number of ocean color satellite processing algorithms have been developed for monitoring optically complex waters of the GBR [20][21][22][23], we note that none have focused on deriving benthic light availability.
Benthic irradiance, here denoted as E b (λ), is essentially regulated by three factors.First, the downwelling irradiance incident on the ocean surface, E s (λ,0 + ) (W m −2 ), which is subject to atmospheric processes (e.g., gaseous absorption and molecular scattering) that diminishes light as it makes it way to the sea surface.Second, the absorption and scattering processes within the water column [14,24,25] that depends on wavelength and have magnitudes and shapes that are proportional to the relative concentration of optically-active constituent matter [26] including: phytoplankton, non-algal particles (NAP) and colored dissolved organic matter (CDOM).It is important to note that IOPs are impartial to the direction and magnitude of the ambient light field [27].In coastal regions, IOPs can be regulated by a number of drivers including, but not limited to: run-off and resuspension of fine sediments, nutrients and organic matter due to wave-induced vertical mixing, human activity such as dredging and dredge-spoil placement on the seafloor, and biological processes related to phytoplankton growth.Collectively, these drivers make the inversion of sensor-observed remote sensing reflectance, R rs (λ) (sr −1 ), to derive IOP values, a non-trivial problem.Third, and also important is the water column depth which defines the distance of light penetration.
In optically shallow waters like those of the GBR, light reflected off the seafloor is another factor that limits the application of contemporary ocean color algorithms developed for deep, open-ocean waters [22,28].Fortunately, a GBR-specific IOP inversion model (Shallow Water Inversion Model, SWIM) has recently been developed [22] that improves the derivation of spectral estimates of IOPs, specifically, total absorption coefficients, a(λ) (m −1 ), and total backscattering coefficients, b b (λ) (m −1 ), for any pixel within the GBR where R rs (λ) are available.We note that R rs (λ) is a fundamental sensor-observed radiometric variable used in ocean color algorithms.SWIM is described in detail elsewhere (see [22]).In brief, SWIM takes advantage of two mapped GBR-specific ancillary datasets, namely: (i) a high-resolution shelf bathymetry map [29] matching the resolution of the remote sensing datasets, and (ii) a bottom substrate brightness (albedo) map [30].Using these ancillary datasets, SWIM can generally correct for optically shallow effects.Note however, that for more optically complex waters, the skill of the SWIM model may still be challenged in deriving IOP values.Nonetheless, the IOP estimates that can be obtained via SWIM are very valuable and sufficient for our current purpose.By using SWIM-derived IOPs, it is thus possible to spectrally characterize near-daily GBR-wide light attenuation and, in turn, estimate benthic light availability (i.e., PAR) using ocean color remote sensing.
A contemporary approach for estimating PAR at depth is to use the broadband, watercolumn averaged diffuse attenuation coefficient of PAR, K d (PAR) (m −1 ), and the Lambert-Beer's Law which assumes K d (PAR) is constant over the optical pathlength (depth).However, in practice, K d (PAR) exhibits strong depth dependence, particularly in the upper water column, even in well-mixed waters [31,32].Thus, using a K d (PAR)-based method may lead to inaccurate estimates of PAR at depth.We note that usable solar radiation (USR), a spectrally integrated broadband term (integrated over the 400 -560 nm range), has recently been proposed as a promising alternative to PAR [31] as its diffuse attenuation coefficient, K d (USR) (m −1 ), exhibits less depth dependence.However, most ecological studies report their findings (e.g., benthic light stress thresholds) in terms of PAR.Here, we described a remote sensing approach to estimate benthic PAR that does not utilize K d (PAR).
In this study, we demonstrate a simple physics-based ocean color remote sensing model to estimate two spectrally integrated benthic light variables, namely: (i) instantaneous PAR at the seafloor (hereby denoted as bPAR i ) (μmol photons m −2 s −1 ), and (ii) daily integrated PAR at the seafloor (hereby denoted as bPAR d ) (mol photons m −2 d −1 ).Estimates are derived from (i) E s (λ,0 + ), (ii) high resolution shelf bathymetry, and (iii) spectral estimates of IOPs obtained using SWIM [22].Here, E s (λ) and IOPs are derived from data collected by NASA's Moderate Resolution Imaging Spectroradiometer aboard the Aqua spacecraft (MODISA).We then validate the model by comparing concurrent satellite-derived and observed in situ benthic PAR values from four test sites within the GBR.Finally, results, as well as limitations and future applications of the model within the GBR region are discussed.

Benthic irradiance model description
To derive E b (λ), the total amount of subsurface downwelling planar irradiance, E d (λ,0 -) (W m −2 ), transmitted to a depth z (m) under a clear sky, we used the Lambert-Beer's Law: where K d (λ) (m −1 ) is the spectral diffuse attenuation coefficient.We emphasize that in our approach, spectral dependence was kept throughout all calculations.Satellite-derived E s (λ,0 + ) (W m −2 ) represents the downwelling planar irradiance when the Sun is at zenith and just above the air-sea surface.Here, E s (λ,0 + ) was obtained by employing the standard NASA atmospheric transmittance model as described in Mobley, et al. [33] which assumes that the cloud/surface system is non-reflecting and non-absorbing.Because E s (λ,0 + ) is estimated just above the sea surface, it is necessary to propagate this parameter across the air-sea interface to derive E d (λ,0 -) [25,34].Here, we followed Baker and Smith [35] and used: where t g is the global transmittance of the air-water interface [36].Essentially, the total downwelling global irradiance just above the interface can be decomposed into direct ("sun") and diffuse ("sky") components, both of which are functions of wavelength and the corresponding solar zenith angle in air [35] and t g is equal to where ρ sky and ρ sun are the sea surface reflectance of the diffuse and the direct irradiances, respectively.The value of ρ sky is roughly equivalent to 0.066 for cloud-free conditions [35] while the ρ sun can be approximated by employing the Fresnel reflectance equation for unpolarized light [27] written as: where θ solz (degrees) is the solar zenith angle in air at the time of satellite overpass and θ t (degrees) is the solar zenith angle of light transmitted downward in the water calculated as: where n w is the real refractive index of seawater approximately equal to 1.34 [27].For cases where θ solz ≠ 0.0 (e.g., sun not overhead), Eq. ( 3) was used.Otherwise, for normally incident light, the reflectance of the direct irradiance was calculated as: We note that in this study we have assumed cloud-free, clear sky conditions for all calculations and neither considered the effect of a wind-roughened surface in our air-sea transmittance calculations.Further details on deriving the transmitted irradiance across the interface are described elsewhere [35,36].
An IOP-based method for deriving K d (λ) [37] was employed as this approach is more robust in optically complex coastal waters relative to simple empirical methods [38][39][40][41].The IOP-based K d (λ) model we have used in this study is described and evaluated in detail elsewhere [37,42] and has a general form: where a(λ) and b b (λ) are the coefficients of spectral total absorption and spectral total backscattering, respectively.An interesting feature of this model [37] is that it allows estimation of the attenuation coefficient for any wavelength with high accuracy, given the coefficients for a(λ) and b b (λ) for any region-of-interest can be obtained with suitable accuracy.For our purpose, we obtained location-and time-specific estimates of a(λ) and b b (λ) from the SWIM algorithm and implemented Eq. ( 6) for each pixel location within our test sites to obtain K d (λ).We note however, that the K d (λ) model above can utilize IOPs derived from other IOP algorithms and is hence not tied to the SWIM algorithm.

Instantaneous benthic PAR
Values of bPAR i were calculated by spectrally integrating the values of E b (λ) as: 700 400 ( ) .
Note that the spectral range (as further detailed later) is not continuous and corresponds to the MODISA ocean color band centers.The integration was performed in Python 3.6.3using the numerical python (numpy) library's "trapz" function that implements the trapezoidal method to numerically approximate the area under the curve using definite integral over discrete intervals (i.e., here defined by the band centers of MODISA).

Daily integrated benthic PAR
To obtain bPAR d from sunrise to sunset for any given day and pixel location, we first estimated the instantaneous benthic PAR at noon, bPAR n (μmol photons m −2 s −1 ), which gives the maximum potential PAR irradiance at the seafloor, using: where the subsurface downwelling planar irradiance is multiplied by the cosine of θ solz (the term enclosed in square brackets).Assuming clear skies (i.e., no clouds) and that the cloud/sky system is non-reflecting and non-absorbing, the magnitude of E s (λ,0 + ) can be estimated as a function of θ solz and is approximately sinusoidal in relation to the daylight hours and symmetric about noon [25].Further, we can compute the solar zenith angle for a given solar hour angle, φ, as: ( ) cos sin( )sin( ) cos( )cos( )cos( ) where, for a given pixel on a single day latitude, lat (radians), the solar declination angle, δ (radians) is computed as: 0.006918 0.399912 cos( ) 0.070257 sin( ) 0.006758cos(2 ) 0.000907 sin(2 ) 0.002697 cos(3 ) 0.00148sin(3 ) The term γ (radians) is the date expressed as an angle: 2 * ( 1) 365 jday where jday is the Julian day (0-365 in regular year, or 0-366 in leap year).Thus, bPAR d can be expressed as the double integral: where φ rise and φ noon = 0 denote hour angles at sunrise and noon, respectively.The integral with respect to solar hour angle in Eq. ( 12) is computed numerically using a trapezoidal rule by discretizing the solution into fifty φ steps ranging from φ rise to φ noon = 0, where φ rise is computed as: ( ) 3. Data and model evaluation

Satellite data processing
MODISA raw Level-1A (L1A) data files spanning 2002 -2018 were downloaded from NASA (https://oceancolor.gsfc.nasa.gov/)and spatially subsetted to encompass four test sites within the GBR (black filled circles in Fig. 1).These data files were then processed to geophysical Level-2 (L2) data products using NASA's Ocean Color Software Suite (OCSSW) processing code that is distributed as part of the SeaDAS display and analysis software package [44].L2 data products were produced for ten visible MODISA bands centered on: 412, 443, 469, 488, 531, 547, 555, 645, 667, and 678 nm and stored in NetCDF4 files.These L2 data products included: E s (λ,0 + ), θ solz , and SWIM-derived IOPs, a(λ) and b b (λ).NASA's standard ocean color atmospheric correction was used to produce all L2 data products [45][46][47].The data processing yielded nearly 5500 daily MODISA L2 files for each of the four test sites.

GBR bathymetry data
Bathymetric data used to define the water column depth (datum: mean sea level (MSL)) were extracted from a gridded high-resolution bathymetry and digital elevation model (DEM) of the GBR, 3D-GBR [29].This DEM dataset has a grid pixel resolution of ~100 m x 100 m that can resolve fine-scale details of the bottom topography of the GBR reefs and inter-reefal systems (see Fig. 1).The full GBR bathymetric dataset is available from Australia's eAtlas website and environmental research data mapping and management system at https://eatlas.org.au/data/uuid/200aba6b-6fb6-443e-b84b-86b0bbdb53ac.

In situ data collection and processing
The Australian Integrated Marine Observing System (IMOS) maintains an array of observing sub-facilities [48].One of the regional Shelf Mooring Array sub-facilities is located on the coast of Queensland, which is part of the Queensland IMOS (Q-IMOS) node established in 2007 and currently operated by the Australian Institute of Marine Science (AIMS).To conduct validation for this study, four existing mooring stations: Yongala, Palm Passage, Myrmidon, and Heron Island South (referred to as Heron throughout the manuscript) were expanded in 2015/2016 to also include PAR sensors, all located >1 km away from reefs and other geometrically shallow features (Table 1).These four mooring stations were therefore selected as in situ benthic PAR validation sites.
It is worth noting that the sites used for validation represent spatially and optically diverse water bodies.The Yongala mooring site is located inshore with shallow lagoon waters often subject to wind-driven resuspension of bottom sediments and terrestrial influence.Hence, it encompasses both turbid and moderately clear seawater conditions and strong variation in optical properties (>20-fold range in PAR values).It is also exposed to seasonal runoff from the nearby Burdekin River.The Palm Passage and Myrmidon test sites are both located on the edge of the continental shelf, and although they are seasonally exposed to the influence of the East Australian Current (EAC), the southwest Pacific poleward western boundary current, these two stations most typically represent clear oceanic optical conditions (four-to six-fold range in PAR values), typical of oligotrophic waters.Meanwhile, the Heron station represents relatively stable intermediate in-water optical conditions (clarity and variability) although it may also be influenced by other important local oceanographic features (e.g., eddy-driven intrusion and productivity, tidally-driven flows, seasonal Trichodesmium blooms), yet not much more than the effects of the EAC oceanic intrusions experienced by the two deep oceanic test sites.Overall, the selected test regions comprise a range of optical water types within the GBR lagoon, and hence, suitable to test the capability of the newly developed benthic irradiance model over an optical dynamic range found in this region.To measure in situ benthic PAR at Yongala, a Seabird Scientific SBE16PLUS v2 SEACAT profiler with upward facing WETLabs Environmental Characterization Optics (ECO) PARSB sensor was deployed 0.5 m above the bottom substrate at ~30 m water depth.At the other three mooring sites, WETLabs ECO PARSB sensors were clamp-mounted to a permanent mooring wire at a fixed distance from the seafloor at each deployment period.
PAR sensors collected 5-second data bursts every 15 minutes (480 benthic PAR records per day) at all stations except Yongala where single measurements were collected every 15 minutes (96 benthic PAR records per day).Deployments and servicing of sensors occurred approximately every six months.After each deployment period, data were downloaded and the instrument's optical component was checked, characterized and tested to ensure the quality and validity of the data between deployments.For each recovery period, the data was analyzed and quality controlled such that: (i) data records from the beginning and end of each deployment were excluded to ensure that only stable PAR measurements were included in subsequent analysis or within a complete day cycle in the case of bPAR d validation, (ii) data points when instrument failure was experienced (e.g., battery problems) were also excluded (e.g., Palm Passage deployments 2 and 3), and (iii) night-time values were forced to zero by applying an offset based on the dark count readings of the sensor for each deployment period.Table 1 lists the deployment details and total number of in situ PAR data remaining after quality controls were applied.

Algorithm validation approach
The model performance was evaluated via matchups between: (i) concurrent instantaneous parameters, bPAR i and in situ measured instantaneous PAR near the seafloor (Yongala) or in the water column (Myrmidon, Palm Passage and Heron), (ii) concurrent daily-integrated parameters, bPAR d and the daily-integrated in situ measurements, and (iii) all four test sites combined against overall satellite retrieval performance for both instantaneous and daily integrated cases.
For clarity, because the PAR loggers, except in Yongala, were mounted in mid-water rather than on the bottom, the model was also run for the depths that complemented the in situ data for each site and deployment period.Consequently, for brevity, we use the term 'benthic PAR' hereto for both PAR near the benthos (i.e., Yongala) or at specific water column depth (i.e., clamp-mounted logger depths at Palm Passage, Myrmidon and Heron).

Matchups of satellite-derived to in situ data pairs
For each L2 file containing results of the benthic irradiance model (either the instantaneous bPAR i and daily integrated bPAR d ) corresponding to an in situ measurement, satellite-derived values were extracted for a 3 x 3 pixel box centered at the mooring station locations (see Table 1 for coordinates).Next, image pixels were discarded if any of the pixels were associated with the following NASA L2 quality flags: LAND, CLDICE, HILT, HIGLINT and STRAYLIGHT [49] for quality assurance.Extracted pixels were also discarded using two additional quality control procedures: (i) if any of the 9 pixels within the test box had been masked [50], and (ii) if the coefficient of variation (CV, standard deviation/mean) between the nine pixels of the satellite-derived bPAR i or bPAR d was greater than 15%, indicating low intra-pixel stability of the modeled parameter over relevant time-scale [51].Finally, common statistical parameters including minimum, maximum, mean, standard deviation, and median values were calculated for each nine-pixel box as a final check; however, we only reported and used median values in further analysis.Matchups between concurrent data pairs of satellite-derived and in situ measured bPAR were then performed using quality controlled median bPAR i (and bPAR d ) timeseries compared against the in situ bPAR timeseries.To qualify as 'concurrent', we included in situ data from ±3 hours around the MODISA satellite overpass time similar to previous studies [52].Matchups were conducted using R Studio (version 1.1.456)running R version 3.5.2[53] (codes available upon request).

Regression analysis and other performance metrics
Type-II linear regression analysis was used to compare the concurrent satellite-derived and in situ measured bPAR values per test site and on the combined dataset to assess how well the model performed over the dynamic range of benthic PAR in the GBR.The regression slope, intercept, and coefficient of determination (r 2 ) were calculated on log 10 -transformed data.Reduced major axis regression analysis was completed using lmodel2 package [54] in R version 3.5.2.Model type II linear regression was chosen due to assumed inherent error in both the satellite and in situ datasets (e.g., because pixels from satellite imagery are not in exactly the same location and time as the in situ data, etc.).Lastly, additional metrics, including bias and mean absolute error (MAE), were calculated on log 10 -transformed data as measures of overall algorithm bias (i.e., systematic error) and accuracy (i.e., average model prediction error) that were back-transformed out of logarithmic scale before interpreting the results.Following Seegers, et al. [51], bias and MAE were calculated as: 10 10 where bPAR sat is either bPAR i or bPAR d , bPAR insitu is the corresponding in situ measured benthic PAR, and n is the number of observations.The resulting error metrics are in multiplicative space and dimensionless.As such, multiplicative bias values of 1.0 indicates no bias, and bias <1.0 denotes negative bias (e.g., underestimation of the observed values).
Whilst for a multiplicative MAE, lower values (or closer to 1.0) would indicate better model performance.

Evaluation of model components: the diffuse attenuation model coefficients and spectral IOP-based approach (Kd(λ) vs. Kd(PAR))
It is worth noting that the K d (λ) model coefficients (0.005, 4.18, 0.52, and −10.18) in Eq. ( 6) were derived empirically from a dataset simulated with radiative transfer code [37].These coefficients best estimate layer-averaged K d (λ) for the photic zone depth (i.e. the layer of water over which irradiance is diminished to 10% of surface incident values).Because we wish to estimate irradiance at the seafloor, the Lee, et al. [37] K d (λ) model may be limited when the geometric depth is not equivalent to the photic depth.Hence, to evaluate the limitations of the K d (λ) model, we conducted a cursory comparison study for two optical scenarios: (i) 'oligotrophic' (chlorophyll-a concentration of 0.05 mg m −3 ) and (ii) 'mesotrophic' (chlorophyll-a concentration of 0.5 mg m −3 ); where the IOPs were known.We then compared layer-averaged K d (λ) values computed using the Lee, et al. [37] model against the K d (λ) values more accurately simulated with Hydrolight-Ecolight 5.1 (HE5) radiative transfer code [43].We used absolute percent difference (APD) as metric for this comparison.
We next evaluated the skill of our spectrally-resolved bPAR i model and two other approaches based on the method of Morel, et al. [41].These alternative methods use K d (PAR).In one model, K d (PAR) is modelled as a function of known/derived chlorophyll-a pigment concentration.We refer to this as the "Chlorophyll-based" approach.In the other approach, K d (PAR) is modelled as a function of known/derived K d (490).We refer to this as the "K d (490)-based" approach.Our spectrally-resolved bPAR i model is referred to as the "IOP-based" approach.In this case study we also used HES radiative transfer modelling to accurately simulate bPAR i values to which we could compare modelled bPAR i .For this radiative transfer modelling exercise, a large set of synthesized IOPs were used as HE5 inputs and simulations were run over 14 depths (spanning 3 -30 m), with fixed θ solz = 30° and a sand seafloor albedo.These simulated data were generated by McKinna, et al. [28] and further details can be found therein.Also, note that we used a consistent value of t g for transmitting surface irradiance values, either E s (λ,0 + ) or PAR, across the air-water interface.We note that for this radiative transfer study, model inputs for each of the three approaches (Chlorophyll-based, K d (490)-based, and IOP-based) were taken directly from the radiative transfer simulation outputs (i.e., no satellite algorithms were employed in this brief case study).

Temporal evaluation
Performance metrics described above (see section 3.4.2) are robust indicators of the absolute accuracy of the model but do not provide additional information on temporal structure or variability of the parameter being considered against the validation dataset used.To assess whether the benthic irradiance model produces temporally stable estimates of benthic PAR (e.g., no unreasonable extreme values or unrealistic trends in the time series), we examined a regional box, the Burdekin region, within the central GBR that encompasses three of the four test sites considered in the study (inset map on Fig. 1) by implementing the model to MODISA data from July 2002 to December 2018.

Simple depth sensitivity analysis
So far, we have used the MSL Beaman 3D-GBR bathymetry [29] as a proxy for watercolumn depth.However, one of the possible sources of error of our benthic irradiance model could be driven by changes in water-column depth due to tidal fluctuations.The tides in the GBR are predominantly semidiurnal (tidal cycle with two high and two low tides of approximately equal size each lunar day) and with regional differences in maximum amplitudes ranging between 2.5 to 7 m [55].To assess the importance of correcting the water depth for tidal effects, we also conducted a bPAR model run using in situ depth (i.e., pressure data measured using seafloor-mounted SEABIRD SBE16 + V2 SEACAT CTD with pressure sensor at the Yongala site) as z in Eq. ( 1).The depth sensitivity analysis (DSA), thus compared two bPAR model parameterizations: (i) the original model using MSL Beaman 3D-GBR bathymetry referred to as "3D-GBR" model, and (ii) a model using in situ CTD depth measurements referred to as "pressure" model, to assess and evaluate differences, if any.As neither of these two model parameterizations was assumed initially superior, difference was quantified using unbiased percent difference (UPD) calculated as per [50].

Spectral IOP-based K d model structure
Figures 2 and 3 show the various layer-averaged K d (λ) spectra calculated for four different layer depths and two optical water scenarios.The solid lines in the upper panels of Figs. 2 and 3 respectively show the radiative transfer-modelled layer-averaged K d (λ) values for oligotrophic and mesotrophic optical scenarios, while the black dotted lines (also in upper panels of Figs. 2 and 3) show the model layer-averaged K d (λ) calculated via Eq.( 6).For both examples, layer-averaged K d (λ) values for a given θ solz are generally consistent with each other over the spectral range 400 -600 nm, with the exception of the 30 m layer in the mesotrophic scenario.It can also be seen that as θ solz increases, the magnitude of K d (λ) also increases slightly whilst the shapes stay much the same.Plots of APD show that the analytical K d (λ) model typically agrees to within 10% of radiative transfer values over 400 -600 nm.We note that for the 30 m depth mesotrophic scenario, radiative transfer-simulated K d (λ) consistently deviates from the model of Lee, et al. [37] in both the blue (400 -450 nm) and red regions (600 -700 nm) for all solar zenith values.
This brief analysis, whilst not exhaustive, nonetheless provides some confidence that the Lee, et al. [37] model, although tuned for the photic depth layer, can estimate K d (λ) to within 10% of true values over a range of geometric depths.Indeed, Lee, et al. [37] model for K d (λ) has been previously tested in optically shallow coral reef waters of the West Florida Shelf by Barnes, et al. [50].Using in situ validation data, Barnes, et al. [50] showed the Lee, et al. [37] model had good skill at estimating K d (λ) at wavelengths shorter than 500 nm, which is consistent with this cursory case study.

Instantaneous benthic PAR matchup analysis
We acquired 1134 matchup pairs (i.e., days of observations) between satellite-derived bPAR i and in situ measured data across the test sites when only the L2 quality flags were used as exclusion criterion (Table 2).After the application of the two additional exclusion criteria (i.e., adjacent to a masked pixel in the 3 x 3 pixel box and CV > 0.15), 696 data pairs remained.These quality-controlled data pairs comprised the dataset that was subsequently used in the regression analysis.The highest number of valid concurrent data pairs were found for Yongala and Heron, with 348 and 210 days of observations considered for matchups, respectively.Palm Passage had the lowest number of valid concurrent data pairs at 45 days, as this test site had only one deployment period with good data because of in situ instrument failures experienced during the other deployment periods.Consequently, the number of concurrent matchups between the test sites varied by almost a factor of eight between sites with the most and the least data pairs.Nevertheless, across all test sites, only high-quality concurrent data pairs were included in the validation analyses, hence, the veracity of the model performance evaluation should not be diminished.
In situ measurements of instantaneous benthic PAR varied by almost two orders of magnitude ranging from 13 to 906 μmol m −2 s −1 , indicating that the validation dataset spanned a wide range of environmental conditions (Table 2).The satellite-derived bPAR i values showed comparable magnitudes spanning 22 to 672 μmol m −2 s −1 .Despite the wide range of inherent optical diversity across the test regions, regression analyses showed strong agreement between concurrent satellite-derived and in situ bPAR data pairs (Fig. 5, Table 2), indicating that our approach was able to produce realistic estimate of bPAR values across a range of complex and optically shallow coastal waters.The model performed strongest at Heron and was only slightly weaker at Yongala and the two deep oceanic test sites, Palm Passage and Myrmidon.Yongala and Heron had smaller dynamic ranges of benthic PAR over the period of observation compared to Palm Passage and Myrmidon, and may somewhat explain the better agreement found in the two former test sites.
Note that the capability of our newly developed benthic irradiance model may be limited under certain environmental conditions.For example, tropical cyclones can lead to extreme optical conditions via mixing and re-suspension of near bottom sediments which might not be fully captured by remotely sensed datasets.Excessive rainfall events during the wet season can also cause major flooding of river catchments [56].During such events, freshwater plumes discharge from rivers and disperse outward across the GBR lagoon, leading to distinct stratification in near-shore regions.Buoyant freshwater flood plumes are typically NAP-and CDOM-laden, thus the water column becomes optically inhomogeneous and the bPAR model's assumption of vertically homogenous IOPs is violated.

Table 2. Matchup statistics for instantaneous satellite-derived (bPAR i ) and concurrent in situ bPAR observations for the four test sites and all sites combined (denoted as ALL). The number of concurrent valid data pairs (days of observation) after all quality criteria
were applied and used in regression analyses is indicated as n.The number of valid data points when only the L2 flags were used as the exclusion criterion is given by (n 0 ), shown to emphasize the importance of additional exclusion criteria as described in section 3. 4 The bPAR model's limitation under optically extreme conditions can be seen at the Yongala test site during early April 2017 when several extremely low in situ bPAR values (filled triangles in Figs.5(a), 5(c), and 5(f)) resulted in a weaker overall agreement between the derived and actual measurements (i.e., lower slope 0.83 and higher intercept = 0.41, hence less unity, but high r 2 = 0.73) (see last row of Table 2).The in situ bPAR values, the minimum of which were an order of magnitude smaller than satellite-derived values during this extreme event period, coincided with Tropical Cyclone Debbie that affected the Queensland coast in late March 2017 [57] and led to low water clarity (due to high terrestrial river runoff) in the vicinity of the Yongala site long after the event [58].In Fig. 5(a), wet season (color-coded red) validation data for the Yongala site, which is seasonally dominated by the Buderkin River plume, were biased slightly high; again suggesting that the bPAR model is underperforming during optically challenging events.Although the scope of the current study did not allow to further explore this limitation, an additional investigation is warranted.
Further analysis conducted on the combined dataset also echoed the overall strong agreement between derived and in situ measured bPAR values (Figs.5(c) and 5(f), Table 2).Across all test sites, the calculated bias was small and ranged between 0.82 (-18%) at Palm Passage to 1.15 (+15%) at Yongala which respectively suggest an underestimation and overestimation of the observed PAR values by an average of ±16.5%.MAE values were very similar for all four test sites.Again, Heron had the smallest MAE of 1.22 while Palm Passage had the highest MAE of 1.31 indicating an average error of satellite-derived bPAR i ranging between 22% to 31%, respectively.The overall model error of the combined satellite-derived bPAR i values was about 25% (ALL dataset MAE = 1.25).This average model prediction error is much smaller and well within the 35% uncertainty and accuracy goal set for the retrievals of chlorophyll-a concentration which is a widely-used ocean color variable ( [59] as mentioned in [60,61]).Collective consideration of the bias and MAE suggest agreement between satellite-derived bPAR i in situ was best at the optically variable site Heron where many validation points existed, and lowest at the optically clear-water sites Palm Passage and Myrmidon, where relatively fewer data points were available.Given the overall strong agreement and small error of the derived values, we suggest that the larger error noted in Palm Passage and Myrmidon was more related to the quality of the in situ validation datasets and/or possibly related to fewer data points resulting in higher variability in the regression statistics rather than to a diminished ability of the model in these water bodies.

Daily integrated benthic PAR matchup analysis
The same in situ validation datasets used for the instantaneous matchup analysis were integrated to obtain daily in situ benthic PAR values and compared with concurrent satellitederived bPAR d .The dynamic range of the in situ and satellite-derived daily values were comparable at 0.4 to 23 mol photons m −2 d −1 and 0.6 to 17 mol photons m −2 d −1 , respectively.Results of the Type II linear regression analysis on the daily integrated datasets also showed an overall strong positive relationship (Fig. 6, Table 3).Again, agreement between derived and in situ values is most notable at the Heron test site with small negative bias (-4%) and highest accuracy (smallest MAE of 21%).Although Yongala showed a relatively lower coefficient of determination (r 2 = 0.65) amongst the test sites, the other regression statistics as well as bias and accuracy metrics indicate better model performance compared to the two deep test sites.As previously discussed, relatively weaker agreement between the satellite-derived and the in situ data at Yongala may be attributed to our bPAR model underperforming under more complex and/or stratified in-water optical conditions in this shallow inshore test site.More specifically, this effect may be compounded when estimating daily-integrated PAR due to intra-daily temporal variability in physical factors that drive nearshore IOPs.For example, clear water, offshore pixels will have optical variability that occurs at longer time-scales, being distant from terrestrial sources of particles and other pollutants, such that temporally extrapolating the IOPs derived at the time of the MODISA satellite overpass for the whole day, is a reasonable assumption and one that indeed showed strong agreement with observed values.However, for more dynamic inshore waters that may have IOPs with variability occurring at shorter time-scales or higher frequency (e.g., hours), it becomes more challenging to temporally extrapolate the IOPs over the whole day.On the other hand, while the agreement between satellite-derived and in situ values was strong for the deep oceanic test sites, the bias and MAE at Yongala showed better results than at Palm Passage or Myrmidon.Again, we suggest that these results may be related to the quality of the validation data obtained for these test sites.As for the combined dataset, satellite retrievals showed strong positive agreement with an almost negligible systematic error (bias of +1%) and high accuracy (low MAE of 25%) (see Figs. 6(c) and 6(f) and Table 3).Although we consider bPAR d to be an ecologically important parameter due to its role as an autotrophic requirement (e.g., energy source for photosynthesis), we also acknowledge that perhaps a more complete measure of light relevant for marine photosynthetic organisms is a PAR value estimated from the scalar irradiance [25,27] wherein light in every direction is considered (e.g., PAR measured by a spherical collector).However, as discussed earlier, the available validation dataset for this study is only collected for visible light measured on a plane (i.e., plane irradiance) which allowed us to only consider the downwelling plane irradiance in our model development.Indeed, if in situ scalar irradiance data becomes available in the future, this potential improvement could be considered.
While we acknowledge that our approach has several limitations and has room for further improvement, the baseline results obtained via this newly developed approach have demonstrated that our simple model can realistically estimate daily benthic PAR values within the GBR with good skill.Our new GBR-specific ocean color benthic irradiance algorithm can thus be utilized for ecological studies within this important and world-heritage listed region by providing estimates of benthic light values.

Time series case study
The 16+ year time-series of satellite-derived benthic PAR values for the Yongala site exhibits strong seasonality (Fig. 7).Lower values (purple to blue spots below the median lines) coincide at the end of austral wet season months (which extends from November to April the following calendar year, denoted as gray areas on both panels of Fig. 7) and beginning of austral dry season (May to October).Relatively higher values (green to red spots above the median lines) were derived for most of the austral dry season towards the proper austral wet season.This temporal pattern agrees with the variability of satellite-derived photic depth data previously obtained for GBR [21].This time-series provides insight into the influence of local events or phenomena that can drive variability of the underwater light field.As an example, decreased benthic light values are seen to coincide with tropical cyclones events (downward arrows in Fig. 7).Indeed, the future near-synoptic, high-temporal resolution output of this model, which will be the first of its kind for the GBR, will provide a critical dataset that will allow us to further explore and understand how light availability varies both in space and time and in relation to the presence of important benthic organisms that is found within the GBR.

Understanding the effects of tide on instantaneous benthic PAR estimates
Both "3D-GBR" and "pressure" models show good temporal coincidence with the observed instantaneous benthic PAR values (Figs.8(a) and 8(b)).For the purpose of this exercise, a large difference between the satellite-derived bPAR i values (for either model) would mean it is necessary to correct the original 3D-GBR bathymetry data for tidal effects.Otherwise, we can reasonably assume minimal tidal influence in the results of our new benthic irradiance model, but sentient of caveats when interpreting The comparison showed a mean UPD of 7.9% between the two models, equivalent to ~7 μmol photons m −2 s −1 (roughly 0.6 mol photons m −2 d −1 ) (Fig. 8(c)).Scatterplot of residuals between the two satellite-derived bPAR i values against in situ values (Fig. 8(d)) also show good agreement.This difference is relatively small compared to published ecologicallyrelevant minimum irradiance value required to maintain autotrophic health, which depending on the benthic organisms considered may vary between 0.4 to 5.1 mol photons m −2 d −1 [1].We concede that the tidal influence will vary depending on a location's known bathymetry and range of daily tidal fluctuations such that effects will be greater in areas with shallow depths and bigger tides.We note, however, that our model does not consider intertidal regions of the GBR or locations where MSL bathymetry is less than 5 m and where influence of the bottom can be significant [50].

Conclusion
Here we have demonstrated that our simple benthic irradiance model can be employed to estimate realistic benthic PAR values within the GBR with good accuracy.We emphasize once again that although our results were reported in terms of a single broadband PAR irradiance, our approach does not utilize K d (PAR).Instead, our model resolved K d (λ) by carrying out all calculations spectrally, where the full spectral variability of E b (λ) in Eq. ( 1) has been retained in the process.Thus, aside from its simplicity, the key strength and novelty to our benthic irradiance model mean it can be easily extended to other ocean color sensors.In particular, we expect that model-derived estimates of bPAR may be improved when applied to data collected by sensors with contiguous spectral bands such as the upcoming NASA's Plankton, Aerosol, Cloud, and ocean Ecosystem (PACE) mission which is scheduled to launch ca.2022.In addition, the bPAR model should also be extendable to other coastal waters for which high resolution bathymetry data and where spectral IOPs can be obtained accurately.
Overall, the validation exercise using four optically diverse test sites proved that our model was able to estimate both bPAR i and bPAR d values from MODISA remote sensing data with good skill as indicated by small bias and low MAE (within ~25%) compared to the observed "true" bPAR values.This suggests that our benthic irradiance model is robust to variation in optical water types found in the optically complex and shallow GBR lagoon.We do note that during extreme seasonal events, such as sediment-laden flood plumes, where acute optical stratification is present (i.e.buoyant turbid freshwater lens) the model's assumption of an optically homogenous water column is violated.Indeed, the validation results revealed that season retrievals of bPAR for the Yongala site, which is adjacent to the Burdekin River mouth, deviated from in situ observations.
We also acknowledge that the baseline bPAR model presented here has some limitations that if addressed, may further improve the model in the future.Firstly, our model assumed clear skies and that our validation data were generally obtained during cloud-free days, yet cloud cover variability modulates E s (λ) [25].This assumption proved to be satisfactory for instantaneous bPAR i estimates but may have significant consequences for deriving daily integrated bPAR d values; particularly in locations where temporal variation in cloud cover may be high.We concede that cloud variability should be considered when computing bPAR d , however, it is challenging to generate a cloud climatology dataset to correct our bPAR d data.Future enhancement may include statistically quantifying the effect that clouds have upon bPAR d by using in situ surface and sub-surface PAR datasets to empirically tune the model, both regionally and seasonally (e.g., as a cloud cover climatology).We note that parallel efforts are underway to include this as a future enhancement of our benthic irradiance model.Secondly, the atmospheric model used to obtain E s (λ) did not include diffuse sky irradiance.Future work to improve the bPAR model may include using a two component "clear sky" and "cloud layer" model similar to Frouin and McPherson [16].Thirdly, contributions of a wind-roughened sea surface to reflectance/transmittance across the interface may also be accounted for following Haltrin, et al. [62,63], provided appropriate ancillary wind data that matches our model input datasets can be obtained.Lastly, the applicability of our model may also be limited in areas where tidal fluctuations are much higher than what our DSA has tried to address.The bottom topography and physical processes that drive the oceanography of the entire GBR are complex and therefore cannot be easily generalized.Fine tuning the water column depth in shallow locations or in regions where tides vary significantly (e.g., Broad Sound located at the southern GBR) may also improve benthic irradiance estimates for these locations.
Despite the above limitations, the model results we have presented are very promising and can be implemented to the entire GBR region using near-daily MODISA data collected over the last 16+ years.This will allow for a spatiotemporally rich benthic PAR dataset to be generated for the first time at this scale.Given the importance of light availability for the health of benthic ecosystems, our model results can now be used to explore the spatial and temporal variability of benthic light in the GBR region (see Fig. 9).The resulting benthic light datasets will provide information needed to assess habitat quality for corals and seagrasses.Specifically, GBR-wide maps of light thresholds (i.e., minimum light requirements for optimal system function) for these organisms can be developed from the benthic PAR values derived using this model and form the basis of a new GBR water quality index based on light.Whilst we have outlined some limitations of the bPAR model, it is important to recognize that we do not necessarily anticipate that the satellite-derived benthic irradiance data presented here will become a standalone solution for monitoring water quality and ecosystem health in the GBR.Instead, we expect the bPAR model will complement existing monitoring tools such as the eReefs hydrodynamic model [64], the Q-IMOS mooring array, and the routine in situ water quality sampling programs.There is great potential for the uptake of the bPAR model in developing additional management criteria for monitoring water quality in the GBR based on light availability, as well as exploring effects of long-term patterns of benthic light availability on corals and/or seagrasses which has not been possible before now.
the results of our model.We also wish to thank Tommy Owens (NASA Ocean Biology Processing Group, OBPG) for initial processing of MODIS time-series data for Yongala, Chris Bartlett and associates (AIMS -IMOS) for the deployment and maintenance of the PAR loggers as well as the provision of the in situ datasets after each deployment, Patricia Menendez (AIMS) for the assistance in merging the satellite-derived and in situ datasets for the matchups, and the five anonymous reviewers for their valuable comments that helped improve this manuscript.

Fig. 1 .
Fig. 1.Location map of the four validation test sites (black filled circles) within the GBR region along the north eastern coast of Australia (inset map).The color gradient indicates depth contours within the 200 m bathymetric shelf.The blue rectangle indicates the boundary of the small regional box, Burdekin region, used for temporal evaluation of the model (as detailed in section 3.6) with the corresponding subset bathymetry showing the complex topographic features in the model region.Red filled circle indicates location of the Burdekin River mouth while gray masked regions indicate land and coral reefs.

Fig. 2 .
Fig. 2. Top row: solid lines represent layer-averaged spectral diffuse attenuation coefficients for four different layer depths for an oligotrophic scenario.Blue, orange, green, and red lines correspond to layer depths of 5, 10, 20, and 30 m, respectively.The dashed line represents layer averaged spectral diffuse attenuation coefficient derived by the Lee, et al. [37] model.Each panel (left-to-right) corresponds to four solar zenith angles (10, 30, 50, and 80°).Bottom row: absolute percent difference (APD) between radiative transfer modelled layer-averaged spectral diffuse attenuation coefficients and values estimated using the Lee, et al. [37] model.The horizontal gray line represents 1% APD.

Fig. 3 .
Fig. 3. Top row: solid lines layer-averaged spectral diffuse attenuation coefficients for four different layer depths for a mesotrophic scenario.Blue, orange, green, and red lines correspond to layer depths of 5, 10, 20, and 30 m, respectively.The dashed line represents layer averaged spectral diffuse attenuation coefficient derived by the Lee, et al. [37] model.Each panel (left-to-right) corresponds to four solar zenith angles (10, 30, 50, and 80°).Bottom row: absolute percent difference (APD) between radiative transfer modelled layer-averaged spectral diffuse attenuation coefficients and values estimated using the Lee, et al. [37] model.The horizontal gray line represents 1% APD.

Figure 4
Figure 4 shows scatter plots comparing modelled bPAR i with values computed from radiative transfer simulations.The results show that the IOP-based model derives bPAR i with root mean square error (RMSE) and mean bias statistics that are better than the Chlorophyllbased and K d (490)-based methods.Thus, this brief analysis demonstrates the benefit of using our spectrally resolved K d (λ)-based model as opposed to a K d (PAR)-based methods when estimating bPAR i .

Fig. 4 .
Fig. 4. Scatter plots that compare known bPAR i (i.e., derived from radiative transfer simulations) with modelled values.Each panel (left-to-right) corresponds to a different model approach: (a) Chlorophyll-based, (b) K d (490)-based, and (c) the IOP-based.The one-to-one line is plotted in black solid line.Data points are color coded by geometric depth.Highest values of bPAR i typically occur in shallow waters less than 5 m in depth.For our radiative transfer simulations, incident surface PAR values were 1679 μmol photons m −2 s −1 .

Fig. 5 .
Fig. 5. Scatterplots of concurrent log-transformed instantaneous satellite-derived (bPAR i ) and in situ bPAR for the four test regions of varying optical properties (a-b, d-e) plotted according to month of observation, and ALL sites combined plotted (c) by month of observation and (f) according to site.Color scale gradient for the months of observations are defined to delineate seasonal contrast between austral wet (November to April of following calendar year, yellowred) and austral dry (May to October, green) seasons.The thin and thick black solid lines indicate the 1:1 line and the reduced major axis regression slope, respectively.Filled triangles in (a, c, and f) indicate the extremely low bPAR values coinciding with a severe tropical cyclone but were excluded in the regression analysis.

Fig. 6 .
Fig. 6.Scatterplots of concurrent log-transformed daily integrated satellite-derived (bPAR d ) and in situ bPAR for the four test regions varying optical properties (a-b, d-e) plotted according to month of observation, and ALL sites combined plotted (c) by month of observation and (f) according to site.Color scale gradient for the months of observations are defined to delineate seasonal contrast between austral wet (November to April of following calendar year, yellow-red) and austral dry (May to October, green) seasons.The thin and thick black solid lines indicate the 1:1 line and the reduced major axis regression slope, respectively.

Fig. 7 .
Fig. 7. Timeseries plot of satellite-derived (a) instantaneous benthic PAR, bPAR i , and (b) daily benthic PAR, bPAR d , for the Yongala test region from July 2002 to December 2018.The color gradient indicate bPAR values while dashed horizontal line indicate the median value of the entire time series data for each parameter, overlaid for reference.The grayed out areas indicate the austral wet season months (extending from November to April of following calendar year).Downward pointed black arrows on the x-axes denote occurrences of selected severe tropical cyclones that hit the eastern Australian coast.Note that vertical axes are in logarithmic scale.

Fig. 8 .
Fig. 8. Depth sensitivity analysis and model results comparison at the Yongala test site.(a) plot of satellite-derived instantaneous benthic PAR values modeled using Beaman 3D-GBR bathymetry as z (red, 3D-GBR model) with in situ instantaneous benthic PAR (black); (b) plot of satellite-derived instantaneous benthic PAR values modeled using in situ pressure data as z (blue, pressure model) with in situ instantaneous benthic PAR (black); (c) unbiased percent difference (UPD) between bPAR i , calculated as: |(bPAR 3D-GBR -bPAR pressure )|/(0.5 * bPAR 3D-GBR + 0.5 * bPAR pressure ) * 100; and (d) scatterplot of the residuals (against in situ values) of the two models where diagonal black solid line denotes the 1:1 line.Gaps in plots a-c indicate periods where there were no in situ data available for validation.

Fig. 9 .
Fig. 9. Map of 2016 annual mean of daily integrated benthic PAR (bPAR d ) for Great Barrier Reef region.Color gradient indicates values in mol photons m −2 d −1 .