Prediction of Changes in Seaﬂoor Depths Based on Time Series of Bathymetry Observations: Dutch North Sea Case

: Guaranteeing safety of navigation within the Netherlands Continental Shelf (NCS), while efﬁciently using its ocean mapping resources, is a key task of Netherlands Hydrographic Service (NLHS) and Rijkswaterstaat (RWS). Resurvey frequencies depend on seaﬂoor dynamics and the aim of this research is to model the seaﬂoor dynamics to predict changes in seaﬂoor depth that would require resurveying. Characterisation of the seaﬂoor dynamics is based on available time series of bathymetry data obtained from the acoustic remote sensing method of both single-beam echosounding (SBES) and multibeam echosounding (MBES). This time series is used to deﬁne a library of mathematical models describing the seaﬂoor dynamics in relation to spatial and temporal changes in depth. An adaptive, functional model selection procedure is developed using a nodal analysis (0D) approach, based on statistical hypothesis testing using a combination of the Overall Model Test (OMT) statistic and Generalised Likelihood Ratio Test (GLRT). This approach ensures that each model has an equal chance of being selected, when more than one hypothesis is plausible for areas that exhibit varying seaﬂoor dynamics. This ensures a more ﬂexible and rigorous decision on the choice of the nominal model assumption. The addition of piecewise linear models to the library offers another characterisation of the trends in the nodal time series. This has led to an optimised model selection procedure and parameterisation of each nodal time series, which is used for the spatial and temporal predictions of the changes in the depths and associated uncertainties. The model selection results show that the models can detect the changes in the seaﬂoor depths with spatial consistency and similarity, particularly in the shoaling areas where tidal sandwaves are present. The predicted changes in depths and uncertainties are translated into a probability risk-alert map by evaluating the probabilities of an indicator variable exceeding a certain decision threshold. This research can further support the decision-making process when optimising resurvey frequencies.


