Extracting Remotely Sensed Water Quality Parameters from Shallow Intertidal Estuaries

Abstract

Sentinel-2 imagery is potentially ideal for providing a rapid assessment of the ecological condition of estuarine water due to its high temporal and spatial resolution and coverage. However, for optically shallow waters, the problem of isolating the effect of seabed reflectance from the influence of water properties makes it difficult to use the observed surface reflectance to monitor water quality. In this study, we adopt a methodology based on Lyzenga’s model to estimate water quality properties such as the dominant wavelength and diffuse attenuation coefficient (K_d) of shallow estuarine waters. Lyzenga models the observed reflectance (R) using four parameters: total water depth (z), sea-bed reflectance (R_b), water reflectance (R_w) and K_d. If R_b is known a priori and multiple observations of R are available from different total water depths, we show that Lyzenga’s model can be used to estimate the values of the remaining two parameters, K_d and R_w. Observations of R from different water depths can either be taken from the same image at different proximal locations in the estuary (“spatial method”) or from the same pixel observed at different tidal stages (“temporal method”), both assuming homogeneous seabed and water reflectance properties. Tests in our case study estuary show that K_d and R_w can be estimated at water depths less than 6.4 m. We also show that the proximity restriction for the reflectance correction with the temporal method limits outcomes to monthly or seasonal resolution, and the correction with the spatial method performs best at a spatial resolution of 60 m. The K_d extracted from the blue band correlates well with the observed K_d for photosynthetically active radiation (PAR) (r² = 0.66) (although the relationship is likely to be estuary-specific). The methodology provides a foundation for future work assessing rates of primary production in shallow estuaries on large scales.

1. Introduction

Human effects on estuarine water quality include increased levels of nutrients and suspended solids and reduced light penetration, all with ecological flow-on effects. The optical properties of estuarine water determine light availability and can control both submerged benthic and pelagic primary production by shortening the depth of the photic zone in turbid water. It can also alter the vertical heat stratification, influencing nutrient availability and the photosynthetic activity of aquatic plants [1,2]. The impacts of reduced microphytobenthos (MPB) and macrophyte (e.g., seagrass) productivity are diminished food supply and nursery habitat, as well as changes to ecological interaction networks, which may trigger tipping points in ecological states and contribute to eutrophication cycles [3,4,5,6]. Considering the relationship between MPB/seagrass and their surrounding aquatic environment, there is a clear need to monitor the optical properties and turbidity levels of estuarine water at scale to inform sound management practices that will prevent potential habitat loss.

To monitor water optical properties at large temporal and spatial scales, remote sensing is widely used to assess water column turbidity, coloured dissolved organic (CDOM) and chlorophyll-a (chl-a) concentrations by building empirical or semi-empirical regression models between remotely sensed colour bands and in situ data [7,8]. Recent work has also used water colour evaluated by dominant wavelength as an intuitive way of monitoring water quality (e.g., Lehmann et al. [9]), since it can reflect the water colour as perceived by the human eye in a numerical way. Woerd and Wernand [10] proposed empirical algorithms to estimate the hue (the dominant wavelength) and intensity of light using a weighted linear sum of remote sensing reflectance blue, green and red bands, which made it possible to transform satellite signals into colour science data [11]. These algorithms have been extended to other sensors deployed on Landsat and Sentinel-2 [12]. However, there are few studies using satellite-derived watercolour as a compound indicator of water quality in shallow coastal waters due to the difficulty of obtaining true water reflectance.

The accuracy of using remotely sensed observations for water quality and water colour monitoring depends on the methods of obtaining the true remote sensing reflectance from observations. In optically shallow water, the seafloor contributes significantly to the water-leaving signals because photons reflected by the seabed contribute to the measured radiance [13,14]. This reflected radiance is a function of the seabed properties (e.g., presence of vegetation, grain size) and causes biases (i.e., a non-water column signal) in remote sensors from which estimates of water quality are derived. Internally reflected light can be ignored because it is not the dominant source of additional upwelling radiance [15]. Seabed reflectance can be removed in intertidal regions if the bottom types are known, and light is assumed to attenuate exponentially with depth during submergence [16]. The seabed reflectance in visible bands can be estimated by the particle size when the sediments are exposed [17]. Based on Lyzenga’s theories, Maritorena et al. [18] proposed an equation to demonstrate the relationship between observed reflectance, bottom reflectance and deep-water reflectance, in which the bottom is assumed to be a Lambertian object:

$$R=R_w\left(1-e^{-2k_{d}z}\right)+R_b{e}^{-2k_{d}z}=\left(R_b-R_w\right)e^{-2K_{d}z}+R_w$$

(1)

where R is the reflectance of a visible band observed by sensors just above the water surface (spectral dependency is omitted for simplicity of notation); R_b is seabed reflectance; R_w is water reflectance at infinite depth; z is water depth and K_d is (the band-specific) diffuse attenuation coefficient. The spectral analysis software Hydrolight provides several theoretical bed reflectance characteristic models, which can be used in addition to in situ water quality measurements and Lyzenga’s model to extract K_d and R_w from imagery [19]. Recent extensions to the Hydrolight capability include wind speed and zenith angle in methods to estimate the true water surface reflectance [20].

