Measuring the effects of bivalve mariculture on water quality in northern New Zealand using 15 years of MODIS-Aqua satellite observations

This study demonstrates that satellite information can be used to quantify the effect of a large coastal mussel farm on water quality. Our study focussed on a 1400 ha area of maricultured green-lipped mussels Perna canaliculus in the Firth of Thames in northern New Zealand. Observations by NASA’s MODIS-Aqua sensor were used to estimate the concentration of chlorophyll a (chl a, as a proxy for phytoplankton), turbidity and sea surface temperature at 500 m resolution. We used approximately 890 single-pass satellite images, each with <50% cloud cover, spanning 2002 to 2017. For each image, the area over the mussel farm was removed and then filled in using kriging with parameters derived from semivariogram analysis of each individual image using linear and exponential models. Three control areas were used to test whether the differences between the kriged data and the observed data in the farm area were significant. We found that, over the 15 yr of the observation, the farm had a significant effect on chl a of −1.6% (reduction) averaged over an area of 1.5 times the area of the farm itself. No significant effect of the farm on turbidity was detected. The farm caused a seasonal pattern of winter warming (+0.03°C) and summer cooling (−0.11°C). The method described here could potentially be applied to other satellite data and other large (multi-kilometre) mariculture facilities to observe their effects on water quality.


INTRODUCTION
Global aquaculture yield is forecast to increase by 2.3% yr -1 between 2017 and 2026 to >100 million t yr −1 (OECD-FAO 2017).Already, aquaculture provides 90% as much human food as wild-caught fisheries, and it is likely to overtake it in the next decade (Bostock et al. 2010).Aquaculture is also big business, with a worldwide value exceeding US$30 billion yr −1 (Kapetsky & Aguilar-Manjarrez 2007).Mariculture of molluscs provides about a quarter of total aquaculture production worldwide, and is dominated by coastal farms for filter-feeding bivalves, especially oysters, scallops, clams and mus -sels (Kapetsky & Aguilar-Manjarrez 2007, Bostock et al. 2010).
Along with the increasing scale and prevalence of mariculture worldwide comes an increasing urgency in managing its effects on the marine environment (Soto et al. 2007, Filgueira et al. 2014).In New Zealand as elsewhere, there is an increasing focus on sustainable, ecosystem-based management of aquaculture (Zeldis et al. 2005, Soto et al. 2007, Bostock et al. 2010, Pinkerton, 2017a).One of the key environmental issues associated with coastal bivalve mariculture is its potential effect on water quality (Zeldis 2005, Soto et al. 2007, Bostock et al. 2010), both in terms of the output of material from the farmed spe-cies to the environment (Soto et al. 2007) and the depletion of phytoplankton by filter-feeding species (Gall et al. 2003, Grant et al. 2007, Zeldis et al. 2013).Managing the environmental effects of mariculture on water quality within an ecosystem context implies long-term monitoring at appropriate time and space scales, which is expensive and complex (Soto et al. 2007).Coastal environments have high spatial and temporal variability which means that in situ sampling alone is unable to adequately assess changes in coastal water quality over large areas and long time periods (IOCCG 2000, Ogilvie et al. 2000, Ogilvie 2000, Grant et al. 2007).
Using hyperspectral airborne remote sensing, Grant et al. (2007) were the first to demonstrate that remote sensing could observe the effect of bivalve mariculture on water quality as a complement to in situ monitoring.Subsequent studies have used data from satellites to quantify changes to water quality by filter feeders (e.g.MERIS: Gernez et al. 2014;Sea-WiFS: Rowe et al. 2015, 2017;MODIS: Warner & Lesht 2015;Sentinel-2: Gernez et al. 2017).Satellite remote sensing has also been used to help identify suitable sites for mariculture (e.g Kapetsky & Aguilar-Manjarrez 2013, Gernez et al. 2014).
Here, we present a new method that uses images from MODIS-Aqua to measure the effect on water quality of a large, coastal mussel farm in northern New Zealand.The farm area studied here is the Wilson Bay Marine Farm Zone (WBMFZ) in the Firth of Thames, a bay of the Hauraki Gulf in northern New Zealand (Fig. 1; Zeldis et al. 2015).This is New Zealand's largest area of contiguous mussel farms (totaling ~1000 ha of farmed space).Environmental resource management of the WBMFZ has in cluded an adaptive management approach called 'Limits of Acceptable Change' (Turner & Felsing 2005, Zeldis et al. 2005), which uses in situ water quality sampling for chlorophyll a (chl a) as a proxy to quantify changes in phytoplankton biomass.The monitoring programme compares monthly chl a samples at a single 'impact' site at the centre of each of Areas A and B (Fig. 1), with control sites surrounding the WBMFZ, located 8−13 km away (Zeldis et al. 2005).While this monitoring programme has been successful in assessing the farm impacts relative to chlorophyll depletion limits (Zeldis et al. 2005, J. Stenton-Dozey pers.comm.), it is limited by its low spatial resolution.In particular, this method is not capable of resolving the true spatial extent of depletion, nor the intensity of depletion at areas other than at its single within-farm impact site.
Satellite remote sensing can help by increasing the spatial and temporal coverage of water quality information (IOCCG 2000, Kapetsky & Aguilar-Manjarrez  The larger area around the farms is called the 'skirt' which is the area used in the present study for kriging (see 'Methods' for details).An oceanographic mooring has been maintained at 36°45.6' S, 175°18.0'E since 1998 2007, Hoepffner et al. 2008)).Satellite observations of ocean colour are now available every day at all latitudes (although they cannot see through clouds) and at reasonable spatial resolution (better than 1 km).Several challenges remain in using remote sensing to assist the development and management of sustainable marine aquaculture (Schaeffer et al. 2013).First, the accuracy of satellite observations of water quality is uncertain, especially in near-shore environments which tend to be optically complex (i.e. the colour of the water is affected by many, independently varying components such as phytoplankton, sediment and coloured dissolved organic matter, CDOM; IOCCG 2000).Second, optical properties can change over time (both on short and long time scales), which potentially impacts the utility of satellite data for monitoring long-term change (Dierssen 2010).Third, the spatial resolution of satellite data can limit their utility for observing the effects of smaller marine farms.
In this study, we used satellite data from MODIS-Aqua at 500 m resolution to directly observe the effect of the large mussel farm in WBMFZ on 3 important water quality variables: (1) the concentration of chl a; (2) turbidity, which responds to particulate material in the water column; and (3) sea surface temperature as an indicator of water column mixing.

