Dynamic stochastic modeling for adaptive sampling of environmental variables using an AUV

Discharge of mine tailings significantly impacts the ecological status of the sea. Methods to efficiently monitor the extent of dispersion is essential to protect sensitive areas. By combining underwater robotic sampling with ocean models, we can choose informative sampling sites and adaptively change the robot’s path based on in situ measurements to optimally map the tailings distribution near a seafill. This paper creates a stochastic spatio-temporal proxy model of dispersal dynamics using training data from complex numerical models. The proxy model consists of a spatio-temporal Gaussian process model based on an advection–diffusion stochastic partial differential equation. Informative sampling sites are chosen based on predictions from the proxy model using an objective function favoring areas with high uncertainty and high expected tailings concentrations. A simulation study and data from real-life experiments are presented.


Introduction
The world's oceans are heavily impacted by human activities, such as pollution and overfishing (Halpern et al., 2008). One example is the harmful effects of mine tailings, drill cuttings, and mud discharge. Along the Norwegian continental shelf, there are sensitive species like cold-water coral reefs and sponges which are negatively impacted by such deposition (Trannum et al., 2010), but still, there exist several active mining seafills and new seafills are in the planning phase. Improved methods for monitoring the extent of the contamination from seafills are important to ensure that sensitive areas are protected. In this paper, we use a seafill located in Fraenfjorden (see Fig. 1) as a case study where we validate our proposed method for monitoring the tailings distribution close to the outlet.
When mapping the tailings distribution, models of ocean dynamics can assist us with helpful information. However, building such models is challenging because of the oceans' large-scale nonlinear processes and high spatiotemporal variability. Existing models use numerical methods to improve their accuracy (Griffies et al., 2000) continuously, but they still have errors due to simplifications, limitations in numerical implementation, uncertain parameters, and initial and boundary conditions. Thus, data assimilation using in situ measurements can help improve the model accuracy. Measurements are commonly obtained using remote sensing, ships, or buoys, but this data is expensive to collect, making it hard to observe the environment in detail. Recently, robotic vehicles like autonomous underwater vehicles (AUVs) have become more robust and less expensive, and by equipping them with sensors, one can collect measurements (sample) more efficiently. Still, the possible observations are sparsely distributed, leading to significant gaps in our understanding of the ocean (Stewart, 2008). This measurement limitation prompts the need for intelligent and targeted observation designs, and an important question is when and where to sample.
This paper addresses this question and suggests an adaptive sampling method utilizing prior data from ocean models and real-time in situ measurements from an AUV. Because the already existing ocean dynamics models have a high computational load, making them unfit for running on embedded robotic systems, an onboard low-complexity proxy model is created so that real-time adaption of the mission can be obtained. The proxy model represents the state of the ocean at the time and can be updated when new information is available. We use a spatio-temporal Gaussian process (GP) with a kernel describing spatial dependencies and a temporal model built using an advection-diffusion stochastic partial differential equation (SPDE). The linear SPDE model and the Gaussian measurement model are written as a state-space model, and the Kalman filter is used to update the state recursively over time. Two existing numerical models, SINMOD (ocean model) and DREAM (particle dynamics model) are used to train the compact GP model creating a proxy model of the tailing distribution. Sensor readings on the AUV can then be used to update the proxy model onboard in real-time. An objective function for waypoint selection is presented to maximize the value of information from the measurements. This function ensures that the area is explored by choosing locations assumed to be information-rich and also considers the travel length limitations of the AUV, leading to an adaptive sampling strategy.
In our case study, the method is tested in simulation using the Fraenfjorden area (Norway), a fjord containing a seafill for mine tailings. The fjord can be seen in Fig. 1, where the location is marked with a red circle. We compare the proposed method to four simpler strategies in the simulation study. As a proof of concept, real-life experiments are carried out in the same area.
The paper is organized as follows: Sect. 2 presents background information on the Fraenfjorden area, the ocean models, and adaptive sampling. In Sect. 3, the development of the proxy model is presented. Section 4 presents a combined state space model explaining how the data assimilation is done before the sampling strategies are discussed in Sect. 5. Section 6 contains simulation results comparing our proposed method to four simpler strategies, and results from experimental field trials are presented in Sect. 7. Section 8 discusses the results before Sect. 9 gives a conclusion and future work.

Related work and contribution
With the improved technology and robustness of AUVs, adaptive sampling using these robotic platforms has had increased focus lately, and a review of the existing adaptive sampling methods for AUV missions can be found in Hwang et al. (2019). Most adaptive sampling methods utilize GPs (Cressie & Wikle, 2011), which are widely used to create non-parametric and computationally efficient models. They are popular not only for adaptive sampling but used in various applications like neuroimaging (Hyun et al., 2016), machine learning (Rasmussen & Williams, 2005) and sensor placement (Krause et al., 2008). Combining GPs with robotic vehicles is amongst others explored in Zhang et al. (2012) where an AUV is used to track an upwelling front, and in Luo and Sycara (2018), which explore a method of multiple vehicles using a mixture of GP models. Foss et al. (2022) use a spatio-temporal model for the same case we are studying here, but with another objective function relying on specifying a threshold for the particle concentration.
This work presents an implementation of an adaptive sampling method for environmental sensing of mine tailings close to a seafill. We show how to utilize prior information from ocean models to create GP models. Similar approaches, also utilizing prior knowledge to learn GP models, can be found in Das et al. (2015) who use an AUV to collect samples for ex-situ analysis, selecting the sampling locations based on previous mission data and maximizing a utility function, and in Ma et al. (2018), who create a model of a spatio-temporal environment by learning and refining the hyperparameters of a GP kernel using prior ocean model data. In our approach, the ocean model data are used to develop a non-stationary correlation kernel as in Berget et al. (2018) and Fossum et al. (2018) and set the parameters in a temporal model based on the stochastic advection-diffusion partial differential equation (Sigrist et al., 2015). This allows the model to capture critical dynamic phenomena such as tides. Estimated solutions of the SPDE can be written in a linear form so that the Kalman filter can be used to update the model state. An objective function based on variance reduction, hotspot detection, and AUV operation time constraints is used; see also Berget et al. (2018). The approach sequentially selects waypoints based on the current model state, and the method is tested in both simulations and field experiments near a mining seafill in a Norwegian fjord.
The main contributions of this paper are: -Propose an adaptive sampling method for environmental sensing of mine tailings. -Use ocean model data to set the parameters in a linear temporal model based on a stochastic advectiondiffusion partial differential equation. -Present results from real-life experiments using the proposed adaptive sampling method.