The primary objective of this study was to develop a methodology based only on Lyzenga’s model (without the need for in situ water quality measurements) to remove the seabed reflectance and isolate the water reflectance and diffuse attenuation coefficient for optically shallow estuarine water. The corrected water reflectance could then be used to extract true watercolour (dominant wavelength) and to evaluate the relevance of this information (as well as K_d) to monitoring estuarine ecological state. We also show how an established regression between the median sediment particle size and reflectance could be used to extend the corrections to shallow subtidal waters.

2. Materials and Methods

2.1. Study Area and Required Data

The study was conducted in Tauranga Harbour (37°39′S, 176°11′E, New Zealand), including its many sub-estuaries (Figure 1a). The harbour is large (242 km²), predominantly shallow (60% intertidal), with a semi-diurnal spring and neap tidal range of 1.62 m and 1.24 m, respectively [21]. Seagrass (Zostera muelleri) only occurs in the intertidal regions (22.5% of the emerged area), with the remaining seafloor consisting of unconsolidated soft sediments. The seagrass biomass has a seasonal fluctuation where it gradually decreases in winter and reaches a minimum value in early spring [22]. Fine sediment accumulates in the upper reaches of the estuary, and the main freshwater source of sediment is the Wairoa River [23]. Therefore, these two factors have long-term and seasonal effects on water colour signatures and influence seabed characteristics.

This study used the Sentinel-2 Level-2A data product downloaded from the European Space Agency. In this product, satellite-measured top-of-atmosphere radiances have been corrected for atmospheric effects to yield surface reflectance. We used the bands (blue, green and red) at their native 10m pixel resolution. Fifty-one Level-2A images with cloud coverage less than 5% were selected for the Tauranga Harbour region. These images were divided into the northern and southern harbour at the natural sand barrier (dashed white line, Figure 1) to compare the difference in water optical quality between basins.

Water depth at each pixel at the time of data acquisition was derived from a merged elevation-bathymetry map and tide levels. The merged map was created by combining a digital elevation model (collected by airborne light detection and ranging (LiDAR)) of the intertidal regions with subtidal bathymetric data (single and multi-beam) assembled by the Port of Tauranga. The bathymetry has an accuracy of ±0.13 m [24]. Sediment grain size observations (Figure 1) were made into a spatial map with inverse distance weighting by Rullens et al. [25] (Figure 1b–d). Photosynthetically-active radiation (PAR) measurements used in Flowers et al. (in prep.) were also collected at the sites shown in Figure 1a.

Figure 1. Optically deep-water regions (where no corrections were needed) and in situ measurement (sediment and PAR) locations in Tauranga Harbour. (a) The in situ sediment data in the intertidal (for regression) and subtidal (for extrapolation) were provided by Clark et al. [26] (and Bay of Plenty Council) and Ellis et al. [27], respectively (see Methods). The PAR data were provided by Flowers et al. (in prep.). The location of the spatial transect shown in the figure in Section 2.2.1 is marked with a black circle. The black dashed line is the natural barrier that isolates Tauranga Harbour into the northern and southern harbours. The selected transect was close to the northern entrance of the harbour. (b–d) Maps of sediment content percentages (gravel, sand and mud) from Rullens et al. [25]. The background image is from Sentinel-2.

Figure 1. Optically deep-water regions (where no corrections were needed) and in situ measurement (sediment and PAR) locations in Tauranga Harbour. (a) The in situ sediment data in the intertidal (for regression) and subtidal (for extrapolation) were provided by Clark et al. [26] (and Bay of Plenty Council) and Ellis et al. [27], respectively (see Methods). The PAR data were provided by Flowers et al. (in prep.). The location of the spatial transect shown in the figure in Section 2.2.1 is marked with a black circle. The black dashed line is the natural barrier that isolates Tauranga Harbour into the northern and southern harbours. The selected transect was close to the northern entrance of the harbour. (b–d) Maps of sediment content percentages (gravel, sand and mud) from Rullens et al. [25]. The background image is from Sentinel-2.

2.2. Methods

Equation (1) shows that the reflectance measured at the water surface contains a contribution of R_b, which depends on K_d and z. To solve Equation (1), R_w, R_b and K_d must be approximated. For R_b, we either measured it directly when the sediments were exposed at low tide (intermittently exposed, Figure 2), or we used the relationship between R_b and particle size to infer R_b in permanently inundated regions (Figure 2) (details in Section 2.2.1). With R_b constrained, we then optimised Equation (1) using observed reflectance collected over multiple water depths, assuming that R_b, R_w, and K_d are constants over these water depths (detail in Section 2.2.2). In the deep channels, z was large enough relative to K_d that R~R_w (deep water in Figure 2).