Study area
The Firth of Thames is an estuary which covers 1100 km 2 with a mean depth of 14 m (Fig. 1), and is subject to inflows of water from the Waihou and Piako Rivers to the south (combined catchment area of 4200 km 2 ), and mixing to the north with the larger Hauraki Gulf which opens onto the New Zealand continental shelf and the southwest Pacific Ocean (Zeldis et al. 2015).Chl a in the Firth of Thames varies between about 0.5 and 5 mg m −3 , and average daily primary productivity rates are ~1 g C m −2 d −1 (Gall & Zeldis 2011).Mixed layer depths are about 20 m for most of the year (Gall & Zeldis 2011, Zeldis et al. 2015).Annual riverine input of sediment is estimated to be 195 000 t yr −1 (Hicks et al. 2011).Tides in the Firth of Thames are semi-diurnal, with average spring and neap tidal ranges of 3.2−3.5 m and 2.0− 2.2 m, respectively (Healy 2002).Over a tidal cycle, currents within the Firth of Thames are typically 0.3 m s −1 , although speeds of up to 1.4 m s −1 have been recorded (Proctor & Greig 1989, Black et al. 2000); residual tidal currents are weak (0.01− 0.02 m s −1 ; Broekhuizen et al. 2002).In terms of bio-optics, Area A is intermittently strongly affected by sediment and CDOM from riverine run-off leading to 'case 2' conditions where water colour is not principally affected by phytoplankton (Morel & Prieur 1977).Satellite data over the farm (this study) suggest total suspended matter concentrations of 2.8 g m −3 (2.1−3.5 g m −3 ; median, lower and upper quartiles) and chl a concentrations of 3.6 g m −3 (2.8− 4.6 g m −3 ).
In the WBMFZ aquaculture facility, mussel spat are seeded onto vertical lines, and grow by feeding on naturally occurring marine phytoplankton.No additional feed or enrichment is used.Harvesting occurs approximately 18 mo to 2 yr after seeding, and stocking densities have no consistent seasonal pattern.Development of Area A (5.8 × 2.4 km) started in 2001 and its current (2016) harvest was 16 400 t yr −1 of mussels.Area B (5.8 × 1.9 km) development started in 2014 and it currently produces approximately 4500 t yr −1 (J.Stenton-Dozey pers.comm.).Because Area A is presently much more fully developed than Area B, this paper is focussed on the former, which will dominate the farm effects on water quality.Based on the density of mussel lines, stocking densities and nominal clearance rates, Broekhuizen et al. (2004) estimated that phytoplankton within Area A could suffer a mussel-induced mortality rate of 30% d −1 but went on to say that in practice it is likely to be much less.
The ecosystem effects of the mussel farm in the Firth of Thames is of concern because the region supports important ecosystem services (MacDiarmid et al. 2013) including spawning and nursery grounds for New Zealand's largest inshore fin-fishery (snapper Pagrus auratus; Zeldis & Francis 1998), and is an important site for shorebirds (wetland of international importance under the Ramsar Convention).

Overview of methods
The methodology was based on the following steps.(1) We brought together satellite data products processed to satellite data Level 2 (individual satellite overpasses) exclusively using case 2 algorithms with strict quality control to eliminate clouds and poor quality data (see 'Satellite data processing' below for more details).(2) We investigated the quality of the satellite dataset using in situ measurements (see 'Validation of satellite data' below).
(3) Statistical analysis was applied to each satellite image with > 50% valid data in the WBMFZ area (Fig. 2).We were concerned with estimating the effect of the farm on water quality in the 'core' area.This core area is centred on Area A but could be smaller or larger than the actual farm area, as the effect of the farm on water quality may extend outside the physical boundaries of the farm because of water movement and mixing (e.g.Gall et al. 2003, Broekhuizen et al. 2004, 2005).We call the larger area outside the core the 'skirt'.The area of the skirt was much larger (125 km 2 , 400 pixels) than Area A (14 km 2 , 56 pixels).The core and skirt areas are almost always optically deep (satellite data suggested that 99.8% of pixels over the 15 yr period had optical depths >1).We tested the effect of changing the size of the skirt and using a constant area of skirt irrespective of the size of the core area and found that neither significantly changed the results; hence, we used a fixed skirt size.Semivariogram analysis 532 Fig. 2. Overview of the method used in this paper.(1) Select each individual image with good data coverage in the area around the farm (example file here is the chl a variable from file A20022830215 at 500 m spatial resolution).(2) Only use data in the farm and surrounding 'skirt' area.(3) Carry out semivariogram analysis to obtain parameters for kriging.Colour: density of points; solid line: fitted model; dashed lines: fitted range and sill.(4, 5) Using only the skirt data, apply kriging over the farm area.(6, 7) Compare the kriged data and actual satellite measurements in the farm area to obtain spatial estimates of the effect of the farm was carried out on data in the core and skirt together to obtain semivariogram model parameters for linear and exponential models (see 'Spatial analysis' be low).We masked out the data in the core and used semivariogram parameters and isotropic kriging to fill in the satellite data in the core.
We propose that these kriged data represent the best estimate of water quality in the core if there was no farm present.Under this assumption, comparing the kriged data with the observed data gives a measure of the effect of the farm over the area of the core.
We tested the effect of anisotropy and found this had little effect (see Discussion).Three areas of different sizes were used as controls to investigate the uncertainty in these estimates of the farm effect (Fig. 3).These control sites are the same size and shape as the core and were analysed in exactly the same way.All analyses in this study were carried out using IDL software version 8.5 (Interactive Data Language, Harris Geospatial Solutions).

Satellite data processing
Satellite data used in this study were from NASA's Moderate Resolution Imaging Spectrometer of the Aqua satellite (https://ocean color.gsfc.nasa.gov/data/aqua/).MO DIS data have a spatial resolution of 250− 1000 m in visible spectral bands, frequent overpasses (daily) and long-term tracking to characterise long-term change in sensor sensitivities.Data from MODIS-Aqua used in this study cover the period from July 2002 until February 2017.Level 1A (top of atmosphere, uncalibrated) MODIS-Aqua data were acquired either by file transfer from NASA (data between 2002 and 2007) as full spatialresolution 5 min granules, or as direct broadcast data by the NIWA X-band receiver (after 2007).All direct broadcast data were calibrated and processed using NASA Collection 6 calibration files.In total, we processed the 6949 data files using NASA's SeaDAS v7.2.Data were rejected for land, cloud cover, solar glint, white-cap reflection, atmospheric correction failure or in-water algorithm failure.Satellite data products were calculated at a spatial resolution of 500 × 500 m and derived variables were projected to a transverse Mercator grid.

Sea surface temperature
The MODIS-Aqua daytime sea surface temperature (SST) product was used in this study with no local adjustment, but with 1 km resolution data resampled onto a 500 m resolution grid.Satellite estimates of surface water temperature were excluded based on the MODIS data quality flag 'sstqual' and when the chl a processing failed (potentially indicative of residual cloud contamination).1], which is only lightly stocked and its biogeochemical effect is not considered important).For each site (A, C, D, E), the white area shows the outer extent of the 'skirt' -the area surrounding the farm that was used for kriging inside the site.The 9 black rectangles inside the white areas show the different sizes of inner areas, i.e. the cores, which have areas of 0.75, 1, 1.25, 1.5, 2. 2.5, 3, 4 and 5 times the area of Farm A. The shapes of the skirts and cores were determined based on the tidal ellipses in the region and models of dispersion over a tidal cycle (Gall et al. 2003) Chl a and turbidity Given that the study area is intermittently case 2, the satellite processing must be able to separate the co-occurring phytoplankton, suspended sediment and CDOM.Failure to distinguish between these coloured materials would result in errors in their spatial distributions and could lead to erroneous conclusions as to the effect of the farm on water quality.A variety of atmospheric correction methods appropriate for turbid (case 2) waters are available in SeaDAS for use with MODIS data (Ruddick et al. 2000, Wang & Shi 2007, Bailey et al. 2010) and others are under development (e.g.Lavender et al. 2005).It is not clear which of these performs best in a given situation, and here, we used the near infrared−short-wave infrared (NIR-SWIR) combined method (Wang & Shi 2007) which may be the most robust when turbidity is high and local information on aerosol optical properties is not available.In the approach of Wang & Shi (2007), NIR wavelengths are used to estimate reflectance by atmospheric aerosols over clear waters, and SWIR wavelengths (that are not affected by allochthonous particulates) are used over turbid waters.
Simple empirical methods (e.g.O'Reilly et al. 1998) can robustly estimate chl a concentration in case 1 waters but typically perform poorly in case 2 waters (IOCCG 2000, Pinkerton et al. 2006, Zeldis et al. 2015).This study used the quasi-analytical algorithm (QAA) (Lee et al. 2002(Lee et al. , 2009) ) to estimate particulate backscatter at 555 nm [b bp (555)] and phytoplankton absorption at 488 [a ph (488)].MODIS data were interpolated to 500 m using paired bands where necessary (Franz et al. 2006).Case 2 chl a was obtained from a ph (488) using chlorophyll-specific absorption coefficients measured in New Zealand estuaries (Pinkerton 2017b).We blended the case 2 chl a and the MODISdefault chl a (which is likely to be more robust in case 1 conditions) using a logistic-scaling of b bp (555) (Pinkerton 2017b).Values of b bp (555) were also linearly scaled to calculate turbidity (in normalised turbidity units) using in situ bio-optical measurements in New Zealand estuaries (Pinkerton 2017b).

Validation of satellite data
Satellite data products were compared with in situ measurements to explore data quality.A biogeochemical mooring has been maintained in the Firth of Thames at 40 m depth since 1998, with approximately quarterly servicing (Gall & Zeldis 2011; Fig. 1).The moorings included internally logging thermistors (Onset) and an integrating natural fluorometer (INF-300, Biospherical) at approximately 7 m depth.Following correction for biofouling and averaging over a day to account for diel effects, INF data were used in conjunction with shipboard biogeochemical measurements to estimate daily chl a (Gall & Zeldis 2011, Zeldis et al. 2015).
Time-series plots show that the satellite observations capture long-term variations in the mooring measurements of chl a and surface temperature (Fig. 4).'Match-up' analysis -where the median of 9 satellite observations at the mooring site are compared to daily average measurements from the mooring for the same day -is widely used to assess satellite data quality (Darecki & Stramski 2004, Chang & Gould 2006).Match-up analysis yielded good agreement for temperature (R 2 = 0.95) but poor results for chl a (R 2 = 0.23 in log space; Fig. 5).There were no measurements of turbidity from the mooring, so instead we compared the satellite estimate of turbidity with in situ measurements at a nearby site (Manukau Harbour, northern New Zealand).Unpublished data show that suspended sediment had very similar inherent optical properties in the 2 areas (Pinkerton 2017b) meaning that the accuracies of satellite estimates of turbidity in the 2 areas are likely to be similar.Satellite observations of turbidity in Manukau Harbour were found to be unbiased with respect to in situ measurements, but explained only a small proportion of the variance (R 2 = 0.088 in log space; Fig. 5).This is likely to be the case for satellite data on turbidity in the study area also.
Low R 2 values between satellite and in situ observations of chl a and turbidity are not fundamentally limiting for the method used in this study.However, inaccuracies in the remotely sensed spatial distributions of coloured material will affect the method developed here, and this is considered further in the 'Discussion'.

Spatial analysis
Many statistical methods are potentially capable of analysing the spatial patterns of a 2-dimensional data field, including kriging, spline-based approaches and 'hot-spot' analyses (Krige 1951, Wahba 1990, Getis & Ord 1992).We used kriging to 'fill in' data over the core (or control) areas as the method is well established, simple to apply and computationally efficient, and it produces unbiased estimates (Papritz & Stein 1999).Semivariogram analysis was used to determine the appropriate parameters (nugget, sill, range) in the kriging model.The semivariogram γ (h) is half the average squared difference in a variable (z) measured at 2 points (x i and x i +h) separated by a distance h (Matheron 1963); Eq. ( 1).The semivariogram measures the way in which pairs of observations in a spatial data field become more different as their separation increases, as has been used to help interpret changes in geophysical data (Isaaks & Srivastava 1989) including satellite images of water quality (Bierman et al. 2011).
(1) Kriging can be used with a number of different covariance models, including linear, exponential, Gaussian and spherical.The covariance data found here (see example in Fig. 6) best fitted linear and exponential models, with some images better fitting the former and others the latter.Consequently, we carried out all analysis using both of these models, and combined the results using the method given in the following subsection.Satellite measurements of chl a and turbidity were log transformed before analysis to make their distributions more normal, but temperature data were not transformed.

Statistical analysis and confidence
For each variable (chl a, temperature, turbidity), there are 8 sets of results.These results correspond to 4 core areas (Area A 'farm', and 3 controls C, D and E), each based on 2 different semivariogram models (linear and exponential).We first define the 'effect' (Eq.2).For chl a and turbidity (which tend to be lognormally distributed), the effect is defined as the median over the core area of the ratios between observed data (z observed ) at location x i and the kriged data at the same location (z kriged ) minus 1.The chl a and turbidity effects are dimensionless.For temperature data (which are closer to a normal distribution), the effect is defined as the median of the differences in the core area between observed data and kriged data.The temperature effect has units of °C.Using the median rather than the mean removes the sensitivity of the measure to extreme values.Negative effects indicate a 'hole' (i.e.reduction in chl a due to the farm, or to cooling), whereas positive effects indicate enrichment due to the farm or warming.
Next, we define the best estimate of the effect of the farm (written here as A μ ) as the mean of the overall effects estimated using the linear and exponential kriging models (A lin , A exp , respectively) minus the mean of the effects found at controls C, D and E (Eq. 3).Upper and lower confidence bounds are estimated as A μ ± 2A σ , where A σ is estimated as the standard deviation of the effect at the control sites by the 2 kriging models (Eq.4): (4)

Trend analysis
Temporal trends in water quality variables were determined using a Mann-Kendall test (Hipel & Mc Leod 1994).Estimates of farm effect from all individual satellite overpasses were used for the trend analysis with no time compositing.Data were deseasonalised by collapsing measurements into a single annual period, smoothing using a moving 2 mo window and interpolating based on day of year.Statistical significance was assessed using the unadjusted Mann-Kendall statistic Z MK (Hipel & McLeod 1994).Trends were identified as significant when the p-value was < 0.01 and the 95% confidence interval of the Sen slope did not intercept 0 (Larned et al. 2015).

RESULTS
The final analysis was based on data from approximately 890 MODIS-Aqua files, with, overall, more than 23 000 pixels in the area of the WBMFZ for each product (Fig. 7).Across all 3 variables, the effects observed at the control sites were always clustered about 0 and were little affected by the size of the core (chl a: 0.003 ± 0.007 (dimensionless); temperature: −0.002 ± 0.008°C; turbidity: −0.0001 ± 0.008 (dimensionless); mean ± SD across core sizes 0.75 to 5 times Area A).In other words, the analysis method developed here behaved in a stable way over different sizes of core at the control sites and tended to 0 effect.Almost without exception, the effect due to the farm was greater in magnitude than the effect at any of the 3 control sites (Fig. 7).The effects on all 3 variables due to the farm, predicted by the linear model, were greater than the effects determined by the exponential model.The direction of the effect of the The effect of aquaculture at the farm on chl a was significant when the area of effect (core size) was 1.5 times the actual area of the farm (t (3) = 4.12, p = 0.026).The effect of the farm on chl a at this spatial scale was −1.6%.At scales smaller than the farm area or larger than twice the farm area, the effect of the farm on chl a was not significant (t (3) < 2.88, p > 0.064).In contrast, the annual-average cooling effect (−0.03°C) of the farm was significant at all spatial scales of analysis (t (3) > 4.42, p < 0.021).No significant effect of the farm on turbidity was found (t (3) < 3.62, p > 0.036).
Because the satellite data extended across years between 2002 and 2017 and had good coverage in all months, we used this analysis to explore seasonal and interannual effects of the farm on water quality (Fig. 8).We carried out this analysis using a core size of 1.5× the area of the farm based on the result for chl a (Fig. 7).For chl a, we found that the effect of the farm was significant in January, May, June, August and September (t (3) > 4.0, p < 0.027), and not significant otherwise (t (3) < 0.90, p > 0.054).In contrast, the effect of the farm on surface water temperature had a strong seasonal variation, with maximum warming in June (+ 0.03°C, t (3) = 8.8, p = 0.003) and maximum cooling in January (−0.11°C, t (3) = 18.5, p < 0.001).The effect of the farm on turbidity was positive (+ 2.0%) in May and June (t (3) > 4.4, p < 0.021), but not significant at other times of the year or in the annual average (t (3) < 2.9, p > 0.064).
In terms of long-term trends between 2002 and 2017 (Fig. 9), we found that the effect of the farm on chl a had a significant trend (Z MK = −2.3,p = 0.02) with a Sen slope of 0.12% yr −1 (increasing effect of the farm).Trends were not significant in temperature (Z MK = −0.67,p = 0.50) or turbidity (Z MK = 1.2, p = 0.25).
The computing resources required for this analysis were reasonable but not prohibitive.The analysis of each image (which involved data mapping, semivariogram analysis and kriging) was completed in about 0.34 s, and each set of 890 files took ~5 min.With 288 runs (4 variables, 2 semivariogram models, 4 areas [Area A and 3 controls] and 9 spatial scales), the analysis took about 24 h of processing resource.Sensitivity analysis (including for anisotropy) took ~5 d of processing.

DISCUSSION
We used satellite data to observe the effect of a large mussel farm on key water quality parameters, namely the concentration of chl a (as a proxy for phytoplankton abundance), turbidity and SST.Satellite observations in the vicinity of the farm were removed, filled in by kriging and then compared with the observed data in this core area.Overall, our results suggested that the effects of aquaculture in WBMFZ Area A on water quality were small.We found that chl a was reduced by 1.6% (long-term median value) over an area of about 1.5× the actual area of the farm.We detected no significant effect of the farm on turbidity.
The size of the effect found here is smaller than suggested by both field sampling and modelling.A previous synoptic study of the effect of the farm on phytoplankton using towed equipment suggested chl a depletion at the farm scale of between 5 and 10% (Gall et al. 2003).Previous numerical modelling studies (Broekhuizen et al. 2004, 2005, Stenton-Dozey et al. 2005) (b,d,f) Mean effect of the farm and 95% confidence intervals (grey area; see 'Methods' and 'Results' for more details).Dashed lines: 0 effect 8−10% within an area extending a few km beyond the boundaries of Area A. The modelling also suggested that the farm could lead to mild enhancement of phytoplankton growth in the summer, something that was not detected in the present study or by field sampling (Gall et al. 2003).
We also found an effect of the mussel farm on surface water temperature, with a consistent pattern of winter warming and summer cooling of up to about 0.1°C.It is likely that the tidal currents moving past the physical structure of the mussel farm break up the surface thermal microlayer, and this caused a change in surface water temperature around the farm.This is consistent with Gall et al. (2003), who found evidence of significant reduction in tidal current speed (20−25%) in the lee of the WBMFZ Area A. In summer, the surface microlayer will be warmer than the underlying water, so the farm's stirring effect leads to a cool anomaly.In winter, when the air temperature is cooler than the water, a cool surface microlayer is formed, and the disturbance of this is likely to lead to the observed warm anomalies around the farm.
The small change in surface temperature detected here suggests that the physical disturbance caused by the farm structure is not sufficient to mix up water from below the pycnocline (typically located at ~20 m depth) to the surface.If this was the case, we would  ,d,f) by year for core area of 1.5 × the area of the farm, and 3 variables: a,b: chlorophyll a concentration; c,d: sea surface temperature; e,f: turbidity.The grey areas show the 95% confidence intervals estimated as A μ ± 2A σ , where A σ is estimated as the standard deviation of the effect at the control sites by the 2 kriging models (Eq.4).
The dashed lines show 0 effect expect to see much larger effects of the farm on surface temperature, as the pycnocline in the Firth of Thames is associated with a 2°C vertical temperature difference over several months (Zel dis et al. 2015).This in turn argues for the reduction in chl a around the farm being due to phytoplankton depletion as a result of filter feeding of mussels rather than being a dilution effect (e.g.changes to the vertical distribution of phytoplankton).
There was considerable variability in the size of the farm effect identified by this method based on different satellite images (Fig. 9).The method estimated both positive and negative effects of the farm for all variables, and in about 20% of the images, the effects were estimated as being greater than 0.1 in magnitude.Some of this variability is undoubtedly real and some will be a methodological artefact.Mussel farms can genuinely have positive and negative effects on concentrations of phytoplankton and suspended particulate matter (e.g.Ogilvie 2000), but, in line with other applications of remote sensing, a large amount of satellite data is needed to offset the relatively high uncertainty.
In this context of variability, it is important to note the use of 3 control sites to test for significance.The way in which the effect varied with the size of the core agreed with our expectations.When the core area for analysis was smaller than the actual area of effect, values in the skirt would be affected by the aquaculture activity and this would tend to make the estimated effect too small.If the core area chosen for analysis was much bigger than the actual area of effect of the farm, the median effect estimated would again tend to be too low.In the latter case, we would also expect the potential for spurious results to in crease so that the uncertainty bounds (based on the control sites) would likely increase.We found that our results corresponded closely to this expected pattern for chl a (Fig. 7).The size of the effect of the farm increased, then decreased as the spatial size of the core increased, and the confidence bounds became wider.The maximum effect on chl a was observed at an area of 1.5 times the area of the farm (WBMFZ Area A), which suggested that this was the appropriate spatial scale for analysis.To further test that the effects are real, we repeated the analysis switching round the results of Farm A with each of the control areas in turn.In each case, the size of the effect did not vary in an interpretable and consistent way as the size of the core changed.In addition, the effect identified from the satellite data was consistent with numerical modelling of the potential chl a depletion 'halo' of the farm (Gall et al. 2003;Fig. 10).The shapes of the effects were very similar, although the magnitude of the depletion in chl a predicted by the analysis of Gall et al. (2003) was greater than the reduction in chl a around the farm observed here.
The same pattern in how the effect of the farm changes with core size was broadly seen for the temperature and turbidity variables but was not so clear.For example, the confidence bounds for temperature increased as the size of the core increased, but the greatest effect for temperature was found at 4× the size of Area A, not at 1.5×.For turbidity, the confidence bounds narrowed as the core size increased beyond 3×.Phytoplankton have a capacity for re -growth following depletion by filter feeding (i.e.chl a will demonstrate non-conservative behaviour), whereas temperature and turbidity effects may be transported conservatively by water movement.Mean hydraulic residence times in the farm area are likely to be on the order of 1 d (Gall et al. 2003), similar to phytoplankton intergenerational times.This may account for the smaller spatial scale of effect seen for chl a relative to SST and turbidity (Fig. 7).
We discount the possibility that changes observed in the satellite data in the vicinity of the farm were due to the above-surface farm infrastructure (such as floats and platforms) affecting the spectrum of waterleaving radiance and hence the derived satellite products.(1) The proportion of the area covered by floats in the farm area is small (< 0.5% of the area).
(2) The reflectivity of the floats is similar to that of water, so their individual effect on water-leaving radiances will be modest.(3) The changes in water quality detected in our analysis extend outside the outer limits of the farm superstructure.There was some suggestion in our results of an increasing effect of the farm on all variables, with a significant trend identified in the effect of the farm on chl a (Sen slope 0.12% yr −1 ).The stocking density of mussels in the farm has increased over the period analysed: by mid-2003 it was ~66% seeded and in 2016 it was ~92% developed.It is hence plausible that the small increasing trend in chl a depletion follows from the increase in stocking density over time.The trends in the effect of the farm on surface temperature and turbidity were not significant, consistent with little change in the physical structure of the farm over this period and little effect of the mussels on suspended particulate matter concentrations in the upper water column.Changes in the regional hydrography of the Firth of Thames over the study period could lead to spurious trends being identified by this method.However, trends in all variables at the 3 control sites were not significant, so there is no evidence for long-term background environmental change over the course of the study.
Given the difference in results obtained using a linear compared to an exponential kriging model, it is reasonable to ask whether our results are robust to uncertainties in the assumed shape of the semivariogram function.Certainly fitting the parameters (nugget, sill, range) to semivariograms of remotely sensed water quality variables was not straightforward.As can be seen in Fig. 6, the individual semivariance data points were highly scattered and the smoothed lines did not unambiguously follow any of the standard semivariogram forms (linear, exponential, Gaussian, spherical) over the whole distance range.Anisotropic semivariogram analysis may improve the fit somewhat, but agreement with models will always be imperfect.The relatively coarse spatial scale of the satellite data also misses substantial small-scale variability in water quality properties.We developed an iterative fitting method for both the linear and exponential forms, where the semivariance beyond a certain distance was not used.Both the fitted sill and range parameters of the linear models agreed reasonably well with those in the exponential model.The median fitted sills were 0.177 (linear) and 0.182 (exponential), and the median fitted ranges (in pixels) were 17.3 (linear) and 27.2 (exponential).We set the nugget in the exponential model to 0 to ensure reasonable values of range, whereas the median nugget in the fitted linear models was non-negligible (0.005).Carrying out the analysis using a linear model with the nugget set to 0 had little effect on the results.We also tested the effect of using fixed values of range (set to 15, 20 and 25 pix-els) in both the linear and exponential models.The re sults using the fixed ranges bracketed those obtained using the fitted semi variogram parameters as expected.It seems there fore that the difference between the results with exponential and linear models was a result of the way in which the models transmitted spatial variations from the skirt into the core area.
We were not comfortable with using different models for different images, i.e. using a linear model for one image and an exponential model for a different image, as it was usually not clear from the fitting which model best described the semivariance data.Hence, we consider it appropriate to use both linear and exponential models with equal weighting.For future development of this method, it would be useful to develop a better way of combining results from different kriging models.For example, a resampling method such as the jackknife or bootstrap using withheld data could be used to evaluate the performance of different semivariograms and kriging models for each individual image.We also note that kriging can provide estimates of the uncertainties of the predictions.Combining these estimates of uncertainty would allow different estimates of the effects of the farm to be weighted in Eqs. ( 3) and ( 4).
The spatial structure in the water quality parameters are not likely to be isotropic as we have assumed, and there will probably be more variation in water quality in the inshore−offshore direction than in the along-shore direction.We tested the effect of this anisotropy using kriging with anisotropic factors of α = 1.2 and 2.0.A factor of α = 1.2 indicates that the length-scales of variation in water quality parameters are 20% shorter in the inshore−offshore direction than in the along-shore direction, and a factor of α = 2.0 means that the cross-shore length-scale is half the along-shore scale.Qualitatively, we found that the anisotropy had negligible effect on the results; the observed effects of the farm as a function of core area (Fig. 7) with and without anisotropy had the same shapes (the correlation coefficients between the medians and confidence intervals were always > 0.95).However, anisotropy did affect the magnitudes of the estimated effects.We found that anisotropy with α = 1.2 and 2.0 reduced the estimated effect of the farm on chl a by 5 and 24%, respectively (with the change expressed as a proportion of the isotropic effect), reduced the estimated effect of the farm on temperature by 4 and 14%, and increased the estimated effect of the farm on turbidity by + 97 and + 230%.Despite the big increases of anisotropy on the estimated effect of the farm on turbidity, the effect of the farm was still not deemed significant because the confidence intervals still spanned 0 effect.These results are consistent with expectations that anisotropy would be greater for turbidity than for either chl a or temperature because of potentially strong cross-shore gradients in suspended sediment concentration.It is not easy to estimate the appropriate level of anisotropy, which will also vary from day to day.Hence, we conclude that our results are qualitatively and approximately quantitatively robust to the effect of anisotropy, but recommend that fully anisotropic kriging (i.e. with anisotropic factors derived separately from each image) be developed in the future to improve confidence in the sizes of the farm effects.
Although the temperature data from MODIS-Aqua were found to be accurate, much poorer accuracy was achieved for chl a or turbidity (Fig. 5).Similar (poor) results have been found in this kind of 'matchup' analysis for chl a and turbidity in coastal areas elsewhere (e.g.Pinkerton et al. 2003, Darecki & Stramski 2004, Chang & Gould 2006).The mismatch between in situ and satellite measurements of chl a and tur bidity potentially have many causes (Pinkerton et al. 2003), including: (1) water column structure: although the water column at the mooring site in the Firth of Thames is usually well-mixed to below the depth of the mooring sensor, there are occasions when chl a at depth differs substantially from that at the surface (Gall & Zeldis 2011); (2) fundamental differences in the spatial and temporal scales of observation by satellites and moorings (Pinkerton et al. 2003, Gall & Zeldis 2011).We integrated mooring measurements over a 24 h period to better reconcile the scales of observation but, even so, the 'match-up' comparison method is limited as a way of assessing satellite data quality; (3) in situ measurements of chl a by the INF are inexact (Gall & Zeldis 2011); (4) satellite algorithms for chl a and turbidity have high uncertainties in the coastal zone (Garver & Siegel 1997, IOCCG 2000, Lee et al. 2002, Lohrenz et al. 2003).These uncertainties are due to variability in the inherent optical properties of coastal water constituents coupled with inaccuracies in the remote sensing of water-leaving radiance, for example due to imperfect atmospheric correction (Pinkerton et al. 2003).It is hence important to de velop methods of analysis of satellite data that are robust to limitations in the absolute accuracy of measuring water quality from satellites and in our methods to quantify these accuracies.
The method developed here is likely to be relatively insensitive to uncertainties and errors in the absolute values of chl a and turbidity derived from ocean colour satellite information, because it is based on analysis of the spatial patterns in the data.In particular, our method will be insensitive to long-term variations in the absolute accuracy of satellite measurements of water quality over decadal scales, which could potentially compromise traditional methods of monitoring (e.g.Dierssen 2010).However, inaccuracies in the satellite-derived spatial patterns of chl a and turbidity will affect our results and conclusions, as the method developed in this study relies on the spatial structure in the satellite observations.For example, if the effects of sediment, phytoplankton and CDOM on water colour were not properly resolved, inaccuracies in both spatial patterns and in the absolute values of chl a and turbidity could follow.We found that the satellite estimates of chl a and turbidity in the skirt area were effectively independent (R 2 < 0.01).This result suggests that the satellite processing was working effectively, as systematic failure in the algorithm would be expected to lead to correlation between chl a and turbidity.
A useful future study would be to examine applying the method described in the present study to higher spatial-resolution satellite data (e.g.MERIS, LANDSAT, Sentinel-2; e.g.Gernez et al. 2014Gernez et al. , 2017) ) to investigate smaller aquaculture facilities.The developing field of hot spot analysis (Getis & Ord 1992, Harris et al. 2017) also promises to provide better statistical tools which could potentially be used to investigate spatial effects of aquaculture without the need to define skirt and core areas.
In conclusion, while both modelling and in situ monitoring have been used to quantify the effects of the WBMFZ Area A mussel farm on water quality (e.g.Gall et al. 2003, Broekhuizen et al. 2004, 2005, Zeldis 2005), this is the first study to use satellite data to directly observe the effects of such a farm.The approach given here delivers direct information on 3 of the key water quality parameters that are required to monitor and manage the environmental effects of mariculture: phytoplankton, temperature and suspended particulate matter (Kapetsky & Aguilar-Manjarrez 2007, Soto et al. 2007, Ferreira et al. 2010).The information is likely to be relevant to monitoring several more key water quality properties, including nutrients, farm-derived detritus, dissolved oxygen and microbial loading.Remote sensing should be thought of as a complement to in situ observation of the environment (e.g.Gernez et al. 2014Gernez et al. , 2017)), direct monitoring of harvest species (e.g.Filgueira et al. 2014) and bio-physical modelling (Broekhuizen et al. 2005, Guyondet et al. 2013, Rowe et al. 2017).

