Toward adaptive robotic sampling of phytoplankton in the coastal ocean

Gaussian process models embedded on an AUV enabled autonomous tracking and mapping of phytoplankton biomass in three dimensions. Currents, wind, bathymetry, and freshwater runoff are some of the factors that make coastal waters heterogeneous, patchy, and scientifically interesting—where it is challenging to resolve the spatiotemporal variation within the water column. We present methods and results from field experiments using an autonomous underwater vehicle (AUV) with embedded algorithms that focus sampling on features in three dimensions. This was achieved by combining Gaussian process (GP) modeling with onboard robotic autonomy, allowing volumetric measurements to be made at fine scales. Special focus was given to the patchiness of phytoplankton biomass, measured as chlorophyll a (Chla), an important factor for understanding biogeochemical processes, such as primary productivity, in the coastal ocean. During multiple field tests in Runde, Norway, the method was successfully used to identify, map, and track the subsurface chlorophyll a maxima (SCM). Results show that the algorithm was able to estimate the SCM volumetrically, enabling the AUV to track the maximum concentration depth within the volume. These data were subsequently verified and supplemented with remote sensing, time series from a buoy and ship-based measurements from a fast repetition rate fluorometer (FRRf), particle imaging systems, as well as discrete water samples, covering both the large and small scales of the microbial community shaped by coastal dynamics. By bringing together diverse methods from statistics, autonomous control, imaging, and oceanography, the work offers an interdisciplinary perspective in robotic observation of our changing oceans.


INTRODUCTION
Processes controlling the growth and accumulation of phytoplankton are central to nutrient, carbon, and energy cycling; these provide the foundations for marine and human food webs (1,2). Despite their ubiquitous presence across the upper water column in the global ocean, phytoplankton presents a remarkable heterogeneity in spatial distribution (we will refer to this heterogeneity as patchiness) that is observable in microscale [subcentimeter, (3)], kilometer scale (4), and up to more than 600 km (5). Sampling at such a wide range of spatial scales is both challenging and interesting, because patchiness has been shown to strongly influence ecosystem stability (6), diversity (7), and productivity (8). Furthermore, measuring phytoplankton is complicated because observations have their basis in the proxy measurement of peak absorption of light at 675 nm, originating from photosynthetic chlorophyll a (Chla) inside planktonic cells.
With the growing capabilities of robotic vehicles, adaptive coverage of the water column is increasingly viable, although they continue to be challenging to operate in the harsh confines of the water column and the dynamic coastal ocean. Motivated by this, this work describes the use of propelled robotic platforms, known as autonomous underwater vehicles (AUVs), for optimizing mapping of water-column features in three dimensions. Such volumetric surveys of the water column are of substantial scientific interest because they provide better insight for understanding and measuring factors such as primary productivity, a vital indicator for ocean health.
To do so, we concentrated on volumetric sampling of subsurface chlorophyll a maxima (SCM) using adaptive sampling methods. In particular, the goal was to develop adaptive behavior for AUVs that can enable increased sampling resolution of water-column processes in three dimensions by focusing sampling efforts in regions of interest. We limited the survey to a bounding box with a certain size, introducing both spatial and temporal constraints, similar to the enclosure criterion used in (9). Shortening the survey time means that the temporal effects will be minimal, and the error contributions from currents, navigation drift, and other time-dependent effects will be bounded, allowing complex time variability to be mitigated. This then transforms into a purely spatial problem of finding the relevant sampling locations within the surveyed volume, given a scientific context (for our case, the maximum Chla concentration depths). Assuming that only limited previous information is available, obtaining an a priori estimate of the SCM feature has substantial benefits, contrary to conducting an adaptive survey directly using interior [inside the three-dimensional (3D) volume] measurements. This direct approach is ineffective due to two main reasons. First, the AUV will have no forward-looking data to guide adaptation, resulting in coverage that is subject to a "hit-or-miss" trade-off. Given the limited volume of the survey, the adaptation distance available for such an approach is not present. Second, a key driver is the need to survey the most interesting aspects of the feature as fast and accurately as possible to reduce time-dependent effects. Building up an accurate estimate of the interior will be central to achieving this goal, because both the probability of losing track of the feature and the survey time can be reduced.
The method suggested in this work is designed to adhere to these ideas by separating the survey into a dedicated exploration and exploitation phase as follows: The AUV first starts by covering the sides of the bounding box volume in an initial phase dubbed MODE 1, followed by an adaptive survey of the interior volume, dubbed MODE 2, capitalizing on the learned information captured in a Gaussian process (GP) model, during MODE 1 (full details are given in Materials and Methods). Thus, the exploitation phase is a single-calculated maneuver based on the available information, allowing the AUV to survey a region of interest in the water column in a fast and focused manner, toward attaining a quasi-synoptic snapshot of the feature.
Our experiments were conducted in Runde, Norway, in June 2017 (see Fig. 1) with the algorithm embedded on an AUV to detect and to map the SCM, providing a 3D map of the Chla concentration, denoted as [Chla]. Note that we specifically define concentration with [Chla]. This is different from "Chla," which is an abbreviation for chlorophyll a. The measurements were complemented with photophysiological data derived from a fast repetition rate fluorometer (FRRf) and two particle imaging systems (SINTEF) (10) with different magnifications. These were profiled from a winch, together with a CTD (conductivity, temperature, and depth) sensor and water samples taken from the research vessel R/V Gunnerus at stations around Runde. An overview of the involved systems, scales, and platforms is given in fig. S2.
In situ measurements of the size, concentration, and type of particulate material in the water column are critical to resolve spatial distributions of the particles in question while minimizing disruption to them (11). A wide range of techniques could be used for such measurements: optical backscatter, fluorescence, laser diffraction (12), digital holography (13), and telecentric transmittance imaging (10,14), each providing specific information about suspended particulate material. However, it is not possible for a single instrument to obtain measurements of the entire size range of suspended material present in natural seawater, making it necessary to combine multiple techniques to monitor the full range of size and types. By combining measurements, a detailed picture of the water column can be rendered with both biometric and volumetric information, providing fine detail for resolving the in situ features while reducing errors related to aliasing and patchiness.