The spectral information of intertidal vegetation (generally seagrass) in the harbour is different from the soft sediments, which change R_b. Intertidal areas where seagrass is abundant can be mapped using a supervised classification scheme tested by Ha et al. [22]. For simplicity, here we exclude these areas and use non-seagrass regions to extract R_w and K_d. However, seagrass areas could also be used by measuring the R_b of these regions while exposed (Tauranga Harbour has almost no subtidal seagrass).

2.2.1. Estimation of Seabed Reflectance

Our method for estimating R_b depended on whether the region was inundated, and so each pixel in the Tauranga Harbour image collection was classified in a simple way: exposed, intermittently inundated and perpetually inundated (Figure 2). Pixels were classified as exposed or inundated using the normalised deviation water index (NDWI) calculated with Band 3 (Green) and Band 8 (NIR) from the Sentinel-2 data [28]. The variation of the NDWI has been shown to be good for detection of waterbodies or changes in water level, using 0.3 as a threshold [29].

$$\mathrm{NDWI}=\frac{R\left(Green\right)-R\left(NIR\right)}{R\left(Green\right)+R\left(NIR\right)}$$

(2)

For the intermittently exposed sediments, R_b could be derived from images taken at low tide. The curve in Figure 3a shows a synthetic example of R calculated using Equation (1) as a function of water depth ranging from 0 m (exposed) to 1.2 m (inundated). In the exposed region, the depth was set to zero and the reflectance was constant (the red rectangle in Figure 3a), and in the inundated region, the reflectance decreased exponentially to deep water, with the rate of decay dependent on K_d. Note that each visible band (blue, green and red) has a different decay rate and a different R_b. The observed reflectance confirmed the expected theoretical trend (examples shown in Figure 3b); the reflectance of the exposed pixels varied only slightly due to minor changes in the bed substrate, and R decreased asymptotically with water depth.

For the perpetually inundated regions, R_b needs to be estimated to correct R (Figure 2) because it is not possible to obtain direct observations of R_b. To do this, we trialled three relationships (linear, logarithmic and exponential) between the in situ measured median particle size (d) and the reflectance of exposed pixels in each visible band (following Sadeghi et al. [17]), which we then applied to predict R_b using grain size maps for subtidal regions. Although R_b may depend on mineralogical composition as well as particle size, in Tauranga Harbour, the sediments in deeper water are from the same source as those in shallow water, which should restrict the variations in grain composition. The in situ median particle size data (115 samples) were derived from an ecological survey of Tauranga Harbour [27] and monthly measurements collected by the Bay of Plenty Regional Council in the intertidal regions (Figure 1a). Combined with 45 samples from subtidal regions collected by Clark et al. [26], Rullens et al. [25] used these in situ data to create sediment texture maps by using inverse distance weighting interpolation (Figure 1b–d). These sediment texture maps were resampled from 100 m to 10 m (Figure 1b–d) and were used to estimate the median particle size of each pixel in the perpetually inundated regions. In summary, spatial maps of R_b were constructed using direct measurements of seabed reflectance in regions that were exposed at low tide and using inferred values in subtidal regions (inferred using the intertidal relationships between median particle size and reflectance (Figure 2)).

2.2.2. Calculation of R_w and K_d

Having established a protocol to determine R_b across the range of elevations in the harbour, the next step was to solve the optimization problem and extract an estimate of R_w from an observed signal R using Equation (3) (Equation (1) rearranged). A numerical search algorithm was applied to seek the value of R_w that minimised the variance (V) and trend (S) in K_d(i) over a range of depths z_i (assuming that the water properties are homogeneous across these pixels and any systematic differences in K_d at different depths were caused by the incorrect choice of R_w).

$$K_{d\left(i\right)}\left(\lambda \right)=\frac{1}{2z_i}\times \ln \left(\frac{R_b\left(\lambda \right)-R_{w,}\left(\lambda \right)}{R_i\left(\lambda \right)-R_w\left(\lambda \right)}\right)$$

(3)

where R_i is the reflectance at pixel i with water depth z_i in band λ. To be assumed to be homogeneous, the pixels needed to be in close proximity. From here on, we drop the notation λ for clarity, but all analyses were performed separately for each visible band.

Given that proximal pixels (either in space or time) would have different depths in a shallow estuary, either because the depth changed gently over the intertidal region (spatial changes) or because the tide changed through time (temporal changes), the requirement of evaluating Equation (3) at multiple water depths to solve the optimization problem could be satisfied. An example of how S (the slope) and V (the variance) of K_d(i) change over water depth is provided in Figure 4, where each line represented a different “guess” of the value of R_w. The correct R_w in this example was 0.031 because it provided the same K_d at all depths while minimising V. Once R_w is known, the corresponding K_d value can also be calculated by using Equation (3).