Background
This section presents definitions and background information on ocean models, the spatio-temporal model, and data assimilation.

The Fraenfjorden area
In Norway, several active and planned seafills have implications for the marine environment (Ramirez-Llodra et al., 2015;Morello et al., 2016). The tailings typically consist of particulates and residual chemicals from processing, and this waste will smother the sea bed close to the sea fill. In addition, some of the waste may be transported by ocean currents away from the outlet, thus affecting a larger area of the sea bed and the water column. As the case study for this paper, we focus on a seafill located in Fraenfjorden, Norway (see Fig. 1). A calcium carbonate processing plant in the factory has been operating for several decades. This processing of ore from mining causes fine-ground waste (tailings) (Kvassnes & Iversen, 2013), disposed of on the ocean floor, in a seafill close to the factory.
The outlets in Fraenfjorden are located at approximately 20 m depth, and tailings are continuously discharged with a time-varying discharge rate. A typical discharge rate from the outlets is around 400 m 3 /h. The discharge alternates between two pipes, illustrated by black lines in Fig. 1. The operational discharge for the data in the simulation study and the field experiments is on the western pipe. This outlet is marked as 2.3 in the figure.

Ocean models
Numerical oceanographic models are used in this paper to build a prior belief of the state of the ocean and to predict the movement of particles with the currents. Ocean models provide information on salinity, currents, temperature, density, and pressure and describe the state of the ocean at a given time. The models use a set of hydrodynamic and thermodynamic equations called the primitive equations, and these are commonly solved numerically using finite difference methods. This makes the models computer intensive; hence, high-resolution predictions can only be computed for limited areas, depending on the available computer resources. Input to ocean models typically includes tides, sea-level pressure, wind, heat exchange, bathymetry, and freshwater runoff. Ocean model outputs are uncertain due to limits in model resolution, limits to our knowledge of the processes resolved by the model, errors in initial values, boundary conditions and parameter values, and inaccuracies in numerical implementations (Lermusiaux et al., 2006). It is, therefore, necessary to validate the output against observations, and data assimilation can be used to correct the model state estimates. For monitoring of releases of mine tailings or other particles or dissolved material, there is a need to develop enabling technology for efficient and targeted ocean sampling. Observations from different platforms such as buoys, ship-based sampling, or robotic sampling from, for example, AUVs can be used to evaluate or, through data assimilation, improve the model estimates.
In this paper, we use data from two specific models to build the stochastic onboard model, SINMOD (ocean dynamics model) (Slagstad & McClimans, 2005;Wassmann et al., 2006) and DREAM (particle transportation model) (Rye et al., 1998(Rye et al., , 2008. Nepstad et al. (2020) described the SIN-MOD and DREAM setups for Fraenfjorden and evaluated the models against field measurements. The relation between SINMOD, DREAM, and the onboard model can be seen in Fig. 2. More details on the two ocean models are given below.

SINMOD: ocean dynamics model
SINMOD is a numerical ocean model system that connects and simulates physical and biological processes in the ocean (Slagstad & McClimans, 2005). The model is based on the Navier-Stokes equations and uses a nesting technique where high resolution models obtain their boundary conditions from larger model domains with lower resolution. SINMOD is established in configurations with horizontal resolutions ranging from 20 km to 32 m. Input to the model Fig. 2 The relation between SINMOD, DREAM, and the adaptive sampling framework. SINMOD provides ocean currents data as input to the particle transportation model DREAM. Both SINMOD and DREAM provides data used to build the prior belief of the onboard model. Ocean current predictions are used in the temporal part of the model and particle concentration predictions on the boundary are used in the updating of the model includes atmospheric forcing, freshwater runoff, and boundary conditions (temperature, salinity, tidal forcing, currents). Ocean dynamics data from SINMOD is used in this case study as input to DREAM and the onboard proxy model (see Fig. 2). Atmospheric forcing for daily forecast simulations was obtained from NOAA's Global Forecast System with 0.5 • resolution (for large-scale model domains) and from the Norwegian Meteorological Institute's MetCoOp Ensemble Prediction System with 2.5 km resolution (for local model domains).

DREAM: particle transportation model
DREAM is a Lagrangian particle transport model which can be used to simulate the behavior and fate of marine pollutants, including particulate discharges from drilling operations (Rye et al., 1998(Rye et al., , 2008. It provides time series of concentrations of released materials in the water column and deposition of these materials onto the sea floor. Input to the DREAM model includes hydrodynamic data delivered by SINMOD and information about the release (amount, rate, densities, grain size distribution). DREAM is often used as a decision support tool for managing operational discharges to the marine environment. This paper uses data from DREAM to set the initial belief for the mine tailings distribution on the onboard proxy model. This data is also used to estimate parameters in the SPDE model, including the GP kernel.

Adaptive sampling
Observing the ocean by taking measurements or samples is crucial to increasing our understanding of ocean processes. Ocean sampling can be seen as the action of sensing or measuring one or more physical properties in the ocean at a specific point in space and time. Traditional methods include mounting sensors to buoys which limits the sampling to only one location, or immersing instruments in the water from boats, commonly using a predetermined plan of sampling locations. Exploiting robotic platforms and especially AUVs can intensify ocean sampling, giving more data and better insight.
We need more sophisticated sampling strategies to exploit the benefits that robotic platforms have regarding autonomy and increased mobility. To obtain the most scientifically relevant measurements effectively, adaptive approaches need to be developed. In contrast to static/pre-scripted schemes, adaptive/data-driven strategies can react to measured conditions, having access to both prior and in situ information; selecting sampling locations depending on past observations taken during exploration.
The benefits of adaptive sampling will depend on the situation. There are numerous studies on the value of additional (sequential) sampling efforts in different contexts, say clinical trials in medical research (Jennison & Turnbull, 1999) or geoscience data (Eidsvik et al., 2015). In our case, there would be limited to gain by adaptation if a priori information makes it possible to pre-script a high-value AUV plan before deployment. However, if there is high prior variability, the ability to adapt the initial plan one way or the other, based on onboard measurements, would bring added value by guiding the AUV to more relevant places.

Modeling
This section describes the advection-diffusion SPDE model and its associated spatio-temporal GP model representation. We also outline parameters in the model.

Spatio-temporal onboard model
To select optimal sampling locations in a dynamic domain, information about the spatial conditions at all times is essential. Since the computational demands using ocean models such as SINMOD and DREAM are too high for real-time onboard modeling, we suggest a stochastic spatio-temporal proxy model based on an advection-diffusion SPDE. When this SPDE is discretized, it represents a spatio-temporal GP model. Apart from having a smaller computational footprint, GPs are conventional tools for statistical modeling of spatial data and have been widely adopted in oceanographic applications . In this paper, the proxy model is used to approximate the distribution of the particle concentration near the seafill in Fraenfjorden.
There exist various types of spatio-temporal GP models, see e.g. Chapter 6 of Cressie and Wikle (2011) for a review. For statisticians, it has been important to understand the spatio-temporal covariance function of such processes, which is crucial in spatio-temporal Kriging. In the simplest case, the covariance function is separable in space and time, but there has also been much work on non-separable covariance functions (Gneiting, 2002). Stochastic differential and difference equations in the spatial domain fit naturally into state-space models, which enable spatio-temporal Kalman filtering using in-situ observations. The spatial version of the autoregressive process is one such example of a linear difference equation. So are integro-difference equations which weight spatial variables from the previous time step to predict the spatial variables at the current time step (Storvik et al., 2002). Recently, extensions of the advection-diffusion equation have been developed (Sigrist et al., 2015;Richardson, 2017). These models provide more interpretability in the spatial weighting as well as opportunities for embedding time-dependent drift. We use one of these models in the current paper. The advection-diffusion equation is also endorsed in the physical modeling community, and in our setting the advection can be specified from the tidal currents and the diffusion coefficients can be set from the ocean models SIN-MOD and DREAM.
We consider a real-valued stochastic process {X (t, s), t ≥ 0, s ∈ D}, representing the uncertain particle concentration on log scale in (east,north) location s = (s e , s n ) ∈ D ⊂ R 2 at time t.
The spatio-temporal process of transport of a medium in a fluid can be expressed by the advection-diffusion SPDE (Sigrist et al., 2015) with c(t, s) = (c e (t, s), c n (t, s)) T is the drift vector field for advection, D the diffusion matrix, (t, s) a noise term, ∇ the gradient operator ∂ ∂s e , ∂ ∂s n and λ ∈ [−1, 0] a damping constant controlling the autoregressive relationship between state vectors (Richardson, 2017). The noise term (t, s) is assumed to be a zero-mean Gaussian noise term that is uncorrelated in the temporal domain but spatially correlated, Corr( (t, s), (t, s )) = R(|s − s |), where R() represents a correlation function for distance input. The source term B(t, s) describes external information, e.g. boundary information.
This paper assumes horizontal, isotropic diffusion leading to a constant diffusion parameter D. We denote the 2-dimensional drift vector as c(t, s) = (c t e , c t n ) T at time step t. Then we get The SPDE is discretized using the upwind differencing scheme to find an approximate solution. We define the length of a time step as dt. Forward differences are used in time, giving We define a regular grid D s with a total of N grid points with resolution ds e in the east direction and ds n in the north direction. The set of spatial locations (east, north) is then For the first derivative in space, we alternate between forward and backward differences based on the direction of the drift vector c t .
In the eastward direction, we have where the upper sign is used for c t e < 0 and the lower is used when c t e ≥ 0. Similarly, in the northward direction, we have where the upper version is used for c t n < 0 and the lower is used when c t n ≥ 0. The second derivatives are discretized as follows, By scaling the time discretization, a change in time can be denoted t +1 instead of t +dt. We have {t = t 0 + j ×dt; j = 0, 1, 2 . . . , T }. When inserting all the discretized terms into eq. (1), the scheme can be written as a linear process on the form where X t is the size N ×1 vector of particle concentrations at time step t for all spatial locations, A t is the N × N propagator matrix including information about advection, diffusion and damping, B t is a matrix adding external information, in our case boundary conditions, and W t is a positive definite variance-covariance matrix for the noise term t .
The particle concentrations at time t = 0 are assumed to be Gaussian distributed with mean vector For the boundary conditions in the SPDE, we use Dirichlet BCs on the east side of our survey area. This can be considered as a constant forcing along the east side and is chosen because of the constant flow of mine tailings that are let out on this side of the rectangle. We assume no drift across the boundary on the remaining sides of the survey area and therefore use Neumann BCs. The boundary conditions are set using predicted particle concentration data from DREAM.

Parameter specification in the SPDE model
The advection and diffusion parameters in the SPDE (1) must be specified. The velocity field is set from physical assumptions, where c t (s) = (c t e (s), c t n (s)) T as the forcing of the convection in our case is treated as the currents in the ocean. The values are set using predicted currents from the ocean model SINMOD. The isotropic diffusion D is set to be the same as for the computations in the DREAM model D = 0.1 m 2 /s. The initial distribution is mainly estimated from the DREAM data. Since ocean processes are affected by currents, wind patterns, bathymetry, and freshwater runoff, we will have elevated variability in some locations. Thus, we choose a non-stationary prior model (Jun & Stein, 2008) where the mean and variance in each location are specified empirically from the training data available by DREAM: Here, μ i = 1 M M m=1 y i,m is the empirical mean concentration in location s i over all M realizations of training data.
The prior covariance matrix Σ 0 is then given by where Σ i j = σ i σ j R i j and R i j represents the spatial correlation. Hence, the diagonal of the covariance matrix in equation (10) contains the location-wise variances σ 2 i , and the off-diagonal elements describe the covariance between the locations. The covariance function is chosen by comparing the empirical covariance of the ocean model training data with known valid parametric covariance functions. The Matern (3/2) kernel (Matérn, 2013) is fitted to such covariance function here; where h i j = |s i − s j | is the Euclidean distance between two locations s i and s j and φ > 0 is a constant meta-parameter regulating the correlation decay with the distance. The best value for φ is specified using training data. To specify covariance model parameters and visualize how the spatial correlation decays as a function of the distance between points in the domain, it is common in geostatistics to use the variogram of the data, see e.g. Cressie and Wikle (2011). The variogram works on difference X t (s j ) − X t (s i ) between spatial data. Doing so, it tends to be more robust to background trends than the covariance function. The theoretical relation between the variogram γ () and the covariance is (12) Figure 3 shows results from the variogram analysis. The empirical variogram (blue) is computed by averaging the square differences of variable pairs having spatial distance within bins. We use regular bins of distances (first axis). Three different theoretical variograms were fitted to the training data. All of them fit the empirical variogram adequately, but the Matern variogram (green) was chosen here because it is closer to the empirical variogram for very small spatial distances. This indicates that a Matern spatial field has a smoothness that fits the DREAM data better. We see that the variogram increases with distance and then flattens out. Here, the effective spatial correlation length is about 500 m. For the Matern correlation function in equation (11) this means that The advection-diffusion model parameters are further set from the physical assumptions in the ocean model. From visual inspection, this GP model appears to give reasonable simulations of the field but with more patchy randomness than what is seen in the ocean model. As commonly seen in physical models data, the variogram computed from the ocean models data starts to fluctuate after a while. In many cases this can be related to physical principles. Notably, the GP model is a statistical proxy model that can be run onboard but is potentially missing some large-scale oceanographic features. The ocean model is a numerical description with much skill but potentially getting systematic bias due to misspecified inputs growing with time. They both have their merits.
The covariance matrix W t in the noise term of the equation is set to have the same correlation structure as the covariance matrix for the initial state Σ 0 , giving W t = V Σ 0 , where V is a constant scaling parameter. Considering this in relation to the original continuous-time SPDE, the discretization in time would also impact this variance, but in our work, we use a fixed time discretization identical to the one used onboard. There does not seem to be any damping over time in the ocean model data, but for stability reasons (Richardson, 2017) we specify a minimal level λ = −0.0010.

State space model and data assimilation
The state of concentrations is assumed to be modeled by equation (8). The aim of data assimilation is to combine observations with this prior process knowledge to obtain an optimal sequential posterior or filtering estimate of the current state (Sar, 2013). In this study, the prior knowledge is represented by the initial onboard model, and observations are taken from sensors on the AUV.
We define our observation model as where Y t is a length m vector of observations at time t, G t is a time-varying observation matrix of size m × N that maps the state variables to the observations, and ω t represent the observational error which is modeled as a zero-mean Gaussian with standard deviation τ , and I representing the identity matrix. Each row of G t has a one-entry on the vector position associated with an observed grid cell, and the other entries are zero. This method constraints the sampling locations to be at the predefined grid cells. Continuous path sampling, as in Choi and How (2010), could improve the sampling method and should be considered in future work. Since the SPDE model defines a GP with a linear process representation, the spatial Kalman filter can be used to update our state at every time step. Letting Z t = (Y 1 , . . . , Y t ) represent all observations and prior information up to time t, the filtering distribution is X t |Z t ∼ N (μ t|t , Σ t|t ) with mean μ t|t and variance-covariance Σ t|t defined by This recursion starts by the initial distribution with mean μ 0 and Σ 0 , where X 0 ∼ N (μ 0 , Σ 0 ). The boundary conditions and the model parameters going into the matrices and vectors in equation (14) are all specified as described in Sect. 3.2.

Sampling strategy
An objective function is suggested to obtain an informative path for the AUV. The function is created based on three criteria, as in Berget et al. (2019): 1. Locations with high variance are preferred 2. Locations with high predicted concentration are preferred 3. Locations leading to a suitable travel length within the next time interval for the AUV are preferred.
The first criterion is chosen because sampling in areas with a high model variance will reduce the total variance, hence creating a more accurate model. This criterion also ensures that the AUV travels to unexplored areas. The second criterion makes the method adaptive. When studying the simulation results of the particle transport from the complex model, it is clear that the variability is highest where there is a high concentration of particles. Hence, the criterion is inspired by this observation and assumes that locations with high predicted concentration will be rich with information. The last criterion comes from the travel length limitations of the AUV. When choosing the next sample location, it is also essential that the AUV does not travel too far, assuring that the risk of the vehicle drifting outside the survey area is kept low.
The suggested objective function is then created by having a term for the first two criteria. At time step t for location s i , the objective function is given by where the constant parameters θ = [θ 1 , θ 2 ] define the weighting for each criterion. Deciding the weighting of criteria 1 and 2 depends on design choices. When the goal is to obtain a good map of the overall region, then criterion 1 should be dominating. However, if hotspots are the main focus, criterion 2 should dominate. In our case, the parameters θ are tuned in simulation to obtain a good balance between the two criteria. The chosen sampling location s t at time step t is determined as the location that maximizes the objective function f t (s i ) for s i ∈ [s 1 , . . . , s N ]. Given the previous sampling location s t−1 , we choose the next sampling location as The constraints on the objective function ensure that the third criterion in the objective function is adhered to by choosing d min and d max such that a suitable travel length is obtained. |s t−1 − s i | is the Euclidean distance between the previous sampling location and s i . Equation (16) is solved using bruteforce enumeration; the value of the objective function (15) is found for all the grid cells within the region defined by the constraints in (16). The next sampling waypoint s t is chosen as the grid cell with the highest objective value. In this paper, ocean currents are included in the spatio-temporal proxy model, but they are not considered in the AUV navigation plans for selecting sampling locations. Rather, the AUV is set to surface often to avoid risk and reduce navigation uncertainty, as described in the experimental section. In a larger-scale mission, it would certainly be beneficial to include more of the ocean currents in the sampling strat-egy. Making the AUV fly with the currents may increase endurance and reduce travel time between waypoints. This topic is amongst others explored in Pereira et al. (2013) and should be considered in future work.

Summary of the suggested approach
We now provide an overview of the suggested approach. At the beginning of the deployment, the GP model is initialized.
The following steps are executed: -Kernel hyperparameters are estimated using prediction and hindcast data from the particle transportation model DREAM. (Described in Sect. 3.2.) -The initial GP state defined by μ 0 and Σ 0 are found.
(Described in Sect. 3.2.) -Ocean current prediction data from SINMOD is used to define the advection fields c t , which is used together with the diffusion constant D to define transition matrices A t . The sampling method starts out either with the initial GP model or with the latest updated GP model. The AUV performs adaptive sampling choosing the next waypoint as described in Sect. 5. The onboard model of the AUV is then updated; the field is moved with the currents, and new information from observations is added. Updating is executed whenever a waypoint is reached, using data from observations taken along the AUV transect in the assimilation. An overview of the adaptive sampling method is given in Algorithm 1.

Simulation study
A simulation study compares the proposed model with four simpler strategies. We vary both the process model and the sampling strategy. The temporal model based on the advection-diffusion SPDE is compared to a stationary pro- cess model where the next time step is equal to the previous time step, as in Berget et al. (2018). Different sampling strategies are tested with sampling strategies based on the suggested objective function, a predefined sampling strategy following a lawnmower pattern, and a method that uses only the SPDE model with no sampling (no observations are taken). The SPDE and Stationary process model are the same for all sampling strategies using the same parameter values given in the text. For the strategy using the objective function we choose the travel distance restrictions d min = 180 meters and d max = 220 meters and set the time between each sample to be 5 min. This gives a suitable travel speed for the AUV, and a constant time step simplifies the updating of the temporal model. When using the lawnmower pattern, the distance between each sampling location is set to 200 ms, and the same time step of 5 min is used. The parameters in the objective function are set to θ 1 = 15 and θ 2 = 23. The different strategies are shown in Table 1. Ocean model data from two consecutive days, 1st and 2nd of April 2013, is used in the simulation study. The data from April 1st is used to train the SPDE parameters as described in sections 3 and 4. Data from April 2nd is used as the test data, meaning that the predicted field on April 2nd is our ground truth and the field from which observations (samples) are drawn. We use realizations of the GP prior as the initial belief and run 100 replicate runs with the presented sampling methods in Table 1. We set the measurement error standard deviation τ = 0.5 log(ppb). We simulate for a total of 6 h, aiming to map the field on April 2nd from 15:00 to 21:00.

Results and discussion
In this section, we present results from the simulation study and discuss the performance of the methods.

Performance measures
To evaluate the different methods, three performance measures are studied. The mean absolute error (MAE) is computed as the mean of the absolute error between the predicted concentration μ t|t and the true field from DREAM (April 2nd) for all the grid cells and time steps. The mean standard deviation (MSD) is computed similarly but then focuses on the standard deviation in the proxy model in all grid cells at all time steps. Lastly, we compute the mean objective value (see Eq. (15)) (MO) in all grid cells at every time step. Table  2 shows these performance measures computed as the mean over the total time steps. Figure 4 shows how the mean objective value (MO) (see Eq. (15)) evolves through time, plotting the mean value for every time step.
Considering the MAE values in Table 2, we see that SPDEObj performs best in this case. Including a nonstationary process model increases the method's performance, as this is the case using both the Obj sampling strategy and the Lawn sampling strategy. The objective function also increases the method's performance compared to the preplanned lawnmower strategy. Considering the MSD and MO measures, the SPDE method stands out with significantly higher values than the other methods. This is due to the method doing no sampling at all, and hence there is no way for this method to decrease the proxy models' variance, affecting the standard deviation and the objective values.
The effect of increasing variance with no sampling can also be seen in Fig. 4. Again, the SPDE method stands out with steadily increasing values while the other methods sta-bilize at a mean objective value between 200 and 300. For both lawnmower methods (SPDELawn and Lawn), periodical fluctuations in the objective function can be observed. This occurs when the sampling is done close to the east end of the survey area, where the assimilation with new data affects a minor part of the modeled field compared to adding information in the center. Gener-ally, the Obj and SPDEObj methods have a lower overall mean objective function compared with the other methods, indicating an intelligent sampling design.
In Fig. 4, we see that the Lawn method has a lower MO value at some time steps compared to the SPDELawn method. This occurs at times when both Lawn approaches get high MO. Using the lawnmower pattern for sampling (see Fig. 5a) is the reason for this unexpected behavior. Sampling is performed in boxes, and some locations in the middle of the boxes get little influence from measurements. The SPDELawn method solves this by having the temporal model move neighboring sample values into this box. The Lawn method does not have a solution for this, and the values inside this box will be similar to the initial values. Since the initial mine tailings concentration in our example (see Fig. 7a) has low values all over the area, the tailing concentration inside this box will stay low throughout the whole simulation for the Lawn method. This will affect the MO value, leading to artificially low values, even though the performance (here, better illustrated by the MAE) can be poor. Figure 5 shows the AUV path and waypoints for the different sampling strategies plotted on top of the standard deviation at every hour of the simulation study. The same lawnmower pattern is used for SPDELawn and Lawn and can be seen in the first column of the figure. The path chosen by the Obj and the SPDEObj method can be seen in respectively the second and third columns. All the methods cover the area quite well, but both Obj and SPDEObj focus on the eastern area. This region is closer to the outlet, and the true mine tailings distribution (seen in the first column of Fig. 7) shows elevated concentrations in this area. Hence, more frequent sampling here is expected because of the proposed objective function (see Eq. (15)). Figure 6 gives insight into how the GP methods behave, showing the distribution in one single grid cell over time using SPDEObj. The blue line shows the method's mean value, the green line is the true value, and the two red dotted lines are the 95% confidence interval. The point we are considering is 340 ms east and 200 ms north in the survey area. When sampling and assimilation are done close to this point, we see a reduction in the confidence interval, and the mean value tends to be corrected towards the actual value.

Sampling design and field predictions
In Fig. 7 we focus on the two methods Obj and SPDEObj. The updated model mean for every hour is plotted in columns 2 and 3, and the true field (the field where the observations are drawn from) can be seen in the first column. Comparing the methods with the true field, they capture the trend quite well. The display illustrates an important difference between the two methods; because of its stationary process model, the Obj method inherits the structure from the GP through- out the whole mission, relying on observations to change its predicted field. This is not the case with SPDEObj, which has a non-stationary process model that moves the field with the currents. Although data from two different days are used as training and test data, both datasets are from the same ocean models. This may cause the comparison with ground-truth to be optimistic. Still, based on the simulation results, it is reasonable to assume that SPDEObj can map the field quite well. The MAE scores of the two methods are indicated for different times (Fig. 7). With the relatively high dynamical uncertainty, the MAE increases with time. The MAE is generally lower for the SPDEObj approach.
Simulation studies with different AUV velocities have shown that the SPDEObj method's advantages are higher when the AUV takes fewer samples in the survey area over the same time period. Naturally, lower AUV velocities led to fewer samples (comparable to having a bigger survey area), while higher velocities led to an increase in the number of samples. When doubling the number of samples, the difference in the performance between the approaches was minimal, making it hard to decide which method was best. However, when operating with half the number of samples (compared to the simulation study above), the differences became more evident, and the SPDEObj method showed improved results. Hence, it is reasonable to think that the difference between the methods will increase with a larger survey area.

Adaptive versus pre-planned designs
We now give an example comparing an the adaptive approach (SPDEObj) with a pre-planned method (SPDELawn) (see Fig. 8). The example shows the predicted results at the end of a mission run of total 5 h. The true field is seen in Fig. 8a. Figure 8b, c show the model mean and variance for the adaptive approach, and Fig. 8d, e present the model state for the pre-planned approach. The AUV path describing the sam-pling design for one hour before the end time is shown in the variance plots.
In the example, the adaptive approach gives a better estimate of the true field than the pre-planned method. From the true field, we observe increased concentrations on the eastern side. The adaptive approach captures this information, and because of elevated values in the area, continues to investigate the eastern side. The pre-planned lawnmower approach, on the other hand, is set to sample in the western area at the Fig. 8 Example run comparing a pre-planned approach with an adaptive approach. The true field after five hours of sampling can be seen in Fig. 8a. Figure 8b, c show the model mean and variance for the adaptive approach, and Fig. 8d, e present the model state for the pre-planned approach. The AUV path showing the sampling design for one hour before the current time is shown in the variance plots time, and the samples taken are not as informative. Adaptive approaches excel when conditions are not as expected. In such cases, the flexibility of an adaptive approach is beneficial since the sampling design is automatically adjusted. A pre-planned approach will always stick to the pre-scripted sampling design.

Experimental results
As a proof of concept, field experiments were conducted in Fraenfjorden. Two experiments were carried out on December 17th and 18th, 2020. The operational discharge outlet at the time of the experiments is marked as 2.1 in Fig. 1, and a rectangular operational area of size 1150 m × 250 m (also marked in the figure) was used. We focus on the depth layer at around 22-25 ms since this is a layer with a lot of variation and a layer where we, according to the DREAM and SIN-MOD models, will find high concentrations of mine tailings. This is verified by our sensor data as well. Figure 11 shows one transect (waypoint-to-waypoint) for the AUV including sensor data values. It is clear from the figure that a layer with increased turbidity is reached.

Robotic platform and implementation
The robotic platform in our experiments consisted of a Light AUV from OceanScan (see Fig. 9) equipped with a Wetlabs EcoPuck sensor measuring total suspended matter (TSM), which is a measure of optical back-scattering in the form of an attenuation coefficient. The measured sensor values of TSM indicate turbidity, and this observation value was used to assimilate our model. It was converted to a measure of (log) concentration to adjust the value to our onboard model. This conversion is not straightforward, and we need to determine the uncertainty τ . The Wetlabs EcoPuck sen- Fig. 9 The light AUV (Harald) produced by OceanScan was used in our experiments. The length of the AUV is 240 cm, and it is equipped with multiple sensors like CTD, Dissolved Oxygen, and Fluorometer Fig. 10 The architecture of the embedded system (on top) and the simulator (bottom). The Unified Navigation Environment (DUNE) runs onboard the AUV, and the agent architecture T-REX sits on top of this, allowing the embedding of multiple decision processes such as planning sor user manual suggests an attenuation error of about 4 % at attenuation coefficient 1 m −1 . With concentrations varying between approximately 3 and 8 log(ppb), assuming measurement standard deviation of 4% gives τ between 0.12 and 0.32 log(ppb). However, both navigation errors and vehicle drift can lead to uncertainty in the sampling location. Hence, we set τ = 0.5 log(ppb) for the experiments. Figure 10 shows the layout of the agent architecture used in the experimental survey. The Unified Navigation Environment (DUNE) (Pinto et al., 2012) runs onboard the AUV and is used for control, navigation, vehicle supervision, communication, and interaction with actuators. On top of this sits the agent architecture T-REX (Teleo-Reactive EXecutive) (Rajan & Py, 2012), which enables an adaptive mission. T-REX allows the embedding of multiple complex decision processes (including planning). The communication between DUNE and T-REX was handled by the LSTS toolchain (Pinto et al., 2013), which provides the back-seat driver API to DUNE. This allows external controllers, such as T-REX, to provide desired tasks for our platform while receiving progress updates on the current state. Our sampling algorithm was written as a python script and was implemented as a reactor in the T-REX framework. A reactor is a component of T-REX acting as an internal control loop and is capable of producing goals that the planner integrates into a series of actions (e.g., Goto and Arrive_at). The set of all those actions forms a plan built while ensuring operational constraints of the mission (e.g., the vehicle should never dive deeper than a certain depth or leave a defined area.) The plan is then sent to DUNE, which handles low-level control and execution. Before deploying the AUV, our method was tested in a simulated environment similar to our embedded system. In simulations, sensor readings from the AUV were replaced with data from SINMOD and DREAM, and an AUV simulator in DUNE was used to simulate the behavior of the AUV. The layout of the simulator is shown in Fig. 10.

AUV path
Inside the operational area, a 2-dimensional grid of size 58 × 12 was defined, giving N = 684. The distance between each node was approximately 20 m in both grid directions, and each grid point was a possible waypoint for the AUV. The AUV was set to surface whenever reaching a new waypoint to reduce the uncertainty caused by navigation errors and vehicle drift and for safety reasons. The AUV did a yoyo maneuver between 20 and 26 ms between each waypoint. As in the simulation study, we set d min = 180 meters and d max = 220 meters. The AUV speed is set to 1.2m/s, and as opposed to the simulation study, where a constant time of 5 min between two waypoints was used, we now let the time vary. The time used between two waypoints is found using the start and end times of the AUV transect. The temporal model and the updating are done when a waypoint is reached, updating the model with the amount of time used for that particular transect. Figure 11 shows one transect (waypoint-to-waypoint) from our survey, starting with a dive until reaching the selected depth layer, doing the yoyo maneuver, and then surfacing in a corkscrew maneuver. The corkscrew maneuver was chosen to obtain a sampling value as close as possible to the waypoint. The yoyo maneuver was chosen to increase the possibility of finding the layer of elevated mine tailings concentration. The sensor reading used for updating the model was chosen as the maximum value in the layer 22-25 ms (the layer marked with red surfaces) for each of the separate yoyos. For the transect in Fig. 11, this would result in 6 updating values. Because the transect length could vary for different waypoints, the number of yoyos was not fixed but mostly alternated between 5 or 6 yoyos per transect, giving 5 or 6 updating values per transect.

Results
The survey area is located near the coast in a fjord arm, so the currents are highly influenced by the tidal cycle. To get a better insight into how the ocean currents behave at the time of the missions, we consider the tidal cycle at the time in Fig. 12. Our temporal model considers the tides' varia- tions, and we choose to focus on two missions occurring at different intervals in the tidal cycle. The first mission is before a tidal high, and the second is when the tides turn at the top. The period of the two missions is marked in red in the figure. Hourly predicted ocean currents between 09:00 and 13:00 from the SINMOD model can be seen in Fig. 13. As expected, the currents are highly influenced by the tides showing an inward flow (eastward) before the tidal high (09:00-11:00) and weaker ocean currents during the tidal high (11:00-13:00). These predictions are used as input to our onboard model as the forcing in the SPDE model. 8.3.1 Mission 1: 08:53-10:44 The first mission was conducted between 08:53 and 10:44 on December 17th, before a tidal high. The estimated current from the SINMOD model can be seen in Figs. 13a, 13b, 13c. Figure 14 presents the onboard model results. The left column shows the updated mean of the onboard model, representing the mine tailings log(concentration) for selected time steps. The right column shows the onboard model standard deviation at selected time steps and the AUV path up until the given time. The red line represents the AUV path, and the red dots are the waypoints selected by our sampling strategy. The results show increased tailing concentration in the eastern area, especially at the beginning of the mission (see Fig. 14a,14c and 14e). We observe a decay in the concentration in this area after 1.5 h (see Fig. 14g). This is Fig. 13 Hourly ocean current predictions from SINMOD between 09:00 and 13.00. The currents are highly influenced by the tides showing an inward flow (eastward) before the tidal high (09:00-11:00) and weaker ocean currents during the tidal high (11:00-13:00) when the ocean currents get weaker due to the tidal cycle. In the beginning, the chosen sampling strategy seems to focus on exploring the area, covering large parts of the whole rectangle. After the exploration, the focus switches to the eastern area, where high mine tailing concentrations have been observed.

Mission 2: 11:52-13:24
The second mission was conducted during a tidal high. Figure 15 shows results from this survey. The predicted SIN-MOD ocean currents are much weaker than for mission 1 (see Fig. 13d and 13e). Figure 15 presents the onboard model results. The left column shows the updated mean of the onboard model, and the right column shows the standard deviation together with the selected AUV path. This mission does not contain areas of clearly increased tailings concentration, and we see that the AUV chooses a different sampling strategy. Since no area shows increased concentration, there is no extra focus on any special area. Instead, the waypoints are chosen to cover as much as possible of the whole survey area during the whole mission.

Discussion
The chosen sampling paths (see second column of Fig. 14 and  15) follow two different strategies. In the beginning, they are quite similar, focusing on exploring the area, but after the exploration, we see differences in waypoint selection. The sampling path in mission 1 (Fig. 14) focuses on the eastern area where increased concentration is observed, while the sampling path in mission 2, where there are no areas with increased tailings concentration, continues to cover the whole region also after the initial exploration. For both surveys, the standard deviation is reduced at the end of the mission, indicating that the method manages to explore the area. The selection of different strategies shows that the method is adaptive and will change depending on prior data and observations and that the sampling is focused on areas with high tailings concentration. The AUV data are undoubtedly useful for refining the DREAM results.
It is hard to rate the field experiments' performance since there is no ground truth for this case. Still, the model mean concentration (first column of Fig. 14 and 15 coincide well with our assumptions on how the tailings concentration behave. Since the outlet is located on the east side of our area, we assume elevated concentration in the east. It is also reasonable to assume that the concentration of mine tailings in the survey area will be higher when the currents flow westward, as is the case for mission 1. We have also compared the results with the ocean model DREAM data and observed the same tendencies. However, this is not an accurate comparison since the proxy model is trained with data from this model.

Closing remarks
In this work, we propose a method for adaptive sampling of ocean processes with an AUV. Prior information from ocean models is combined with information from in situ measurements to obtain the best possible sampling strategy. A spatio-temporal GP proxy model with a non-stationary covariance function and a process model based on finite difference solutions of an advection-diffusion SPDE are built. This is a low-complexity model which can run onboard the AUV. We can choose informative sampling locations by using an objective function favoring locations with high uncertainty and high predicted particle concentration. The main contribution of the paper is to have a dynamic (spatio-temporal) model onboard a dynamic agent. This enables reliable uncertainty quantification in spatial predictions while actively searching for 'hot-spots' in space and time. The method is tested in simulations and compared with simpler models and sampling strategies. A simulation study shows that the SPDE process model gives improved prediction results in a dynamic environment, and the sampling strategy using the suggested objective function performs better than the sampling strategy using a simpler lawnmower design. Since there is no easy way to obtain a ground truth for field experiments, it is difficult to determine how well the method performs in the field. Still, our results agree well with our assumptions of the mine tailings distribution at the time, and the selection of different sampling strategies for two separate missions illustrates the method's adaptivity. Future work includes expanding the model to 3D (with depth coordinates), adding current estimation from measurements onboard the AUV, including ocean currents in the waypoint selection method, and exploration using multiple AUVs. It is also interesting to test such methods for extended time periods, potentially combining the AUV measurements with satellite data or data from drones.
Funding Open access funding provided by NTNU Norwegian University of Science and Technology (incl St. Olavs Hospital -Trondheim University Hospital). This work is supported by the norges forskningsråd (Grand Nos. 223254, 267793) Data availability The datasets generated during and/or analysed during the current study are available from the corresponding author on reasonable request.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.