Related work
Methods for tracking and sampling biomass distributions using autonomous platforms are still in their infancy. However, methodological studies are regularly reported, and several experiments have been conducted in the recent past. Moline et al. (15) introduced AUV-based sampling for ecophysiological work. Focusing on harmful algal blooms (HABs), Robbins et al. (16) combined AUV surveys with remote sensing toward attaining 3D information about bloom conditions. Similarly, Zhang et al. (17) used a longendurance AUV to sample the peak Chla layer along horizontal cross sections for HAB monitoring using onboard gene processing equipment. Ryan et al. (18) used hidden Markov models to enable online estimation of water-column properties, enabling water samples to be triggered under certain upwelling conditions in Monterey Bay. The same region was explored in (9), where a GPS-tracked Lagrangian drifter was used in combination with AUVs for coordinated sampling of dynamic oceanic features, trying to follow the sides of a volume of water as it moved with subsurface currents. More decision-based efforts in sampling were presented in (19), where GPs were used in an online optimal choice approach to adaptively decide where to trigger water samples in the water column based on biomass concentration. In (20), cooperative tracking of phytoplankton blooms using buoyancy-driven gliders driven by ocean model predictions was articulated. This was developed further in (21), where the sampling resolution was modified by adjusting the undulation angle of a glider while accounting for ocean currents, with the aim to characterize and sample blooms. Graham et al. (22) discussed problems related to estimating spatiotemporal correlation in time-varying fields and compared Euclidean and Lagrangian approaches using AUV and drifter information to perform spatial interpolation. Feature tracking of patches and plumes was also discussed in (23), where a plan-based policy was learned, tested, and evaluated using simulated patches of Chla and then sampled in 2D in the open ocean. In (24), the combined use of satellite measurement, AUV, and research vessel data was used to describe advection of phytoplankton cells from open waters transported under the sea ice in the Arctic linking Chla presence, photosynthesis, and biodiversity. Spatial and temporal decorrelation scales in phytoplankton bloom magnitude were reviewed in (25), focusing on the U.S. West Coast. Sahlin et al. (26) gave an overview of different 3D interpolation strategies for oceanic data, with the conclusion that kriging is the preferred method for studying the marine pelagic environment.
Our work presents a practical real-world implementation of a sampling approach for the upper water column, capable of adapting to distributions in three dimensions. Such uncertain environments require a trade-off between exploration and exploitation; the proposed method presents a solution by splitting the survey into two dedicated phases within a limited volume, bounding time-variability effects. Under these conditions, the AUV can map the region of interest in the water column in an efficient quasi-synoptic manner. We demonstrated the method in a field experiment together with additional sampling on board a ship, showing the potential of collaborative mapping of planktonic distributions in a dynamic coastal environment.
Sampling Chla in the water column Planktonic patchiness can be a result of complex interplay of physical (stirring and mixing), chemical (nutrient availability, CO 2 , O 2 , and trace A B metals), and biological (growth, movement, and mortality) processes, which can rapidly change in marine ecosystems and vary across multiple spatiotemporal scales (27). These interacting processes create a broad spectrum of patchiness, which can occur in small vertical layers [SCM layers (SCMLs)] (28) or in bigger plumes (29,30) that can only be partially visible in remote-sensing imagery (Fig. 1). Phenomena such as eddies can interrupt these layers by upwelling or downwelling and potentially create deep Chla maxima layers (DCMs) (31). Details about the difference between SCMLs and DCMs can be found in (32).
Knowing where and when to sample is vital to accurately resolve processes. Furthermore, data collection will, in some cases (especially Lagrangian buoys and floats), itself be influenced by the same forces, such as currents and winds, in addition to the physical limitations of the platform (speed, battery, depth rating, etc.) affecting the spatiotemporal coverage. Capturing synoptic information is therefore a major challenge. Remotesensing data can only obtain coverage at the very surface, omitting valuable information in the depth dimension and, consequently, features such as the SCM. Furthermore, remote-sensing products rarely reach the subkilometer resolution, which is vital for evaluation of oceanographic phenomena that control coastal primary productivity at fine scales (33).