Fig. 1 .
Fig. 1.Wilson Bay Marine Farm Zone (WBMFZ) in the New Zealand Firth of Thames, at the southern extreme of the Hauraki Gulf.The locations of the aquaculture farms in WBMFZ Areas A and B are shown as rectangles.The larger area around the farms is called the 'skirt' which is the area used in the present study for kriging (see 'Methods' for details).An oceanographic mooring has been maintained at 36°45.6' S, 175°18.0'E since 1998

Fig. 3 .
Fig. 3. (a) Farm A; (b−d) 3 controls (labelled C, D and E to distinguish them from Farm B [see Fig.1], which is only lightly stocked and its biogeochemical effect is not considered important).For each site (A, C, D, E), the white area shows the outer extent of the 'skirt' -the area surrounding the farm that was used for kriging inside the site.The 9 black rectangles inside the white areas show the different sizes of inner areas, i.e. the cores, which have areas of 0.75, 1, 1.25, 1.5, 2. 2.5, 3, 4 and 5 times the area of Farm A. The shapes of the skirts and cores were determined based on the tidal ellipses in the region and models of dispersion over a tidal cycle(Gall et al. 2003)

Fig. 4 .
Fig. 4. Biogeochemical variables from MODIS satellite (grey diamonds and lines) between 2002 and 2017.(a) Chlorophyll a concentration; (b) surface temperature.In situ measurements of chl a concentration and surface temperature from the mooring are shown as black dots for comparison

Fig. 5 .
Fig. 5. Validating key satellite products by match-up analysis: (a) chlorophyll a concentration; (b) sea surface temperature; (c) turbidity.Solid black line: y-on-x regression line; elipses: 80 and 95% confidence limits; dashed line: 1:1 line.Also shown: no. of points (N), and coefficients of determination (R 2 ) in linear and log-log space.Match-up analysis yielded good agreement for temperature, but poor results for chl a and turbidity; see 'Methods' for more information

Fig. 6 .
Fig. 6.Example of semivariogram analysis using the chlorophyll a concentration data from image A20021940220.(a) Semivariogram showing the change in semicovariance between pairs of points as the spatial distance between them changes.Red (blue) colours indicate high (low) density of points.Black dots are individual semivariogram points.The black line was fitted by median smoothing over a distance of h = 3.(b) Fitting the kriging parameters.The grey line is the smoothed semivariogram (SV) data, to which linear (lin) and exponential (exp) models are fitted (black lines).The fitted nugget (N), sill (S) and range (R) are shown.Note that the nugget is forced to 0 for the exponential model Fig. 7. Results of changing the size of the core area for each variable: (a,b) chlorophyll a concentration; (c,d) sea surface temperature (SST); (e,f) turbidity.In each case, the median effect across all months and all years is shown on the y-axis.'Area' on the x-axis is the size of the core area in the analysis expressed as a proportion of the size of Farm A. (a,c,e) Individual results for 2 different models (lin: linear; exp: exponential), for the effect of WBMFZ Area A (thick lines) and for control sites C, D and E (thin lines).(b,d,f) Mean effect of the farm and 95% confidence intervals (grey area; see 'Methods' and 'Results' for more details).Dashed lines: 0 effect

Fig. 8 .
Fig.8.Effects of farm (a,c,e) by month and (b,d,f) by year for core area of 1.5 × the area of the farm, and 3 variables: a,b: chlorophyll a concentration; c,d: sea surface temperature; e,f: turbidity.The grey areas show the 95% confidence intervals estimated as A μ ± 2A σ , where A σ is estimated as the standard deviation of the effect at the control sites by the 2 kriging models (Eq.4).The dashed lines show 0 effect

540Fig. 9 .
Fig. 9. Trend analysis of effects of the farm on (a) chlorophyll a concentration; (b) sea surface temperature; (c) turbidity showing all individual satellite estimates of effect, deseasonalised as described in 'Methods'.The Mann-Kendall trend (red line), number of points (N), Mann-Kendall statistic (Z) and significance (p) are shown.In no case was the trend significant at the 1% level 541

Fig. 10 .
Fig. 10.Spatial distribution of chlorophyll a depletion.Values shown are the ratio of chl a with farm present to chl a with no farm; thus, values <1 indicate a reduction in phytoplankton associated with the farm.(a) Depletion 'halo' over a 24 h period based on synoptic surveys around the farm (reproduced from Gall et al. 2003, their Fig.29).Black rectangle: aquaculture farm; white background: 0 predicted depletion.(b) Spatial distribution of chl a reduction derived in the present study using a core size of 2.5× the area of the farm.A bigger area than 1.5× is shown here to illustrate the spatial shape of changes in chl a around the mussel farm.Light grey background indicates that no measurement exists (as this is outside the kriged area); black rectangles: outer limit of the skirt (largest), 1.5× the area of the farm (intermediate) and the area of the farm (smallest rectangle).Note that the colour scales in (a) and (b) differ to emphasize spatial patterns rather than absolute values