Proximal pixels with varying water depths could either be selected from a single image where the water depth varied spatially due to the gently sloping intertidal bathymetry or could be selected from a series of images collected at the same location but at different tides. Hereafter, we called these two approaches the “spatial method” R_ws and “temporal method” R_wt, respectively. In the first case, R_ws retained the original temporal resolution, but the spatial resolution was reduced. To retain cross-shore resolution across the intertidal, we implemented the spatial method using square tiles of pixels with different water depths rather than transects. Therefore, the size of the tiles determined the spatial resolution of R_ws. In the second case (R_wt, reflectance corrected by the temporal method), the spatial resolution was the same as the original images (10 m), but the temporal resolution was sacrificed depending on the time period over which the images were selected and associated with different tides.

In order for the methods to work, the underlying water properties and bed properties should be uniform over the region of proximal pixels. Therefore, when we applied the spatial method, we conducted a manual check on the selected pixels to avoid areas where seagrass and sediment varied substantially or where seagrass and bare sediment occurred together within the same tile. Considering the spectral information of seabed coverage did not vary substantially over short time frames, we expect variation in vegetation (seagrass) will not have a significant effect on the temporal method.

2.2.3. Assessing the Optimal Resolution to Use in Spatial and Temporal Corrections

There is a trade-off between the number of pixels that are used to evaluate R_w and the requirements for homogeneity. More input data should allow Equation (3) to be optimised with greater confidence but may increase the possibility of non-uniform reflectance, in which case it is difficult to detect a minimum in the variance and slope of Equation (3). In total, we tested the spatial method with tile sizes of 40, 60, 80 and 100 m and the temporal method with images from the same month, season and year. Where values were missing due to the algorithm not finding a minimum within a likely range, nearest neighbour interpolation was used to fill gaps in the maps. We assumed that the variation caused by errors in the method should be reduced across the region when the optimal resolution is used for processing. The resolution should maximise the information used to calculate the correction at each location but should still be smaller than large-scale changes in the environment (such as seasonal changes to estuarine characteristics caused by, e.g., changes to seagrass areas, and spatial changes to estuarine characteristics, such as those caused by, e.g., changes to sediment grain sizes). The standard deviation (STD) of corrected reflectance (R_w) in each visible region (using either the spatial or temporal method) from shallow water regions was used as a measure to evaluate the model. In addition, a successful correction should return small missing value rates (n), where the missing value rate is the percentage of pixels where Equation (3) could not be solved (where no minima are found). Therefore, STD and n were assessed to evaluate the best temporal and spatial resolution.

2.2.4. Validating K_d

Due to the difficulty of obtaining direct measurements of K_d at a specific wavelength, we instead used the in situ measured PAR values (E_d) at the same locations but at two different water depths (z) to calculate the K_d (PAR) using Equation (4). We then compared the K_d values in visible bands to K_d (PAR) to verify the model for ecological use in shallow water regions.

$$K_d\left(\mathrm{PAR}\right)=\frac{1}{z_2-z_1}ln\left[E_d(z_2\right)/E_d\left(z_1\right)]),$$

(4)

where E_d is downward irradiance over the wavelengths of PAR.

2.2.5. Calculation of Dominant Wavelength Using R_w

The dominant wavelength, which reflects the natural colour of objects, is the wavelength of pure spectral colour derived as the intersection of the line joining the white point to the point of interest (in the form of chromaticity coordinate pairs (CIE x, CIE y)) at the border of the Commission Internationale de l’Eclairage (CIE) horseshoe chart [11], shown in Figure 5.

The chromaticity coordinates (CIE x, CIE y) are calculated by normalising tristimulus values (X, Y, Z) (Equations (5) and (6)), which can be obtained by using the multiple linear regression relationships developed in Woerd and Wernand [10] (Equations (7)–(9)). The chromaticity coordinates provide a measure of watercolour from the reflectance of blue, green and red bands (separately corrected for the effect of the seabed).

$$x=\frac{X}{X+Y+Z}$$

(5)

$$y=\frac{Y}{X+Y+Z}$$

(6)

$$X=6.423R_w\left(\mathrm{Blue}\right)+53.696R_w\left(\mathrm{Green}\right)+32.028R_w\left(\mathrm{Red}\right)$$

(7)

$$Y=22.289R_w\left(\mathrm{Blue}\right)+65.702R_w\left(\mathrm{Green}\right)+16.808R_w\left(\mathrm{Red}\right)$$

(8)

$$Z=31.101R_w\left(\mathrm{Blue}\right)+1.778R_w\left(\mathrm{Green}\right)+0.015R_w\left(\mathrm{Red}\right)$$

(9)

3. Results

3.1. Estimation of R_b in Shallow Water

For the exposed regions, R_b was measured directly, but in permanently inundated regions, R_b was inferred. One hundred and sixteen in situ samples were used to develop a regression fit between exposed reflectance R_b and seabed particle size (21 samples were excluded due to the presence of seagrass and other vegetation or because the pixel at that location was inundated at the time that the reflectance image was collected). Of the three regression models trailed, the logarithmic regression provided the relationship with the highest r² values (

Table 1. The r² values and coefficients of the logarithmic regression ($R_b=a_1\log \left(d\right)+a_2$ for the blue, green and red bands from 2018 to 2020 (p < 0.05), where d is median particle size (μm).

Year	Band	a1	a2	r2	Standard Deviation of Rb
2018	Blue	−0.008	0.107	0.61	0.015
	Green	−0.010	0.120	0.66	0.018
	Red	−0.006	0.105	0.45	0.018
2019	Blue	−0.008	0.102	0.64	0.013
	Green	−0.007	0.112	0.53	0.016
	Red	−0.007	0.107	0.55	0.015
2020	Blue	−0.006	0.090	0.72	0.015
	Green	−0.008	0.119	0.61	0.020
	Red	−0.006	0.102	0.56	0.016

), explaining between 5% and 11% more variation than other models.

The seabed optical properties did not change greatly over three years, as indicated by the small changes in the standard deviation of reflectance in each band as well as minor inter-annual variation in the fitting coefficients (Table 1 and Figure 6). These regression relationships were applied to estimate R_b in the perpetually inundated regions (i.e., where R_b could not be observed directly). A spatially resolved example of R_b is shown in Figure 7. The R_b values in the visible bands remain similar in the perpetually inundated regions (where the grainsize is similar) and reach lower values (particularly in the red band) in the middle, intertidal regions of the harbour where seagrass is the dominant vegetation.

3.2. The Corrected Reflectance and Assessing the Optimal Resolution

3.2.1. Correction with Temporal Method

In the case of the temporal method (R_wt), four, nine and sixteen Sentinel-2 images taken at low tides with low cloud coverage were selected to be used in the monthly (February 2019), seasonal (Summer 2018–2019) and annual (2019) corrections, respectively. Using only a month’s worth of images meant that conditions varied less, but there were fewer degrees of freedom used to estimate reflectance. Nevertheless, extraction using the monthly and seasonal corrections had a greater success rate with a lower standard deviation (STD) in the optimal values of R_wt (in Figure 8) and n (

Table 2. The STD and n of the R_wt in visible bands after monthly, seasonal and annual correction, respectively.

Band	Evaluator	Monthly	Seasonally	Annually
Blue	STD	0.012	0.011	0.012
	n	17.98%	35.62%	54.06%
Green	STD	0.011	0.010	0.013
	n	16.28%	27.21%	50.18%
Red	STD	0.007	0.006	0.008
	n	15.15%	34.61%	52.69%

). The annually resolved correction performed the worst, with more than 50% of pixels returning missing values in each visible band and a high value of averaged STD after using temporal methods. Therefore, we selected monthly-resolved R_wt to calculate the dominant wavelength and K_d of shallow water in Tauranga Harbour.

Correcting for seabed reflectance using the temporal method changed the signals substantially in most regions, removing the patterns in R caused by the effect of R_b (the uncorrected images showed patterns associated with underwater channels (Figure 8)). Although in regions of large spatial variations (such as near the interface of sediment and water), the method was less effective, some original signatures of the seabed were still evident. The corrected R_wt in shallow water was reduced to a similar range as that of the deep channels (Figure 8), where no correction is needed, which also contributed to a reduction in STD. Sometimes large seagrass meadows in shallow water made it difficult to obtain R_wt because R_b varied too much. The pixels that were difficult to correct (where a solution to Equation (3) could not be found) were either close to a deep channel or surrounding irregular seabed features (missing values in Figure 8).

The dominant wavelength of shallow water decreased by ~20 nm after correction using the temporal (applied monthly) method, from yellow (reflecting the influence of the generally sandy bed) to a greener colour (compare Figure 9b,c). The corrected colour of shallow water was close to that of deeper channels near the entrance, where the bottom was assumed not to contribute to water reflectance (Figure 9d). The exposed sediments had fairly homogeneous optical properties with a dominant wavelength of ~570 nm (Figure 9a).

3.2.2. Correction with Spatial Method

The spatial method (R_ws) decreased the original spatial resolution (10 m) depending on the size of tile used to fit Equation (3). We tested the optimal number of pixels to use in fitting the Sentinel-2 image on 5 February 2019 at high tide. The optimum resolution of the spatial method was 60 m, where the STD of corrected reflectance was also minimal relative to other resolutions with a low missing data rate (

Table 3. The STD and n of the R_ws in visible bands using four different tile sizes with spatial method (40, 60, 80 and 100 m), respectively.

Band	Evaluator	40 m	60 m	80 m	100 m
Blue	STD	0.012	0.008	0.012	0.011
	n	16.35%	19.67%	20.18%	21.89%
Green	STD	0.014	0.013	0.011	0.012
	n	15.72%	17.31%	18.65%	23.54%
Red	STD	0.010	0.009	0.011	0.010
	n	17.19%	18.85%	20.72%	21.29%

). The 40-m or 80-m resolved corrections provide similar results and are useable when a higher or lower resolution of data are needed, respectively. The 100-m resolved correction was poorer, with n of 21.89%, 23.54% and 21.29% in each visible band (Table 3). Therefore, in this study, we used a 60-m resolved correction to remove the seabed reflectance and calculate the dominant wavelength and K_d.

Although the spatial method generally showed good performance in shallow water, the results might be missing or biased in the regions where the input reflectance of the seabed varies substantially. These regions included the interfaces between exposed sediments and shallow water, between shallow water and deep water and between vegetation and sediments or water. In shallow water, the corrected reflectance in each band was significantly reduced, by 13% (blue), 25% (green) and 23% (red), relative to the uncorrected one (Figure 10).

The area of the exposed seabed (reflectance values not needing correction) was small in this example (compared to Figure 9) because high tide covered the sediments (Figure 11a). The colour of shallow water regions after correction was significantly different from that before correction. In the uncorrected image (Figure 11b), the seabed morphology was clearly evident, with the shallow regions differentiated from the green channels; these patterns disappeared after correction (Figure 11c). When the corrected image was replotted using a finer resolution colour scale (Figure 11d), the intersection between channels and shallow water regions might show a slightly erroneously high value of the dominant wavelength. This might result from a drop in the number of uniform data points used in the processing to estimate the true water reflectance. The deep channel was not substantially affected by the seabed reflectance, and so it did not need correction (Figure 11d).

The effect of our bottom correction on dominant wavelengths can be visualised using chromaticity coordinate plots (Figure 12). The spatial and temporal methods showed consistent results, comparing the dominant wavelength before and after correction. The dominant wavelength ranged between 520 nm and 578 nm and between 540 and 580 nm for uncorrected data, and from 521 nm to 574 nm and from 572 nm to 540 nm after correction with the temporal and spatial methods, respectively (Figure 12). The scatter using the spatial method declined partially because the spatial method decreased the spatial resolution (Figure 12a), while results from the temporal method were not affected in the same way (Figure 12b). In both plots, some points collected in shallow water with a dominant wavelength over 570 nm before correction likely showed the effect of seabed reflectance on the true water reflectance, considering their colour was almost the same as exposed sediments. There were clusters of points at the dominant wavelength of 560 nm~570 nm before correction. After correction using both methods, the dominant wavelength of these pixels returned toward the greener region (<565 nm) and the majority of pixels from shallow water cluster around the dominant wavelength of between 540 and 560 nm. The lower limit of watercolour before and after correction is close because seabed reflectance did not need correcting in deeper channels (520 nm~540 nm). Some pixels from the shallow water region showed a similar colour to the deep water after correction, which meant the seabed reflectance was effectively removed.

3.3. The Diffuse Attenuation Coefficient (K_d)

The K_d values derived after application of both spatial and temporal correction methods to the visible bands shared similar features. Here we show K_d (Blue) as an example (Figure 13). The missing values mainly came from pixels in the deep channels (it was only possible to estimate K_d where the seabed can be detected, Equation (1)) and exposed sediments. K_d (Blue) usually had a minimum value (~0.12) in the pixels close to the deep channels, which reflected the probability that turbidity would be low and clarity high in areas that were more flushed with open ocean water. Conversely, the K_d (Blue) tended to be higher at the margins of shallow water and exposed sediments. These regions were usually close to the surrounding subestuaries and river mouths, where sediment was stirred into the water column by wave action, or delivered to the estuary by the river.

In order to assess the ecological relevance of the satellite-derived K_d, K_d in each visible band was then compared to the in situ measured K_d (PAR) (the measurement location is shown in Figure 1). Although the correlation between K_d from each visible band and K_d (PAR) is usually estuary-specific, depending on turbidity/chl-a concentrations, the fairly high correlation (Figure 14) between the in situ K_d (PAR) and K_d (Blue) might indicate that K_d (Blue) could be used as a proxy for K_d (PAR) for routine ecological monitoring over large regions.

3.4. Monthly Variations in Watercolour and K_d

Monthly variations of shallow and deep reflected water colour over a three-year period showed an apparent seasonal fluctuation, whereas by comparison, the exposed sediments had near-stable optical properties using a combination of the temporal (monthly) and spatial method (60 m) (when image number was fewer than four per month, we used the spatial method instead and averaged the results for that month) (Figure 15a,b). The dominant wavelength decreased to a minimum in winter and returned to the maximum in summer, which meant that the watercolour became greener in summer than in winter (the colour scale on the left of Figure 15a,b reflected true colour). The southern Harbour had slightly greener water colour in shallow and deep water regions (Figure 15b), but the overall water colour in Tauranga Harbour was still stable from 2018 to 2020.

We found that diffuse attenuation was mainly affected by location, suggesting that tidally driven variation in light-attenuating material was small relative to tidal conditions. K_d of the southern harbour was higher in the blue, green and red bands than that of the northern harbour (Figure 13c), meaning a generally lower-light state. The red and green bands of both harbours had the largest and smallest values of K_d, respectively. The K_d in the band red also had the greatest standard deviation in both harbours. The seasonal fluctuation of K_d (Blue) was small, but the values tended to increase in spring and autumn (Figure 15c). The tide conditions seemed to have a weak effect on K_d. The K_d value of ebbing water tended to be higher than that of flooding water (Figure 15d).

4. Discussion

In this study, we developed a new application of Lyzenga’s model to provide detailed measurements of watercolour and K_d from satellite reflectance measurements of shallow intertidal estuaries by removing the effect of seabed reflectance. The fraction of reflectance accounted for as the bottom signal by our method was the greatest in shallow water and decreased with water depth in a manner consistent with the law of exponential light attenuation in water (Figure 16). Note that the error depends on K_d (high K_d weakens the signal from the bottom at the surface), and the example shown here is for our case study. Compared to traditional geostatistical methods for creating maps of water quality (e.g., inverse distance weight interpolation applied to in situ measurements, which are more suitable when water quality does not vary spatially substantially), this application can provide a more accurate representation of a large region by accounting for substantial spatial variations in water quality. Therefore, the ecological information provided by this new measurement is useful in these settings to underpin management actions such as setting loading limits for estuaries and monitoring the general health of waters.

We propose that water colour corrected for bottom effects is a useful measure for water quality monitoring because it integrates the effect of suspended particles and coloured dissolved substances in the water column, including chl-a [30], CDOM [31] and suspended particular matters (SPM) [32]. The dominant wavelengths reached a peak (shift toward green) in summer (Figure 14a), which was consistent with the patterns observed in the in situ chl-a concentrations in Tauranga Harbour [33]. Though the concentrations of CDOM, SPM, salinity and suspended sediments varied seasonally [32], chl-a is usually the principal colourant in low-level turbidity estuaries [34]. We found averaged dominant wavelength values of shallow and deep water in the green part of the spectrum (528 nm and 545 nm, respectively) (Figure 15b), confirming the importance of phytoplankton pigments for optical water quality in Tauranga Harbour. Moreover, the dominant wavelength in the yellow to red part of the spectrum in the upper fringes of the inundated regions (Figure 9 and Figure 11) was consistent with CDOM and SPM loading from rivers and subestuaries as well as bottom sediment resuspension driven by waves [9,34,35]. Note that because of the requirement in our method of a set of proximal pixels with uniform reflectance, the bottom correction using the spatial method could not be applied at the edge of subtidal and intertidal regions where high variation of chl-a and CDOM concentrations occur [26]. These both cause significant local variations in bottom reflectance and the water depth correction.

Another ecologically relevant parameter that can be derived from satellites is the diffuse light attenuation rate, which governs the photosynthetic rate of inundated areas and thus the productivity of microphytobenthos and submerged seagrass. The variation of K_d responds to multiple drivers and reflects temperature, the input of sediment and the hydrodynamic regime, which the resolution of our remotely sensed data allow us to explore. Our results showed an increasing trend of K_d in spring, likely associated with phytoplankton blooms triggered by temperature increases, which usually bring elevated chl-a concentrations depending on the cell size and pigments [36] (Figure 13). For the remaining seasons, SPM and CDOM generally dominated the turbidity, which would also affect K_d [37]. In autumn, K_d remained similar to the spring value because CDOM contributes to higher absorption coefficients and compensates for the loss of chl-a due to the degradation of organic matter that added to the CDOM pool at the seafloor [38].

The maps of K_d (Figure 13) showed that location (northern versus southern harbours) and tidal state are the two main drivers of differences in diffuse attenuation in the shallow water regions of the harbour. Considering Tauranga Harbour has several subestuaries, the introduction of extra CDOM from the river mouths might create a gradient toward the estuary mouth (observed in Cussioli et al. [35]). This might explain the higher K_d of pixels at the landward margins compared to those close to the deep channels (Figure 13a,b). Generally, the more urbanised southern harbour had a higher K_d (Figure 13c). About 45% of sediments were introduced to the harbour from the Wairoa River annually, and therefore this region showed a high value of SPM and K_d [39]. The dredged navigational channel in Stella Passage also yielded 25% of sediment loads to the southern harbour and likely had more resuspended sediments from port activities like ship movement [37]. The K_d difference between ebbing and flooding water might result from tidal differences in the resuspension and transport of suspended solids caused by tidal currents and waves (Figure 13c) [40]. K_d associated with ebbing tides tended to be higher than that of flooding tides because the ebbing tidal currents drained off the shallow water regions where wave-induced bed shear stress is higher [35].

The primary production of seagrass and microphytobenthos is not only affected by the total amount of light (light quantity) but also the distribution of energy with wavelength (light quality). Light quality can be affected by the source of sediments [41]. Our results allow us to assess K_d at three different wavelengths (Figure 14), so we can make spatial and temporal estimates of light quality, ultimately allowing understanding of the implications of, for example, terrestrial run-off events.

Limitations

Our study describes a method to determine the bulk characteristics (Figure 15) of an estuary using satellite remote sensing. Our method can be used to complement in situ monitoring for management purposes. Although we did not validate the results due to the difficulty of obtaining direct measurements of K_d, when we made estimates of K_d derived from images from similar adjacent days, the results were consistent (the same was true of the dominant wavelength; details are provided in Figure 12 and Appendix A). The variation of K_d between these images remains small over a short temporal scale and within the same season. In situ measurements of K_d (PAR) by Cussioli et al. [41] also show that K_d does not vary over short timescales. The good correlation between K_d (blue) and K_d (PAR) also supported the potential of using K_d (blue) to estimate K_d (PAR) in ecological studies (Figure 14).

An assumption of the method is that the seabed and water properties do not change substantially between the collections of pixels that are used to optimise Equation (3). To test the effect of these assumptions, we provided additional modelling in Appendix B where random noise was superimposed on known values of K_d and R_w, and then Equation (3) was optimised to determine how well the correct values of K_d and R_w could be retrieved. The magnitude of the error caused by lack of homogeneity was within 10% as long as variations over the set of proximal pixels were less than 10%. The temporal and spatial patterns (or clusters in water and seabed properties [42]) can also increase non-uniformity in the input proximal pixels in our model and increase STD and n [42]. For example, in the temporal methods, if the selected temporal resolution is too long (e.g., annually), long-range variations in the water cycle can contribute to outliers and confuse Lyzenga’s model, and therefore is not an optimal time resolution to use (a shorter one is better).

In the intertidal regions, our estimates of R_b were based on direct satellite observations and were therefore spatially explicit at pixel resolution and capable of detecting variations caused by, e.g., seagrass beds. However, in the subtidal regions, R_b depended on the in situ measurements used to create a model (in our case, the median particle size), which were fairly sparse. Future work will develop better ways to model this and address this weakness. The effect of poor knowledge of R_b was explored in Appendix B, where we showed that underestimates of R_b in shallow water could lead to large errors in R_w and K_d, (a bias of 5% could lead to an error of ~15%), whereas in deeper water (>2 m), the error is much lower. We did not use pixels with seagrass to develop our model. If the method was extended to R_b of seagrass regions, extrapolations of relationships developed for the intertidal to the subtidal may not account for potential species composition differences. In addition, seasonal variations in growth characteristics (e.g., in Ulva and seagrass) [43] and other environmental cycles would not be accounted for in the subtidal R_b. This bias shown in Appendix B would be particularly relevant to seagrass beds due to their lower reflectance in visible bands.

The method will not work well in areas with a very small intertidal range where R_b cannot be measured directly in the shallower waters and used to design a subtidal model. With additional work, it will be possible to develop models based on radiative transfer equations, for example, as provided by software packages such as “Hydrolight”. In the intertidal regions, the errors of the model for estimating R_b for permanently inundated areas originated from the precision of the median particle size extraction and whether the seabed was covered by seagrass. The model was developed using the reflectance of unvegetated exposed pixels because, in our case study, seagrass did not grow subtidally. The existence of subtidal seagrass would create a bias in deeper water; however, the correction becomes quickly unnecessary at these depths. A more sophisticated classification of the sediments may increase the accuracy of the median particle size values used in the regression model. Several machine-learning-based detection techniques (e.g., Ha et al. [22] and Liu et al. [42]) can provide the spatial distribution of vegetation in the exposed regions, which could be used to provide values of R_b for seagrass regions in order to extend our method for use in sites with subtidal seagrass. We also assumed that R_b did not change before and after flooding, an assumption that might also introduce errors because the dry and wet sediments might have different reflectance in each band [44].

This method cannot be applied to images that have been histogram stretched or contrast enhanced because this kind of image manipulation would skew the optical relationship between targets and measurements used to develop Equations (7)–(9). In addition, the spatial method correction requires remote sensing images with a high original spatial resolution. The use of images with a lower resolution, e.g., 20 m or 60 m, increases the mixing of nonuniform reflectance and generates more missing values in the correction process. The images from a medium-resolution spectral imager, such as MODIS and Sentinel-3, are not suitable for the spatial method in these types of estuaries.

5. Conclusions

In this paper, we propose corrections that exploit changes in water depth, either spatially (for good temporal resolution) or temporally (for good spatial resolution), to remove the effect of seabed reflectance on water surface reflectance. For the perpetually inundated regions, the relationship between reflectance in different visible bands and the median particle size were established using exposed information, which was then used (with the inundated sediment texture) to estimate the water bottom reflectance of each underwater pixel needed for correction. The most suitable spatial and temporal resolutions for corrections were 60 m and monthly or seasonally, respectively, within the water depth of 6.4 m. The good correlation between K_d (blue) and K_d (PAR) indicated that satellite-derived information could be used to estimate the production of light-dependent inundated estuarine plants. Future work will be directed towards applying this methodology to other estuaries in New Zealand to evaluate their primary production level and predict the changes if turbidity or water depth fluctuate.