Chla patchiness, distributions, and interpretation
The advent of fluorometry as a proxy for measuring Chla fluorescence, in addition to the progress of remote-sensing approaches, has provided a significant advance in the study of plankton distributions (34,35) providing continuous measurements, where Chla is used as a common "currency" for biomass estimation. Biomass is a broad and practical term used to describe the amount of living material in the water column. We use the concentration of Chla as an indicator of the phytoplankton biomass (organic carbon, expressed as milligrams per cubic meter); hence, the terms [Chla] and phytoplankton biomass, or just biomass, are used interchangeably in this text. Although Chla is used as an indicator of phytoplankton biomass worldwide, variability of cellular Chla content is expected to occur as a function of light intensity, species composition, and nutrient availability (36,37).
The use of Chla fluorescence measured in vivo, i.e., in live cells measured in situ, as a proxy of phytoplankton biomass and patchiness, needs to be interpreted carefully. When phytoplankton cells are exposed to high light levels (for instance, those usually found at the surface ocean, particularly near mid-day during summer), some of the excess of energy absorbed is dissipated as heat (38). Part of this is related to a photoprotective process, called nonphotochemical quenching (NPQ), that occurs as a rapid response to high light and results in the reduction of fluorescence emission at the surface (39). Thus, a subsurface Chla fluorescence maximum observed in a vertical profile may not necessarily reflect biomass accumulation; rather, it can be a spurious representation of patchiness. To circumvent this challenge, we compared in vivo fluorescence Chla collected in situ from the AUVs and FRRf measurements with in vitro Chla, which is derived from filtered water samples and methanol extraction to check whether a subsurface Chla maximum is present at the sampled locations.
At first glance, the complex spatial structure of the phytoplankton distribution occurs as a result of biological and physical processes and interactions in the water column. The spatial distribution of phytoplankton biomass, its intensity, morphology, and scale dependence [such as the open ocean (10 to 100 km, mesoscale) or coastal zones (10 m to 1 km, submesoscale)] is substantially driven by processes, such as turbulent advection, upwelling, convergence, and vertical mixing (27,40). For example, local processes such as upwelling zones can bring deep water nutrients to the surface layer and nurture phytoplankton, creating regional hot spots (with high biomass concentration), at scales ranging from 5 to 10 km (41) or even ≤1 km for complex coastal zones (42). In the open ocean, the same aggregation can range from 70 to 140 km (spatial correlation) in the horizontal plane; vertically, persistent upper water-column stratification may lead to a subsurface maxima, where phytoplankton is concentrated in the bottom of the pycnocline (density gradient) to better use light and nutrients (43). Vertical correlation is much weaker due to stratification effects (26). Figure 1B shows the surface distribution of [Chla] for 2 May 2017 along the Norwegian coast, depicting the spatial patterns at the mesoscale.
Regardless of the physical processes taking place, spatial heterogeneity must be "seeded" from external factors that trigger phytoplankton growth, such as nutrient availability and light levels (4), in addition to intrinsic evolutionary and physiological traits of the species (e.g., buoyancy, vertical migration, and daily cycles of photosynthesis) (3), reinforcing the role of physiology on shaping phytoplankton patchiness. Many fluorometers [such as pulse amplitude modulated (PAM) and FRRf] offer the opportunity to study the health of phytoplankton communities, which allow examination of the physiological state of a community and, when combined with reconstructions from AUV surveys, are powerful assets to investigate the spatial structure of phytoplankton distributions (24).
Zooplankton grazing is another biological process that accounts for patchiness as it relates to phytoplankton mortality. Although extremely important, it is difficult to study the effect of such grazing on patchiness due to the nonlinear nature of biological interactions and the scarcity of observations that constrain models (4). Like phytoplankton, zooplankton themselves are patchily distributed, either occurring in swarms or vertically migrating in the water column, which complicates the interpretation of abundance due to undersamplingthat is, being below the Nyquist sampling rate, so that the dynamics of a phenomena is not captured at a level necessary for reconstructing the original distribution (aliasing). The development of acoustic sensors, such as Acoustic Doppler Velocity Profiler and Acoustic Zooplankton Fish Profiler, and particle counters capable of collecting continuous zooplankton counts provides an opportunity to assess abundance and the phytoplankton influence.

Experimental setting
Our study site is in the coastal zone of the southern Norwegian Sea. In this region, the Norwegian Atlantic Current approaches the coast bringing saline Atlantic water onto the narrow continental shelf (44). The Norwegian Coastal Current, carrying fresh water from the outflow of the Baltic and receiving additional fresh water from land runoff as it flows northward along the coast, also enters the shelf in the Møre region (44,45). Both the volume transport of the Norwegian Coastal Current and the lateral location of the interface between Atlantic and coastal water vary seasonally and are strongly affected by wind (46,47); there is also considerable mesoscale variability in the coastal current system (47,48). Mesoscale eddies play an essential role in ocean dynamics and significantly affect ocean biology and biogeochemistry (49). Although mesoscale eddies, from a few to a few hundred kilometers, are known to have a strong effect on plankton distributions, community structure, and diversity (49), smaller instabilities can locally enhance vertical mixing (50). Such submesoscale processes, on spatial scales of 0.1 to 10 km, can enhance both phytoplankton productivity and its variance on time scales of a few days (50).
Repeated sampled temperature and salinity profiles from Skinnbrokleia (7 km southeast of our AUV surveys) were used to place the AUV observations in a local and seasonal setting. Meteorological and oceanographic (metocean) time series from a stationary buoy at Breisundet (14 km northeast of the AUV surveys) were also used. The buoy, placed in 345-m water depth in an inlet between two islands, measures surface waves, wind, current velocity, salinity, and temperature at selected depths.
In the weeks preceding the field campaign, the hydrographic conditions in the study area were typical of early summer (51,52). Upper layer salinity was low, and a shallow seasonal thermocline had formed in the upper 10 m (seen in data from a monitoring station on 8 June, 16 km southeast from the survey site). The beginning of the cruise, however, coincided with the onset of stronger winds, around 10 m/s with gusts up to 16 m/s, that prevailed throughout the first half of that week (Fig. 2B). The fresh breeze was accompanied by increasing wave height, up to 2 m (as measured at the stationary metocean buoy in Breisundet). Wind stirring during this event caused a deepening of the mixed layer. Time series of temperature and salinity from the metocean buoy revealed that the mixing started immediately with the onset of stronger winds on 17 June and progressed during the following days until eventually temperature and salinity were homogenized down to the 40-m depth level ( Fig. 2A).
The AUV and ship-based CTD profiles taken on the same day showed a 30-m-deep well-mixed layer (Fig. 3). The water-column structure was similar in the AUV survey area and at the ship-based CTD stations some 400 m away. In the thermocline between depths of 30 and 40 m, we found some features that were different, such as local temperature inversions, and that were seen at one site (CTD stations) but not at the other (AUV surveys). The different measurement methods lead to small variations between the profiles: The AUV captured the surface layer from the top centimeter down, whereas the CTD profile excluded the uppermost meter. On the other hand, the CTD profiles reached deeper than the AUV's maximum dive depth of 50 m. The many profiles from the AUV survey form an envelope that indicates the small-scale spatial variability of temperature and salinity. The difference in water mass properties between the CTD sites and the AUV survey is comparable with the variability within the AUV survey area, which, in turn, is comparable in size with the distance between CTD stations and the AUV survey (Fig. 1A).

GP model of the distribution of Chla
We used a spatial statistical model in the form of a GP to model the SCM. The GP formulation is widely used for dealing with spatial applications in the ocean (22,53). Its popularity is largely due to the simple yet flexible modeling formulation, which allows consistent and efficient data conditioning. In contrast to a full numerical ocean model, a GP can operate with limited computational resources on board an AUV and assimilate measurements on the fly.
Measurements were coherently assimilated to render a more detailed representation of the SCM as the survey progressed. The SCM is predicted by a GP in 3D space (latitude, longitude, and depth), by assimilating the identified SCM depth into a 2D grid, of which the detailed steps are explained in Materials and Methods. For environmental applications, a GP is typically discretized to a grid map with a certain resolution describing variation in space. This map is, in essence, a collection of random variables (on the grid) with a multivariate normal probability density function, where the spatial dependence is captured by a covariance model.  important, because this will affect how the measurements are interpolated. In our experiments, the effective correlation range was set to f~1000 m based on the east-west variogram results from (42), investigating the Chla scale dependency for coastal zones using highresolution remote sensing [notably, this implies having correlation across the smallest survey area ( ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 700 2 þ 700 2 p ¼ 990 m)]. Here, both larger (26) and smaller (55) correlation ranges can be defined, depending on the scale of heterogeneity, type of marine system being studied, and previous insight of the water body and ambient conditions. As measurements are assimilated, the 3D mean surface is adjusted toward the derived Chla maximum depth, and prediction variances are reduced (details are provided in Materials and Methods). The SCM peak depth is found for each vertical crossing as the AUV undulates via a "yo-yo" maneuver in the water column, where d t is the AUV location at time t with the corresponding SCM peak. The measurement model for time t is given by where y t;d t is a length m t;dt vector of the survey line observation and the m t;d t Â n matrix G t;d t contains "1" entries only at the designated d t indices and 0 otherwise. The error term v t e N ð0; R t;d t Þ represents measurement noise. The covariance matrix R t;dt is set to a constant diagonal matrix (56), and there is no dependence of measurement error terms over assimilation time. Under Gaussian linear modeling assumptions, the sequential updating of survey data leads to the Gaussian conditional distribution pðxjy 1;d1 ; …; y t;dt Þ, defined by conditional mean m t;d t and covariance P t;d t . These equations are recursive over the data-gathering steps: where the mean and covariance matrix are initialized by m 0;d0 ¼ m and P 0;d 0 ¼ ∑.

Experimental setup
The Alongside AUV measurements, both FRRf and two SilCam camera systems were used to resolve details of the water column, as well as allowing cross-reference and verification of AUV measurements and behavior in addition to the algorithm embedded onboard. These systems were deployed from a winch, together with a CTD. The FRRf measures irradiance and [Chla], in addition to photophysiological parameters of the phytoplankton community. The SilCam systems imaged a particle size range spanning 28 mm to 4 cm in equivalent diameter, enabling characterization of particle size over a range covering three orders of magnitude in diameter. For accurate fluorescence-based Chla determination and calibration of the FRRf, water samples were collected every 10 m in the upper 60 m on the upward CTD casts using 5-liter Niskin bottles mounted on a rosette frame. Water samples were immediately filtered onto 25-mm glass fiber filters (from GF/F Whatman Inc.) and stored in a −20°C freezer until analyzed in the laboratory 1 month later.
[Chla] was extracted in 100% methanol for about 4 hours at 10°C and fluorometrically determined using a Turner Designs fluorometer. This sampling activity was conducted from the research vessel R/V Gunnerus. The survey map and sampling stations are shown in Fig. 1A.

AUV surveys
The adaptive AUV surveys took~1.5 hours to complete (survey 2 covered a larger volume), with a mean vehicle speed of 1.6 m/s over ground; details are provided in Table 1. To follow the GP model estimate of the SCM 3D surface, we explored two strategies: (i) using fixed depth (MODE 1.fixed) and (ii) using an undulating/yo-yo motion, ±2.5 m across the 3D surface (MODE 2.undulating). Although finer discretization is possible, given the uncertainty due to current speed and direction, as well as natural stratification, the SCM depth estimates from the GP surface (~23-m resolution) were discretized into 200-m segments for the AUV to follow. A detailed analysis of the approach is presented focusing on survey 1, using fixed depth control, and comparing that with survey 3, tracking the 3D surface using an undulating/yo-yo pattern. Results from survey 2 are shown in fig. S3.
In Fig. 3B, the first~4 km is used to cover the sides of the volume (MODE 1), gathering an estimate of the SCM, before MODE 2 tracks this estimate by either fixed or undulating/yo-yo behavior. MODE 2, in turn, trails an hourglass pattern in the interior of the volume for the next 4 km (see Materials and Methods for more details). The second half of Fig. 3B shows the AUV adjusting its depth according to the SCM estimate, with periodic surfacing events when reaching the corner points of the survey box. Both survey 1 (Figs. 3A and 4) and survey 3 (Fig. 5A) confirm that the [Chla] concentrations were higher during the tracking phase of MODE 2, with the AUV sampling between the 20and 30-m depth range. These depths correspond to the SCM found from analysis of the water samples taken from the research vessel and shown together in Fig. 6 with FRRf data. The temperature and salinity data in Fig. 3C suggest that the water column was well mixed down to 30 m, where an abrupt gradient is present. This interface influences the distribution of biomass in the water column (58).
Because the first phase (MODE 1) took~40 min to complete, some of the estimates were likely affected by current velocity and direction. It is therefore important to verify that the depth adjustments result in an increase of [Chla]. Elevated levels of [Chla] can be seen throughout the second phase of the missions, but to verify this further, in Figs. 3 (A and  B) and 4 (B and C), an asterisk shows the first part of the SCM tracking and the successful depth adjustment that results in elevated levels of [Chla]. The first part is chosen, because the earliest observations have the highest probability of being affected, due to their high lag between the initial measurements and revisitation, affected by current.
Furthermore, by comparing the volumetric estimate before (Fig. 4A) and after (Fig. 4, B and C) the interior survey, we can see the distribution change in light of the new measurements. The isosurfaces (marked with black lines in Fig. 4) changed, resolving the interior information with finer detail. The change is greatest within the volume, which is expected, because these locations are unsurveyed, previously estimated on the basis of the data from the sides. The main features and distribution are largely similar before and after the initial survey. This coherence supports the credibility of the design premise, limiting the survey size to effectively handle time dynamics, as well as the shared information between the interior and boundary used to focus sampling.
Comparing the results with survey 3 using the undulating/yo-yo depth tracking shows that survey 3 also has elevated concentrations while tracking the SCM. Because the undulating/yo-yo approach is less sensitive to errors in the depth estimate while providing higher sampling density at the SCM, it is more practical for mapping the structure of the SCM but only subject to a certain current regime; with strong currents, a fixed depth approach would likely be preferable because of shorter survey times. Figure 6 shows a vertical profile of the in vitro (extracted) and the in vivo (measured in situ) [Chla] derived from the FRRf deployed at stations 4 and 5 and from AUV survey 1 indicated on the map in Fig. 1A; all FRRf surveys are shown in Table 2. Details on intercalibration can be found in (24). In the in situ assessments (FRRf and AUV), [Chla] is higher  A quenching (reduction) of the fluorescence signal in live cells (in situ) is commonly found in surface waters during summer, when phytoplankton are exposed to high light levels. This occurs to protect the cells' photosystems [a process called NPQ (38)], by converting part of the excess energy to heat. Results from this study show that the in vitro [Chla] presented approximately the same vertical pattern as the in situ measurements, with a subsurface maximum at depths of 20 to 30 m. This suggests that NPQ processes were not evident in this study, possibly because of cloudy conditions. Clouds reduce the amount of light that reaches the surface water; consequently, such conditions do not interfere as much with the in situ fluorescence measurements derived from the sensors.

SilCam observations
In situ [Chla] derived from the FRRf and the AUV provided relevant information regarding phytoplankton patchiness, but there may exist some linkages to biological activity (e.g., grazing) at higher trophic levels and thus larger particle sizes. Results from the two SilCam particle imaging systems provided in situ images of suspended material, allowing insight into the presence and distribution of larger phytoplankton and zooplankton, in addition to other material in the water column, such as marine snow and suspended sediments. Figure 7 shows particle number distribution and sample images recorded within the diameter range of 30 to 10,000 mm at stations 4 and 5 using the two SilCams with different magnifications. With an average exponent (m) of 3 (shown by the dashed lines in Fig. 7), the size distribution follows approximately that of a Junge distribution (59): where n(D) is the number of particles of diameter D, k is a constant that scales according to the particle concentration, and m is the slope of the distribution. In the simplest form, the Junge distribution is an approximation of the particle size distribution, made up of a series of multiple log-normal distributions with varying median sizes that arise from subpopulations of particulate material such as specific types of phytoplankton, zooplankton, or marine snow (60). In Fig. 7, the observed distributions show many points of variation from this average Junge slope, with one peak at just over 100 mm that corresponds to a high abundance of the dinoflagellate genus Ceratium (indicated by the red shaded region), which can be seen in the middle column of  West and between S2 and S1 Low the figure. The copepod genus Calanus is also highly abundant (indicated by the blue shaded region) and can be seen in the rightmost column of the figure. In the montages of particle images (right two columns of Fig. 7), individual particle images are packaged into the plot by randomly sampling from the recorded size distribution so as to maintain the ratio of large to small particles, albeit without being representative of the true concentration of particles, which was more sparse in situ. Most particulate material is present within the top 10 m, and it is where relative abundance of Ceratium appears greatest because the 100-mm peak in the particle number distribution is most prominent. Ceratium is still present in approximately the same concentration between 10 and 30 m but is surrounded by a higher concentration of other material also present within this depth range. This is evidenced by the reduced dominance of the 100-mm peak. Below the thermocline at 30 m, the Ceratium peak in the number distribution is not discernible, and there is less total material in suspension. With such a high dominance of biological material in the water, the reduced concentrations below 30 m, observed by the SilCam, agree well with the reduced [Chla]. The highest load of particulates observed by particle imaging was within the 10-to 30-m depth range, which also corresponds to the SCM. However, the small difference in the concentrations of material at 100 mm between the 0.5-to 10-m and 10-to 30-m depth ranges suggests that other smaller phytoplankton species (possibly not captured by the imaging system), in addition to Ceratium, are contributing to the SCM.

DISCUSSION
We have presented methods for autonomously mapping spatial heterogeneity of phytoplankton biomass in 3D using AUVs by combining GP models and robotic sampling. The method is shown to successfully estimate and track layers of high Chla concentration, focusing sampling efforts and increasing resolution along important features such as the SCM. We used spatial modeling and interpolation to reconstruct the distribution in 3D. This volumetric estimate was then augmented with in situ images of suspended particles, fluorometry, and discrete water samples taken co-temporally from a research vessel. Comparison of infield data shows correspondence between AUV data and behavior, providing a broad and extensive perspective of the pelagic activity. The results demonstrate the complexity of conducting interdisciplinary coastal ecology and support unification of marine data sources to achieve a detailed environmental picture of the water column.

Limitations and future steps
The largest source of sampling uncertainty comes from the effects of currents (speed and direction), making the observations time-dependent. The AUV used in this study is fast enough to find and make use of spatial structure, but not so fast as that it can capture the full spatial field at any time, because an AUV is not a synoptic measurement platform. Ideally, the AUV should therefore try to "stay" with the same water mass, working in a Lagrangian frame of reference to reduce these effects. Lagrangian correction is, however, nontrivial, and, unless a good proxy measurement of advection can be provided (such as from a surface drifter), there is limited value in adding to experimental complexity [see (22) for a more detailed discussion]. A simple and effective measure to account for this is to limit the method/survey area to subkilometer size, setting a bound on these uncertainties similar to the enclosure criterion used in (9). Use of this strategy allows us to use directed Euclidean sampling of the interior volume while increasing sampling resolution and determining SCM variability. This can be seen by looking at the recorded data and comparing the estimated Chla distribution before and after the initial survey, as in Fig. 4, where only a minor current influence was discernible (except at unobserved locations).
Furthermore, the approach is designed with two modes: MODE 1 that is purely exploratory and MODE 2 that is exploitative. Consequently, the most interesting part of the volume is surveyed last. This is important because previous data from MODE 1 allow sampling (the interior) to proceed with higher survey speed and reduced probability of losing track of the feature. The number of yo-yo envelopes is platform dependent, and the spacing, depth, and number of envelopes could be used to further optimize the time spent finding an adequate estimate. In addition, forsaking some surfacing events could help reduce survey time, if more accurate navigation is possible. At the cost of a greatly increased operational complexity, additional assets such as another AUV could be used to increase coverage capacity and sampling density, should sharing of information between assets be available [e.g., (55)].
The AUV does not update the survey depths based on measurements inside the volume due to the previously mentioned factors such as hit-or-miss consequences and time-variability effects. Using the measurements from the interior (within the defined volume) survey can be considered in future work, although the value of this information is limited because of the small survey area and correlation scales present in the coastal ocean. Thus, making use of these data can be challenging. One possible use would, for instance, be to improve the undulating depth from a fixed oscillation to a data-adjusted offset. It is also the case that variability in the ocean is higher along the vertical axis, motivating the use of such capabilities. However, this variability depends on a host of factors and can be hard to estimate correctly without the aforementioned drawbacks. Furthermore, the results from this work show that platforms such as AUVs can make measurements adequate to resolve fine-scale features using an enclosure criterion. This is also apparent in the results presented, which shows spatial coherence throughout the mission (see Fig. 4).
The volumetric results from kriging (such as Fig. 4) are dependent on the accuracy of the correlation parameters used. The values used in this paper are based on the variogram from (42) and can only give a suggestive estimate for conducting similar experiments. Aspects of spatial correlation and use of marine sensor fusion toward description of the water column on the submesoscale require effort that is not within the scope of this work. This includes looking at the time dependence of the data, as well as using previous historical data sources (25,61).
Future integration of particle imaging instruments on board the AUV would greatly change how volumetric surveys can be conducted, enabling adaptation based on taxonomically descriptive data (targeting HABs or other materials such as suspended plastics) as well as reducing the current sensitivity to patchiness arising from ship-based winch deployment.

MATERIALS AND METHODS
The aim of the study is to obtain a fine-resolution synoptic view of a water-column feature, using a single AUV, by concentrating sampling efforts in the water column. In our case, the feature of interest is the SCM. To achieve this, we split sampling into two modes or phases: an initial phase dubbed MODE 1 that is explorative, followed by a second phase (MODE 2) that exploits the information from the first phase to sampling a feature of interest, using the data to plan the survey. By this design, the most interesting part of the volume is surveyed last, allowing the time spent mapping it at high resolution to be minimized. This is a subtle but important point because this reduces both spatial and temporal uncertainty of the most critical measurements. To avoid redundant data collection, we also split sampling of the volume between the modes. MODE 1 surveys the sides of the volume, whereas MODE 2 surveys the interior. MODE 1 uses a spatial model (GP) to estimate the SCM distribution inside the volume so that MODE 2 can follow these estimates subsequently. The following and Fig. 8 show these steps in detail: Step 1) MODE 1 Aim: Estimate the SCM inside the volume using the GP model. Behavior: Survey the sides (red square shown in Fig. 8) of the volume while doing yo-yo envelopes from 0 to 50 m in depth. Tasks: MODE 1.peak-Identify the depth of the Chla peak. MODE 1.assimilate-Assimilate SCM depths into a GP surface.
Step 2) MODE 2 Aim: Use the learned information from MODE 1 to perform a detailed survey. Behavior: Track the GP SCM depth 3D surface inside the AUV survey volume following an hourglass pattern (shown as blue lines in  MODE 2.fixed-Follow the estimate using fixed depth control. MODE 2.undulating-Follow the estimate using undulating depth control. MODE 1 comprises the initial phase, where the location and depths of the SCM are to be identified and assimilated (MODE 1. assimilate) into the GP SCM depth 3D surface. During each full yo-yo envelope (one descent and one ascent), the final SCM depth is to be found and stored as [latitude, longitude, depth, Chla concentration]. It is the task of the subroutine MODE 1.peak to do this analysis. The function will decide the final SCM depth based on the following logic: 1) Segment the concentration measurements into depth bins with average measurements. Collecting five data points in one bin, with an average depth.
2) Find the two depths corresponding to the largest (depth 1 ) and second largest (depth 2 ) concentration spikes.
3) Calculate the depth difference between these spikes. 4) If the depth difference is less than 10 m, the SCM depth is set to the depth of the largest concentration: SCM depth = depth 1 . If not, the SCM depth is set to a weighted mean, as SCM depth = depth 1 * 0.70 + depth 2 * 0.30. This is based on the fact that about 70% of normally distributed values are within 1 SD of the mean. The subsequent subroutine MODE 1.assimilate will integrate the final SCM depth into the GP model (Eq. 3) using the observation location (lat, lon → grid cell) and "lifting" or "lowering" the depth estimate at the surrounding locations/grid cells within the correlation range. These two concepts are encapsulated in Fig. 8.
Once the GP SCM depth 3D surface has been established, the AUV switches to MODE 2 and depth tracking mode, where the goal is to provide a high-resolution volumetric representation of the Chla structure by tracking this 3D surface in either MODE 2.fixed or MODE 2. undulating following a 3D path. This enables an effective and improved way of mapping detail within the SCM. The biomass patchiness can thus be evaluated and augmented with supporting ship-based samples, allowing organism community structure to be described more accurately.

Kriging in 3D
Using kriging (54), we can construct a 3D representation of the volume, conditioned on the data. The kriging surface is the interpolated prediction given as the conditional mean of the GP while providing the conditional prediction covariance in an explicit form.
Spatial correlation is an essential component of kriging, and this is incorporated by a kernel function that relies on the distance between pairs of locations, h ij = |s i − s j |; in this work, this is measured in 3D [x, y, z] coordinates. Compared with the more common 1D or 2D applications of kriging, this means a considerable increase in the size of the covariance matrices required for data conditioningassimilating all the samples from the survey is therefore not practical when doing the interpolation. In our approach, a pruned version of the assimilated 1700 samples from this volume was used, with a square size of 2700 covariance matrix (reflecting the [30 × 30 × 30] volume grid).
The interpolation time on a Dell Latitude E7540 i7 took about 1.5 min with these parameters, resulting in figures such as Fig. 4. To capture the difference between horizontal and vertical correlation dis-tance, we artificially exaggerated the z dimension by a factor of about 20 based on (26). The horizontal correlation was, as previously discussed, based on the east-west variogram from (42). Although there are methods for handling memory and computational time challenges for 3D GPs, our focus was on an applied study using MODE 1 to define a region of interest.