Introduction
The bathymetric surveying necessary for accurate nautical charts and adequate maintenance dredging is a joint responsibility of the Netherlands Hydrographic Service (NLHS) and Rijkswaterstaat (RWS). The seabed mapping is conducted using the active, acoustic remote sensing method of single-beam echosounding (SBES) and multibeam echosounding (MBES). This continuous mapping is a national responsibility and adheres to the SOLAS (Safety of Life at Sea) V convention developed by the International Maritime Organisation (IMO) and is also in accordance with the International Hydrographic Organisation (IHO), Special Publication No. 44 (S44) hydrographic survey standards. The Netherlands has developed and refined its own survey policy, designed to complement the set of standards provided by the IHO S44 standards [1,2].
The Dutch economy is heavily reliant on the activity in the Ports of Rotterdam and Amsterdam and safe navigation in the waterways to these ports, particularly for deeper draught vessels, is considered crucial to prevent risk to vessel grounding. The fundamental knowledge of seabed dynamics in relation to the changes in depths is therefore an important factor to consider for safety of navigation, maintenance of channels and hence for optimising the resurvey policy where the overall aim is to provide timely and accurate predictions of the changes in seafloor depths as a resurvey decision-making tool. This can also ensure more cost-effective and efficient use of the mapping resources given the extensive use of the NCS and the nature of the oceanographic and seafloor properties which make it a complex marine environment.
The NCS is characterised as being a shallow sea, average depth 40 m with a sandy seafloor (typical grain sizes of 0.4 mm, though the spatial distribution of the grain sizes differ) and results in regions of the NCS with varying seabed dynamics i.e., sandbanks which are aligned parallel to the current flow and sandwaves and smaller megaripples which are transverse to the flow [3,4]. Previous research [5][6][7][8], emphasised that characteristics of sandwave features are of the greatest concern in the North Sea since sandwaves vary in amplitudes up to tens of meters, wavelengths in the order 100-1000 m and have migration rates of up to 10 m per year. The migration of the sandwaves are caused by residual currents [4] and tidal asymmetry [9].
For optimisation of resurvey frequencies based on the available time series of bathymetry data, this paper considers an area where tidal sandwaves present challenges to safety of navigation. The library of functional models is characterised to estimate the main parameter of interest, changes in seafloor depth. This required defining the functional and stochastic models using the time series of available bathymetry data. The development of the methodology uses statistical hypothesis testing applied to a library of functional models, to arrive at the most likely model which can be used for predicting the changes in the depths with associated confidence intervals. This method of deformation analysis using time series was also previously studied in [7,[10][11][12]. Studies by [7,12] analysed areas in the 0D (zero-dimensional, per grid node or a nodal analysis), 1D (one-dimensional, line analysis in the directions of maximum and minimum variability) and 2D (planar analysis). A 0D, nodal analysis was also conducted by [8]. These studies concluded that the disadvantage of the 1D and 2D analyses, both required an assumption of the sandwave dynamics in the area of study and hence limited to a fixed area definition.
For this paper, a 0D spatial analysis was chosen to determine the most likely functional model to avoid any spatial assumptions on the seafloor dynamics when defining the library of functional models. An advantage of the choice of 0D analysis is that the patterns related to temporal trends in the changes in depths can be revealed for different nodes (e.g., the crests and the troughs) whereby a distinction can be drawn since they may exhibit different behaviours. In deformation studies such as in [10,13] a 0D approach was also taken and in the case of ground deformation in [14] using Interferometric Synthetic Aperture Radar (InSAR) point scatterers, adding flexibility to the applicability of the 0D approach.
The choice of parameterisation of the spatial and temporal dynamics is achieved by defining a minimal set of model parameters to describe and predict the changes in the depths and associated uncertainties. This primary problem required defining a library of simple, functional models to determine whether changes in depths can be detected by testing mutually exclusive alternatives. Due to the morphological characteristics and migration of sandwaves and hence the time of increase or decrease in the depth within a time series, the parametric models aim to estimate the changes in shallowest likely depth (for the case of safety of navigation which is depicted on nautical charts).
The method is developed to automatically select from the library the most likely model for a reliable prediction. To find a functional model that accurately describes the characteristics of the changing seafloor, it was necessary to take into account the variability of the observations with an a priori stochastic model based on the measurements and gridding process. This functional model selection procedure has been developed based on statistical hypothesis testing and can be viewed as an adapted version of the Detection, Identification and Adaptation (DIA) procedure from [15,16].
Additionally, the functional models selected and the prediction of the changes in depth are used to develop an approach for a risk-alert system for different applications such as safety of navigation, exposure of pipelines and cables on the seafloor [4] and the regeneration of sandwaves after human intervention such as maintenance dredging [17]. This is done by translating the predicted uncertainties of the changes in the depths to estimated probabilities and further evaluating when depth change exceeds a decision threshold. The choice of threshold levels serves as the basis for determining risk-alert levels using quality indicators as a tool for decision-making when improving the planning of resurvey frequencies in the NCS.
For this research, the time series was limited to the preprocessed bathymetry data acquired from SBES and MBES since this combination provides the most accurate measurements, increased coverage and high resolution with regards to MBES, for detecting changes in the depths for the aforementioned applications. This research provides added value to use the methodology developed to assimilate and synthesise depth observations and associated uncertainties derived from other remote sensing techniques such as bathymetric LIDAR (light, detection and ranging) [18], bathymetry retrieval from synthetic aperture radar (SAR) [19] and satellite derived bathymetry (SDB) [20][21][22]. For the NCS however, the efficient use of optical remote sensing techniques such as bathymetric LIDAR and SDB are limited due to reduced water clarity and low turbidity requirements for monitoring bathymetric changes for applications with high accuracy requirements as opposed to MBES.
Although the primary research application for this paper is safety of navigation for optimising resurvey frequencies, other related scientific research in the fields of coastal sea level rise, coastal geomorphology, marine conservation, ocean circulation modelling and marine infrastructure development all rely on the most recent bathymetry data as the basis for their research where monitoring of seafloor changes is essential. Figure 1 shows the present categorized subdivision of the NCS and the corresponding survey frequencies ('Opnameplan') which vary from critical areas, once every 2 years to other areas, once every 4 years, 10 years, 15 years and 25 years, as published within the NLHS organisation in 2017. Due to varying circumstances, the NLHS have not achieved some of these survey frequencies over recent years [7], which makes the prediction and hence the probabilistic risk-alert system using the available time series of bathymetry even more relevant to resurvey planning and decision-making.
The main objectives of this research are (a) the definition of additional functional models which can be used to describe the complex behaviour of the changing seafloor depths especially where sandwave features are present (b) designing a more rigorous approach to select the nominal model to be used in the testing procedure and (c) designing the statistical hypothesis testing procedure to allow for reliable prediction and evaluation of probabilities for risk-alert levels.
The contribution of this research is arranged first with the detailed methodology in Section 2. This section outlines the library of functional models and the description of the stochastic model used for the hypothesis testing procedure developed. The prediction and validation step using the selected functional models is also described here. Section 3 presents the results and discussion of the methodology applied for a site-specific case study area that exhibits tidal sandwave features. Section 4 presents the results of the prediction and uncertainties using the selected models. Section 5 introduces the probability estimation that the indicator variable will exceed a decision threshold, leading to a probability risk map. Section 6 summarises the findings and contributions made.

Preparation of Bathymetry Time Series
Remote sensing technology of SBES and MBES systems provided the time series of bathymetry measurements used for this research. The time series are preprocessed archived data retrieved from the digital Bathymetric Archive System (BAS) of the NLHS and RWS.
NLHS is responsible for the mapping of deeper areas greater than 10 m isobath and RWS, the shallower waterways, nearshore coastal zones and the maintenance of approach channels to port areas [8].
The quality of the archived data is assumed to be at least precise as the IHO S44 standards Order 1 at the 95% confidence level and have been preprocessed and cleaned to remove gross errors and to correct any systematic errors introduced [23]. The high resolution MBES datasets provide an observation at a minimum resolution of 3 m × 5 m. For SBES datasets, there is a minimum line spacing of 50 m and can increase up to 125 m and for this reason interpolation of the SBES surveys for gridded digital elevation model (DEM) generation is necessary. The available time series of bathymetry is in the geodetic datum WGS84 in the UTM Zone 31N projection and the depths are reduced to Lowest Astronomical Tide (LAT).
For the preparation of the time series, there are options for the original data to be gridded on a uniform grid (for e.g., 5 m × 5 m) or using a new quadtree approach [24] for creating a multi-resolution DEM, taking into account the heterogeneity of the data in space and time and the variability of the depths. This allows for optimal use of the time series of data available without loss of details in the final, resulting DEM [25]. Previous research [25] suggested the optimal criteria and standard deviation threshold (measure of variability within a grid cell), σ th values for this application. These suggestions are still flexible depending on the end-users' requirements. The quadtree approach also uses a time ordered approach to store the depths for each epoch or survey moment [25]. Figure 2 shows an example of a single epoch decomposition with σ th = 0.3 m. The input data are therefore gridded so that each grid node represents the mean depth, its variance value and the time of survey defined by the year and month. The variance of each gridded depth value is a combination of the measurement uncertainty (provided by the standards [1]) and the spatial uncertainty (which takes into account the number of observations per grid cell). In [7] the size of the grid cells was limited to 50 m × 50 m and in [8] 25 m × 25 m. For SBES data that requires interpolation, the stochastic interpolation method of Universal Kriging [26] is applied which provides the gridded depth values d and the uncertainty, σ d . Previous studies of [8,27] applied a vertical nodal analysis using linear regression for each node, excluding an uncertainty analysis of the estimated vertical trend as opposed to this contribution. The time series analysis requires the following input:

3.
Observed depths values, d for each grid node, p for each epoch.

Specification of Stochastic Model
The quality of the input bathymetry is a source of uncertainty that can affect the quality and reliability of the DEM and hence the prediction. There are several sources of uncertainties in depth and position measurements since the system of measurement has a combination of its own uncertainties e.g., GPS, gyro and motion sensor. Additionally environmental factors such as sound speed profile, tides and draft. In addition, the uncertainty due to the spatial and interpolation errors needs to be taken into account.

Multibeam Echosounding Variances
The MBES datasets are considered full coverage with a minimum spatial resolution of 5 m × 3 m for each depth measurement, essentially requiring no interpolation. This suggests that the propagation of the MBES uncertainties are a combination of the measurement uncertainties which are derived from the IHO specifications Order 1 and the spatial variability from the gridding technique where the variance of the mean as the estimator is considered.
A main contribution to the DEM uncertainties is thus the measurement uncertainties. This is provided by the IHO standards which suggest the maximum allowable vertical and horizontal uncertainties for each depth measurement resulting from the data collection process [1]. The specifications suggest the maximum, vertical uncertainty for different survey orders depending on the survey priority (which is depth dependent for risk-alert management and safety of navigation) is computed using: Here, a is the coefficient representing the uncertainty that does not vary with depth, b is the coefficient representing the uncertainty that is depth dependent and d is the reduced depth [1]. Sources of errors for the depth measurements and their variances can be used as in [7] where it includes the errors introduced by the echo sounder, sound velocity measurements, static draft measurement, dynamic draft, water level or tidal reduction, heave, pitch and roll. For MBES surveys, variances of the depths are assumed to be at least as precise as the IHO S44 standards of Order 1 uncertainty at the 95% level. Survey Order 1 specifies the depth precision for MBES data as in Equation (3) where the coefficients for computing σ 2 noise in Equation (1) are a = 0.5 m and b = 0.013 respectively. Other errors are the spatial variability of the depths within a specific grid cell size. The grid cell depth is calculated as a (weighted) average of observations in that cell, and uncertainties due to different interpolation or gridding techniques when the data are sparse [28,29]. The spatial uncertainty is therefore attributed to the grid cell size and the number of observations per grid cell. Hence, the grid cell uncertainty footprint can be estimated using: where σ spatial is the spatial uncertainty, σ th is the standard deviation threshold and N, the number of observations within a grid cell. It is assumed that the underlying seafloor is flat or has a stationary mean. This is done during the quadtree decomposition phase as in [25] where the criteria for creating a multi-resolution grid are discussed. Three options are considered when storing the mean depth value for each quadtree grid cell. Spatial homogeneity is assumed when the number of observations is greater than 5 and the sample variance of the depth observations within each grid cell is less than the variance threshold value of σ 2 th = 0.01 m 2 . The total vertical uncertainty of each stored, mean depth value is computed using Equation (2). This estimation is also applied to grid cells where there are less than 5 observations and the grid cell size is considered to be small (less than or equal to 10 m). If, however, there are instances where there are not enough observations within a grid cell, the search area is increased by doubling the current grid size and applying Equation (3).
In general, the combination of propagated errors can be computed using the total propagated uncertainty equation: where we need to use Equations (1) and (2). The kriging variance σ krige is added, but in general for MBES no interpolation is required and hence the value will be zero.

Single-Beam Echosounding Variances
For SBES surveys performed between 1991 and 2003, the total propagated uncertainty in the depth measurement is also a function of various stochastic influences. These stochastic influences include the effect of sound speed velocity, heave, positioning and tidal reduction [7,23] resulting in a posteriori estimate of the variances. For SBES measurements the variances are predominantly due to the interpolation errors, since SBES techniques do not provide full coverage mapping of an area. The geostatistical method of Universal Kriging [30] is used which was adopted from [7] and the DIGIPOL algorithm in [31]. This method models the spatial variability of the seafloor using a Gaussian covariance function. The expected similarity of the observations and the interpolated positions are therefore modelled as a function of distance and direction [23]. The universal kriging provides both the gridded depths and their variances. The total propagated variance for SBES surveys is computed using: as in [7].

Estimated Parameters
The method of Best Linear Unbiased Estimation (BLUE) was used for the estimation of seafloor parameters of the areas selected and the associated uncertainties. Previous studies made assumptions to allow for one-dimensional analyses (1D). The 1D analysis approach required a fixed definition of the areas that are considered uniform in sandwave dynamics and that the depths be gridded using a regular gridded structure. This was a requirement for the 1D analysis used in [7,12] to allow for the characterisation and detection of sandwave dynamics per grid line, in the direction of maximum spatial variability.
For the 0D analysis of an area, there are no a priori assumptions about the areas being homogeneous or exhibiting uniform morphodynamics. Because of this, the size of the areas to be analysed is only limited due to computational time (typically 25-30 min per 3 km × 3 km sub-area). Nodal parameter estimation of seafloor dynamics was studied by [10] and point wise estimation of other deformation parameters applied to geophysical processes and infrastructure using InSAR was studied by [14]. This choice of analysis is based on the available time series data and how we choose to parameterise the changes in depths. This parameterisation is a crucial step that is reliant on the changes that are important in relation to the applications. The advantages of this approach in comparison to the 1D and 2D approaches is that the vertical nodal dynamics will be able to capture the temporal changes in the depth.
The application of deformation analysis for predicting the changes in seafloor depths uses mathematical models as the expression of the relationship between the depths at the gridded locations to the unknown parameters of the model that describe the seafloor. The mathematical model is therefore a combination of the functional and stochastic model: Here d is the m-vector of depth observables, x is the n-vector of unknowns, A is the m × n design matrix describing the assumed linear relation between d and x, Q yy is the m × m variance matrix, e the m-vector of random measurement errors. E{·} and D{·} are the expectation and dispersion operators, respectively. Since it is common to assume the observations are normally distributed, the estimators or library of functional models, with d being a linear function of x is also a general linear Gauss-Markov Model for computing the BLUE of d: The principle of BLUE implies that the estimator must ensure minimum variance and must be unbiased. The BLUE [32] of the unknown parametersx, the adjusted observationŝ d and the residualsê are obtained by: The precision of the estimator is given by: from which can be seen that it depends on the precision of the observations Q yy , the design matrix A and therefore also on the redundancy or the degrees of freedom (DOF) (m − n). The estimated residuals in relation to the variance matrix Q yy affect the decision on the choice of the most likely model and the estimated parameters used for the seafloor estimation.

Functional Model Selection
To characterise the spatial and temporal change of the depths, this behaviour is modelled as a function of the available observations per epoch, t, per grid node. To determine the most appropriate characterisation and hence the model to be used to best describe and predict the behaviour, we apply statistical hypothesis testing in linear models [15]. Prior to applying the hypothesis testing procedure developed, a w-test is performed as a check, for any large errors within the observational data sets, i.e., any outlying survey within the time series [16].
This procedure then starts with a general, goodness of fit approach by defining a library of functional models or multiple hypotheses H (i) where i = 1, 2, 3, 4, .., m where m is the number of hypotheses per time series. This library selection was designed during the exploratory data analysis [33] of the time series of the depths per node (0D) to identify patterns. In addition, empirical studies on the formation and analysis of sandwaves such as [27,34,35] were reviewed. Additionally, the purpose of the predictions was defined i.e., for the purposes of risk-alert management since the extreme changes in the depths pose a concern for safety of navigation, stability of offshore structures and exposure of pipelines and cables on the seafloor.
The seafloor characterisation in time, d s f consists of a library of models using four levels of complexity (refer to Figure 3). The number of hypotheses in the library can be expanded with newly defined characterisations depending on the purpose of each case study. In this instance, the characterisations are kept simple to avoid over-fitting and additional parameter estimations which do not add to the understanding of the behaviour of shoaling depths in NCS. Additionally, the current library considers physically realistic characterisations that can detect and identify significant changes in seafloor depths at any given epoch in the time series. To identify the changes in depth before and after dredging moments, i.e., human intervention, the piecewise linear model characterisation would also prove useful.
The simplest characterisation is represented by a horizontal line where d 0 is constant and is the only unknown parameter. The estimated constant depth is also the mean: Another characterisation, adds the additional slope parameter, creating a constant velocity model with two unknown parameters, initial depth, d 0 and slope, v 0 : The quadratic model has three unknown parameters where a is the acceleration and d 0 and v 0 are the estimated, initial depth and slope parameters respectively: For the piecewise linear characterisation, the model consists of two separate line segments, i.e., for each node, depending on the length of the time series, several piecewise linear models are constructed with a breakpoint at each epoch with the restriction that there is a minimum of two observations at the start and end of each line segment creating the piecewise linear model, refer to Figure 3c. The piecewise linear characterisation consists of four unknown parameters, an initial depth and slope for each line segment respectively, d 0 , v 0 , d 1 and v 1 : for t ≥ t end+1 where, i = 4, . . . , 4 + n p (12) Here n p is the number of piecewise linear models, where n p = m − 2, where m is the number of observations. Figure 3 shows the illustrated representation of the characterisations for each model equation. The total number of hypotheses listed in Table 1, is expanded to include the varying number of piecewise models, n p per time series and therefore the total number of models is more than four per time series.
Similar to previous research studies such as [7,14], this approach does not a priori assume that one hypothesis is more likely than the other, as is the common approach in classical hypothesis testing [15].
To describe the overall model inconsistency separately, for each of the functional models, the general expression of the Overall Model Test (OMT) statistic is first applied and defined as T q(i,j) in Equation (13), where i represents the model number and j the node. The OMT statistic gives the quality of the model approximation by observing the relationship between the residuals,ê, which is the difference between the observations and the adjusted observations using the estimated parameters of the model and the uncertainty of gridded depth values σ 2 d defined by the diagonal of the variance-covariance matrix Q yy .
The OMT statistic is computed separately for all competing models and is compared to a critical value, k α(i,j) . The OMT statistic follows a central chi-squared distribution with (q = m − n) DOF as in Equation (14). This test statistic is computed by defining a critical value k α(i,j) , where α is the probability of falsely rejecting the OMT. This is also referred to as the probability of a false alarm or a Type I error, which should be low to unnecessarily reject a valid hypothesis and therefore the choice of α, or level of significance, is usually chosen to be a small number [16] and can also vary depending on the goal of the analysis [7]. For this contribution the choice of α is 0.001 for step 1(a) as in Figure 4.
Since the library of models vary in complexity and hence has different dimensions, these models will follow different χ 2 distributions. This implies that when comparing the goodness of fit of these models, Equation (13) cannot be used. Instead, to account for the varying complexity of the models (more unknown parameters could result in overfitting), the test statistic in Equation (13) can be normalised by dividing by the critical value, denoted by the statistical significance, α. The test ratio value for each node is defined in Equation (15).
Based on Equation (15) there can be three outcomes for the nodal analysis: 1.
Only one model M (s,j) results in R (i,j) ≤ 1. Then, that model is selected for that particular node.

2.
More than one model can result in R (i,j) ≤ 1. Among the competing models, for each node j, find s = arg min Then M (s,j) is selected for that node. 3.
If no model is accepted based on Equation (16), then none of the models is selected for that node.
The limitations of selecting the most likely model per node, M (s,j) , using these criteria only, are as follows: 1.
The values of two or more of the test ratios can be quite close to each other. Hence the associated models fit almost equally well.

2.
Due to irregularities in the data it can occur that the piecewise linear model is accepted, although in reality the constant velocity model should be selected. Additionally, the length of the time series and the larger number of unknown parameters in the piecewise linear model can result in over-fitting.

3.
The spatial correlation between the nodes is not accounted for and hence the selected models in the neighbouring nodes can show inconsistency, especially in areas of high variability. To avoid misinterpretation, an area-wide assessment is done to make a decision on the choice of H 0 . This is done by selecting the most frequently, acceptable model, according to the OMT statistic, where the choice is between M i=1 and M i=2 as the default models for the initial assumption of general, area-wide dynamics. Between M i=1 and M i=2 , the model resulting in the maximum number of grid nodes is set as the H 0 model. This approach effectively makes the choice of H 0 more adaptable according to the sub-area of interest. This is different from current practice and previous research studies such as [7,14,36,37] where the OMT statistic is applied to a steady-state assumption of H 0 .
Where there are nodes for which the alternative model, H a is preferred instead of the null assumption, H 0 , this could imply that there are errors (disturbances or anomalies). This is another way of testing for large errors in the observations. We assume the m × 1 vector of observations d is normally distributed with known variance matrix Q yy : The following two hypotheses for two competing models are and For the alternative hypothesis, number of additional unknown parameters ∇ (∇ = 0), with vector length q is introduced and is related to the expectation of d by the m × q matrix C d [15]. The General Likelihood Ratio Test (GLRT) provides a decision on whether or not the additional variables in ∇ as in Equation (19) should be taken into account [15].
The GLRT in step 1(b) as in Equation (20) is then applied for each grid node separately using the two competing models where the null assumption H 0 and the alternative model, H a , is set as the assigned model, M (s,j) per grid node after step 1(a). H 0 is rejected in favour of H a if the GLRT value is greater than the critical value, k α(i,j) : The choice of α for step 1(b) is 0.05 and is moderately larger than the choice of α for step 1(a) which was chosen as 0.001. This tradeoff is needed to minimize the effects of a missed detection and to thus maintain the probability of a correct decision after applying the GLRT in step1(b). This implies that if an alternative model, H a is preferred according to the level of significance 0.05 and the degrees of freedom (DOF), then the detection of a deviation from H 0 will be correctly detected and hence H 0 will be rejected when H 0 is false. This is done by keeping the choice of α for step1(b) small to ensure that γ is high. The power of test γ is defined as: The least squares residuals for the null hypothesis, H 0 and the alternative hypotheses, H a are defined as:ê

Prediction and Validation Using Selected Functional Model
Once the final selected model, M (i,j) has been selected as the most likely model using the procedure outlined in Figure 4, this model can be used to describe the evolution of the depths in time, t. The present resurvey frequency map defines the next expected resurvey time, though in practice, it may not always be possible to resurvey according to the planned resurvey frequency. In [7], to predict the depths, the next epoch (in years) is denoted as t S+1 as: where the average period between surveys, T was used as the predicted time and t S as the last available epoch. For this research, the average period T used is two years for the results presented in Section 3. T is a flexible input parameter to ensure the use of short-and long-term prediction times, without the use of predetermined resurvey frequency times.
The predicted depths are represented as an M × 1 vectord p(S+1) as: A p (t) is the known M × N matrix andx is the estimated unknown parameter vector representative of the selected model, M (i,j) . The precision of the predictions is obtained from the variance-covariance matrix Qpp by converting the main diagonal into a column of variances, σ 2 , once the precision of the estimated parameters Qxx is also defined, as: An example of this principle is given in Figure 5 where the 95% and 99% confidence interval around the estimated trend and the predicted trend have been constructed. To provide insight about the reliability and hence the predictive performance of the selected functional models, cross-validation is necessary. The leave-one-out cross-validation (LOOCV) approach [38,39] has limitations due to the short length of the time series and the modelled temporal trends between observations in the time series. To preserve this time ordered dependency during the functional model selection and to ensure a validation procedure, an approach is implemented whereby the time series is extended and hence DEMs and the associated uncertainties are created using simulated observations at the predicted times. To achieve this, the final selected model, M (i,j) was used to solve the forward model using the estimated parameters of M (i,j) . The error distribution, z , for the simulated observations, z (S+1) , is defined as: where σ 2 is defined as the variance of the last observation in the time series. This is derived as the square root of the last element of the variance matrix Q yy . We assume that the simulated observations z (S+1) at different prediction times t S+1 (refer to Equation (23)) can be obtained by simulating the random error z which is added to the predicted observations as in Equation (24), based on the forward model and is given by: With the simulated DEMs, several options can now be applied to assess the crossvalidation error, CV e for an area, i.e., the mean absolute errors for the total number of nodes per area J. First, for the prediction time t S+1 , the average period T used is two years, therefore, the absolute errors between the predictionsd p(S+1) , two years after the last survey moment and simulated observations z (S+1) are computed to obtain the CV e as: Figure 6 gives an example of the simulated observations, z (S+1) at the last survey moment t S and different predicted times, t S+1 where T (years) is [1 2 3 4 5 10 15 30] using the forward model or estimated trend (BLUE) and predicted trend, including the 95% confidence interval. Additionally, since the time series can be extended using the simulated observations and the prediction time, t S+1 is a flexible input parameter, another option can be applied whereby the simulated DEMs are used to further extend the length of the existing time series of the model selection procedure described in Figure 4. This model selection procedure can then be re-executed to evaluate the CV e for an area and prediction time of interest.

Results of Spatial and Temporal Testing
In this section, the spatial and temporal results of the model selection procedure developed are presented for a sub-area within the IJgeul Approach area, West of IJmuiden. In the previous section the methodology for selecting the most likely functional models per grid node and the specification of the stochastic model were outlined. In Section 2, the functional models are different with respect to the specification of each design matrix, A and the stochastic models are the same throughout the testing procedure. The initial assumption of what is known about the general dynamics, that is the null model, H 0 , is selected based on the a priori assumption that each model is equally likely and follow the procedure step 1(a) in Figure 4 of Section 2. The resulting model after step 1(a) which is the most frequently accepted model, gives the corresponding null assumption of the general known dynamics of the area. This result proves to be a good assumption since the results of the model selection show spatial similarities after the GLRT between H 0 vs. M (s,j) is applied as in step 1(b).

Area West of IJmuiden
To test the method using real observations, the area IJgeul shown in Figure 7, located West of IJmuiden, is selected as the case study which was also analysed by [7] since it is an area exhibiting sandwave features. This area is subdivided into 30 sub-areas each of approximately 3 km × 3 km or less. Three sub-areas are the approach anchorage areas, namely AA, AB and AC. A list of the surveys available according to acquisition dates for the whole area is given in Table 2. For positioning, intentional degradation of the GPS signals are referred to as Selective Availability (SA). The older surveys (prior to 1996) used the terrestrial Hyperfix system, which provides the same order of uncertainty as the differential GPS (DGPS) with SA, i.e., in the order of 10 m. All surveys after 1996 used DGPS with SA [7]. According to the resurvey planning, the IJgeul area is resurveyed on average once every 2 years.  Results for Sub-Area NF Figure 8a shows the number of epochs per grid node, which is the length of the time series available per grid node for the sub-area NF. The majority of area NF has a time series length per grid node of 10 available epochs which spans 24 years from 1989 to 2013. The edges of this sub-area show a decrease in the length of the time series because of the size of the buffer areas used to subdivide the whole area into smaller sub-areas in a preprocessing step. However, 97% of the sub-area has a length of time series of 10 which indicates a minimum DOF of 6. Figure 8b shows the computed test ratio, R (i,j) for the selected model M (s,j) in step 1(a). The spatial distribution of R (i,j) shows larger values on the crests of the sandwave features, where there is higher variability in the depths and hence larger residuals after model fitting. Although there are spatial consistencies different models are accepted in step 1(a), mostly M 5 to M 11 . Close to 0% of the total number of nodes are selected for the stable model, M 1 and the constant velocity model M 2 during step 1(a). Figure 8c shows how often a certain model was accepted after step 1(a). However, to avoid misinterpretation due to the limitations outlined in Section 2.4, it is also evaluated which models would be accepted based on the OMT for each node ('OMT statistic only'). From this, it follows that M 2 is the most realistic model to describe the general dynamics of the sub-area, and is therefore selected as the null hypothesis in step 1(b).
After setting the H 0 assumption to M 2 , the GLRT is performed between H 0 vs. M (s,j) . Figure 8d shows the models that are accepted after Equation (20) is applied where H 0 is accepted or rejected in favour of one of the alternative M (s,j) models. The resulting models accepted in Figure 8d show that 93% of the nodes accepted M 2 and the remaining 7% accepted one of the alternative piecewise models, M 5 to M 11 where 4% of those remaining nodes accepted M 11 . Model numbers, M i=5,6,7...11 are characterised as the piecewise linear function and the acceptance of these models is consistent along the crest of the sandwaves where there is higher variability in the depths due to the shape of the sandwave features.
The piecewise linear models are able to identify a significant change in depth, ∆d, where a stable, M 1 , constant velocity, M 2 or quadratic model, M 3 is not adequate to identify the changes between two constrained time intervals, t end and t end+1 . This is further validated through visual inspection of the time series of the observations at selected nodes. Figure 9 shows a zoomed-in section of area NF and the time series of the nodes selected, numbered 1 to 6 in Figure 10.  By visually inspecting the plots in Figure 10, it shows that the accepted models are consistent with the observations in the time series which implies that the change was not misinterpreted as a flat or constantly sloping seafloor. Hence, the testing procedure correctly identifies the ∆d at a particular epoch where the piecewise linear models are accepted in favour of the constant velocity model, M 2 . Figure 11a shows the spatial consistency of the epoch at which the H a models are accepted and the corresponding offsets, ∆d in Figure 11b which are of the same order and spatially consistent as well. Figure 12a,b shows the estimated parameters, ∆d andv of the final selected model.
The results of nodal trends exemplified in Figure 10 also seem to correspond well with a migrating sandwave passing those locations. This is also suggested in Figure 13 where an example of cross-sections through a sandwave feature is extracted in the direction of maximum variability θ x = 59 • in sub-area NF for the different, available epochs. This was done using a 5 m × 5 m gridded DEM for simplicity. The nodal trends of the observations, extracted for different x-locations (local coordinate system) show similar trends as for the models selected in Figure 10. This indicates that the testing procedure developed resulted in selected models that can detect a correct identification of the spatial and temporal change in depths. This then implies that the selected models can be used for predicting the changes in depths, which is described further in Section 4.

Prediction Using the Models Selected: Area NF
Using the most likely models or estimators that were selected using the method outlined in Section 2.4, a prediction with associated uncertainties can be obtained for area NF for t S+1 , two years from the last survey moment. The predictive performance of the selected models were evaluated and the results for the CV e for the sub-area NF as defined in Equation (28) was 0.23 m. Figure 14 shows the relationship between the predicted variables, depth,d p and the predicted slopev at predicted time t S+1 where the variability ofv is greater for the nodes that accept an H a model, and whered p values are on average somewhat smaller than that of the nodes for which the H 0 model was accepted. Using the models that are accepted, we present the spatial results for the prediction at t S+1 , for a projectionT of 2 years after the last epoch, t S . The difference in the depths between the last epoch, t S and t S+1 is computed asĉ: c =d re f −d p (29) and presented as a difference map in meters in Figure 15a. To account for the direction and magnitude of the change in depth, a negative change indicates that the trend in the depths is becoming shallower and a positive change indicates the depths are becoming deeper. The variance ofĉ can be determined using the variance-covariance matrix Qĉĉ as: The main diagonal of Qpp in Equation (25)  Additionally, the σĉ does not provide information on the probability thatĉ will be smaller or greater than a certain value. Here the concept of confidence intervals or regions can be used [32]. Given parameterĉ, and α as the user-defined level of significance, the 100(1 − α)% confidence interval is defined as: where α gives the distance from which the estimatedĉ will be from the true value with a confidence level of 1 − α [32]. Such a confidence interval is often denoted asĉ ± α where commonly used intervals are the 1.96σĉ and 2.58σĉ i.e., the 95% and 99% confidence intervals, respectively. This concept is further explained and applied in the following Section 5 whereĉ is used at the indicator variable to evaluate the probabilities of exceeding a particular threshold value.

Risk-Alert Indicator and Decision Thresholds
In this section, we introduce the problem of estimating the probability that an indicator variable,ĉ (in this case per node, j for a sub-area) will cross an upper and lower decision threshold. This choice of threshold level serves as the basis for determining risk-alert levels using quality indicators as a tool for decision-making when determining the next time of resurvey.
Under the assumption that the estimators selected for the prediction follow a Gaussian distribution, the confidence interval for a given choice of α, denoted as a user-specified probability, can be computed. This provides information on the predicted probability that c is smaller or larger than an upper and lower, threshold value, : Using Equation (32), the 100(1 − α)% confidence interval is defined. The α or significance level is often set to 5% so that the associated confidence interval is 95%. This implies that the predictedĉ will be within a distance of α with a confidence level of (1 − α), as in Equation (32). The choice of α used is 0.05 for the 95% confidence interval. The upper threshold is defined as − and the lower threshold as to distinguish between a shoalingĉ and a deeperĉ respectively.
Advanced geostatistical techniques allow the uncertainty of the estimation to be expressed in terms of risk where [40][41][42] used Disjunctive Kriging and other research used an Indicator Kriging approach based on [43] for evaluating the probability, given the observations, that the indicator variable exceeds a critical threshold. This research contribution of developing a risk-alert indicator uses a simple count of the nodes of the indicator variableĉ where it crosses a decision threshold using the same sub-area of interest, NF. The 95% confidence represent the interval that includes the trueĉ with a probability of 0.95. Varying the value of the for the corresponding α changes the decision threshold at which the probability ofĉ is crossing the upper and lower threshold values.
Values for α are provided by the NLHS as a guide for negative changes in depth (i.e., shoaling depths and the respective decision thresholds for risk-alerts regarding safety of navigation) for different depth ranges as in Table 3. Table 3. Decision thresholds provided as a guide by the NLHS with specific focus on safety of navigation.

Water Depth [m] (LAT)ĉ [m]
An example of estimated probabilities is shown in Figure 16 at the 95% confidence interval for upper and lower threshold values set to [−0.1 m 0.2 m] and [−0.5 m 0.2 m] respectively, where the lower threshold value is kept constant here for comparison and absence of user guidelines. The black and red nodes tend to occur in clusters, suggesting the spatial correlation with the areas that are shoaling and becoming deeper as a result of the shape of sandwave features and its migration. To evaluate other applications such as pipelines and cables on the seafloor, stricter threshold levels may be needed. As the first steps towards determining the upper and lower decision thresholds per sub-area, a range of threshold values for is incorporated and evaluated per node for sub-area NF, shown in Figure 17. The estimation of the probabilities, the choice of the indicator variable used and the adaptable decision thresholds, need to be further evaluated. This will be done to further evaluate the use of the uncertainty of the predictions and hence the approach towards delimiting categorical risk-levels for survey areas of different seafloor dynamics.

Conclusions
An adaptable approach using time series of bathymetry data and a combination of the OMT statistic and the GLRT to obtain the most likely functional model for evaluating a risk-based approach towards optimisation of resurvey frequencies has been developed. Combining the OMT statistic and the GLRT improves the selection of the estimators and accuracy of detection of change in depths.
This library of functional models is built on empirical observations and consists of physically realistic characterisations of the seafloor dynamics and applied to a case study area where sandwaves are dominant, using nodal analysis (0D). The addition of piecewise linear models offers another characterisation of seafloor changes to detect and identify the trends in the nodal time series corresponding to the changes. This hypothesis testing approach ensures that the assumption for the nominal model, H 0 is empirically derived, resulting in an innovative and adaptable approach on the choice of H 0 for the hypothesis testing design and model selection procedure and reliability.
The selection of the most likely functional model allows for short term prediction in changes in seafloor depth with associated uncertainties since the results show spatial consistency where alternative models are selected. This methodology allows for the predictive values to be applied in estimation of the probabilities that the changes in seafloor depths cross an upper and lower decision threshold for the purposes of risk-alert and decision-making for determining the next resurvey moment or dredging intervention. This functional model selection procedure is also relevant to other marine engineering applications where detecting seafloor changes is critical to decision-making, such as buried cables and pipelines on the seafloor. Future research will evaluate further the choice of additional variable indicator(s) and decision thresholds, for risk-alert assessments towards a probability map for decision-making regarding resurveys. The approach developed is applicable to areas that exhibit different seafloor dynamics and when different types of remote sensing data are used or combined. The research contribution in general, is not limited to the use of seafloor depth data.
Funding: This project was part of the larger multidisciplinary SMARTSEA project (number 13275) entitled 'Safe Navigation by optimising seabed monitoring and waterway maintenance using fundamental knowledge for seabed dynamics' which was partly financed by the Netherlands Organisation for Scientific Research (NWO). The data used for this research was provided by the Netherlands Hydrographic Service, Rijkswaterstaat and Deltares.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data used for this research was provided by the Netherlands Hydrographic Service, Rijkswaterstaat and archived by Deltares. Data and further information regarding access to data source can be found at https://publicwiki.deltares.nl/display/OET/Dataset+ documentation+bathymetry+NLHO (accessed in June 2018).

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations and symbols are used in this manuscript: