A multi-collocation method for coastal zone observations with applications to Sentinel-3A altimeter wave height data

In many coastal areas there is an increasing number and variety of observation data available, which are often very heterogeneous in their temporal and spatial sampling characteristics. With the advent of new systems, like the radar altimeter on board the Sentinel-3A satellite, a lot of questions arise concerning the accuracy and added value of different instruments and numerical models. Quantification of errors is a key factor for applications, like data assimilation and forecast improvement. In the past, the triple collocation method to estimate systematic and stochastic errors of measurements and numerical models was successfully applied to different data sets. This method relies on the assumption that three independent data sets provide estimates of the same quantity. In coastal areas with strong gradients even small distances between measurements can lead to larger differences and this assumption can become critical. In this study the triple collocation method is extended in different ways with the specific problems of the coast in mind. In addition to nearest-neighbour approximations considered so far, the presented method allows for use of a large variety of interpolation approaches to take spatial variations in the observed area into account. Observation and numerical model errors can therefore be estimated, even if the distance between the different data sources is too large to assume that they measure the same quantity. If the number of observations is sufficient, the method can also be used to estimate error correlations between certain data source components. As a second novelty, an estimator for the uncertainty in the derived observation errors is derived as a function of the covariance matrices of the input data and the number of available samples.

In the first step, the method is assessed using synthetic observations and Monte Carlo simulations. The technique is then applied to a data set of Sentinel-3A altimeter measurements, in situ wave observations, and numerical wave model data with a focus on the North Sea. Stochastic observation errors for the significant wave height, as well as bias and calibration errors, are derived for the model and the altimeter. The analysis indicates a slight overestimation of altimeter wave heights, which become more pronounced at higher sea states. The smallest stochastic errors are found for the in situ measurements.

Different observation geometries of in situ data and altimeter tracks are furthermore analysed, considering 1-D and 2-D interpolation approaches. For example, the geometry of an altimeter track passing between two in situ wave instruments is considered with model data being available at the in situ locations. It is shown that for a sufficiently large sample, the errors of all data sources, as well as the error correlations of the model, can be estimated with the new method.

1 Introduction

Coastal areas like the German Bight are often characterised by strongly heterogeneous ocean dynamics, typically associated with complicated bathymetry, small-scale coastline features, and river runoff. A few instruments, like high-frequency (HF) radar, are able to capture at least 2-D surface currents with large coverage and high resolution quite nicely. Such systems have a typical range of about 100 km, spatial resolutions on the kilometre scale, and about 20 min sampling (Stanev et al., 2015). However, most instruments provide only point measurements (e.g. buoys) or transects (e.g. satellite altimeter). The combination, and interpretation, of such data is therefore often a challenge. In heterogeneous coastal areas with strong gradients, spatially distributed instruments can observe very different components of the dynamics, even if they are in close proximity.

In the following, this situation is studied in more detail with respect to ocean waves and the significant wave height in particular. Wave height information is of paramount importance for many applications, e.g. shipping, offshore operations, or coastal protection. Although numerical wave forecast models have reached an impressive level of accuracy, there is still room for improvement, in particular in coastal areas with complicated dissipation processes associated with wave breaking and bed friction (Woolf et al., 2002; Reistad et al., 2011; Voorrips et al., 1997; Herbers et al., 2000; Bouws and Komen, 1983; Young et al., 2013; Semedo et al., 2015), as well as with coupling processes between ocean waves, ocean circulation, and the atmosphere (Cavaleri et al., 2018; Staneva et al., 2017; Alari et al., 2016). The focus in this study is on the North Sea, which has interesting ocean wave dynamics mainly caused by the semi-enclosed geometry (Semedo et al., 2015; Voorrips et al., 1997; Boukhanovsky et al., 2007; Staneva et al., 2014). The bathymetry of the considered area with the locations of some in situ wave measurement stations used in the following analysis is shown in Fig. 1.

Figure 1 Bathymetry of the North Sea with the locations of some in situ wave observation instruments considered in this study. The plot shows isobaths for 30, 60, 90, and 120 m water depth.

Traditionally, validations of new data sets are performed by comparing to data from established standard in situ measurements, which are regarded as a reference. As a first step this is acceptable; however, one has to take into account that these reference instruments are affected by measurement errors as well, and the separation of the error contributions from the new data set and the reference instrument is, in general, not possible unless additional information is used. This is easy to see if two data sets, x and y, with uncorrelated additive noise are considered, i.e.

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }x=t+{\mathit{ϵ}}_x,\text{(2)}& {\displaystyle }& {\displaystyle }y=t+{\mathit{ϵ}}_y,\end{array}$$

where t represents the “truth”. If statistics are performed on the difference ξ of x and y, one gets for the mean square error

$$\begin{array}{}\text{(3)}& \mathit{\xi }=\langle (x-y{)}^{\mathrm{2}}\rangle =\langle {\mathit{ϵ}}_x^{\mathrm{2}}\rangle +\langle {\mathit{ϵ}}_y^{\mathrm{2}}\rangle ,\end{array}$$

and it is apparent that it is not possible to derive either $\langle {\mathit{ϵ}}_x^{\mathrm{2}}\rangle$ or $\langle {\mathit{ϵ}}_y^{\mathrm{2}}\rangle$ from ξ alone. The usual approach is therefore to use additional data sets and to make certain a priori assumptions about the errors. If only one data set is added, this leads to the triple collocation method, which has been used and discussed in a number of previous studies (Janssen et al., 2007; Vogelzang and Stoffelen, 2012; Stoffelen, 1998; Caires and Sterl, 2003; McColl et al., 2014). Collocation studies, as presented here, often use a mixture of observations and numerical models. The term “data source” will therefore be used in the following to refer to different types of input data.

In this study the triple collocation approach is extended and adjusted with the special requirements of the coast in mind, where one can usually expect stronger gradients and smaller scale variations than in the open ocean. The objective of the study is to deal with the following four specific issues:

• In the triple collocation method, different information sources within a certain distance are assumed to measure the same quantity, which can be unrealistic in regions with strong gradients, like most coastal areas.

• So far, assumptions about correlation errors were made a priori (Vogelzang and Stoffelen, 2012), but they were not obtained as a result of the collocation process.

• So far, no systematic approach was presented to deal with more than three data sources.

• The quantification of uncertainties concerning estimations of systematic and stochastic data source errors was so far only done based on bootstrap approaches (Caires and Sterl, 2003).

The question about the accuracy of error estimates is of particular concern for new instruments, like Sentinel-3A, for which the amount of available data is still relatively limited. It is also clear that collocation distances are of concern mainly for point measurements or transect observations from satellites. The interpolation of numerical model data to given observation locations is usually less critical if the spatial resolution is appropriate.

The work presented here addresses the issues mentioned above and makes the following main contributions:

• A generalisation of the triple collocation method is introduced, where the truth is not necessarily represented by a single number, but by a more general parameterisation of the truth state that is measured by a group of instruments within a certain distance. The analysis presented here concentrates on 1-D models (i.e. lines), and 2-D models (i.e. planes), but can be easily extended to include more sophisticated approaches.

• In certain configurations, i.e. definitions of truth vectors and spatial distributions of data sources, the approach allows an estimation of cross covariance components of the stochastic errors contained in the considered observations or numerical models.

• The theory includes the definition of a general data source vector, which can contain an arbitrary number of observations and numerical model data.

• Analytical expressions are derived for the estimation errors regarding both systematic calibration errors and stochastic errors of the different data sources.

Like the standard triple collocation method, the extended approach also provides estimates of systematic bias and calibration errors. We will refer to the standard triple collocation method as “TRIPCOL”, and to the multi-collocation as “MULTCOL” in the following.

As an example for the generalised parameterisation of the truth, one can imagine two wave buoys and a satellite altimeter track passing between them. Let us furthermore think about a situation where the wave buoys are too far away from the track to assume that all three instruments measure the same quantity. However, it may be an acceptable assumption that the wave height measured by the altimeter is a linear combination of the wave heights observed by the two buoys. If independent numerical model wave height estimates are available at the buoy locations, the method presented in the following provides a systematic approach to estimate not only the stochastic errors of all data sets, but also the error correlation of the model at the buoy locations.

The present study is supposed to make a contribution to the exploitation of measurements with larger distances, where additional assumptions about the spatial variation in the truth are required. As an illustration, Fig. 2 shows maps of the North Sea with altimeter tracks and collocated buoys with the colour coding referring to the number of obtained collocated data samples within the period April 2016 to August 2017. The data sets will be introduced in more detail in Sect. 3. The plot Fig. 2a shows the situation if a collocation distance of 10 km is assumed as acceptable, whereas Fig. 2b shows the same with a collocation distance of 20 km. One can see that the number of data sets increases rapidly if larger distances are considered.

Figure 2 (a) Map showing Sentinel-3A altimeter tracks together with wave observation platforms with less than 10 km distance to satellite measurements. The colour coding refers to the number of obtained collocated measurements in the period April 2016 to August 2017. (b) The same as (a) with a collocation distance of 20 km.

With regards to the estimation errors, expressions are derived which provide a quantification depending on the covariance matrices of the data sources, and the number of available data samples. These results can give valuable information on the trustworthiness of estimated observation errors, in particular in situations with a small number of samples.

The paper is structured as follows. The multi-collocation method is introduced in Sect. 2. This includes the explanation of the underlying theory for the treatment of the stochastic and systematic errors in Sects. 2.1 and 2.2, as well Monte Carlo simulations to illustrate and verify the method. In Sect. 3 the analysed significant wave heights from in situ stations, Sentinel-3A altimeter, and numerical model wave height data are introduced. As a special case of the multi-collocation method, the triple collocation technique is applied to the wave height data sets in Sect. 4. This includes a new step in the analysis, in which estimation errors are quantified. Section 5 describes the combination of more than three observations taken at a certain distance to estimate measurement errors and error correlations.

2 Multi-collocation method

In this section the multi-collocation method is explained, which includes the triple collocation technique as a special case. In the first step, the approach for the estimation of the stochastic errors is presented, and in the second part systematic bias and calibration errors are considered.

2.1 Symmetric approach

The approach presented in this section to estimate stochastic errors does not require bias-free reference instruments. Calibration errors are not considered in this first step. Let us assume the truth is given by a vector t of dimension n_t, and n_o data sources, $y_{\mathrm{1}},\mathrm{\dots },y_{n_o}$, are related to the truth by

$$\begin{array}{}\text{(4)}& \mathit{y}=\mathbf{A}\mathit{t}+\mathit{ϵ}+\mathit{b}.\end{array}$$

Here, A is an n_o×n_t matrix; ϵ is an n_o-dimensional zero mean Gaussian process, which represents the stochastic data source errors; and b is a vector of length n_o containing the biases of the different data source components. Bold type is used for vectors. The triple collocation method is then a spacial case with n_t=1, n_o=3, and $\mathbf{A}=(\mathrm{1},\mathrm{1},\mathrm{1}{)}^T$. Here and in the following, the symbol T denotes the transpose operation. This case will be considered in Sect. 4 looking at a larger number of in situ observation locations in the North Sea. Using different definitions of the truth vector and the matrix A, various relationships between the truth and the data sources can be formulated with the above approach. In this study, we will concentrate on 1-D and 2-D linear models. It should be emphasised that the truth cannot, in general, be represented by a finite number of parameters. However, it is reasonable to assume that the reality is sufficiently smooth, and hence a Taylor expansion can be applied. The triple collocation method is then a special case, where only the constant term is considered. Depending on the number of available observations, the approach in Eq. (4) allows for the addition of higher-order terms. We will concentrate on linear approximations in this study; however, the method is able to deal with interpolation approaches of higher order if a sufficient number of data sources is available. Conceptually, this issue is related to the topic of representation errors (e.g. Van Leeuwen, 2015). The 1-D case will be considered in the Monte Carlo simulations presented in Sect. 2.5, as well as in Sect. 5.1. The 2-D case will be discussed in Sect. 5.2.

Let us now define a matrix B, which contains a basis of the null-space of A as rows. This can, for example, be obtained by singular value decomposition of A and selecting the eigenvectors corresponding to vanishing eigenvalues. If A has full rank, B is a $(n_o-n_t)\times n_o$ matrix. For the triple collocation method this leads to

$$\begin{array}{}\text{(5)}& \mathbf{B}={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}}}}\left(\begin{array}{ccc}\mathrm{1}& -\mathrm{1}& \mathrm{0}\\ \mathrm{1}& \mathrm{0}& -\mathrm{1}\end{array}\right).\end{array}$$

Multiplying Eq. (4) from the left by B gives

$$\begin{array}{}\text{(6)}& \mathbf{By}=\mathbf{B}\mathit{ϵ}+\mathbf{B}\mathit{b}.\end{array}$$

Averaging over all measurements then leads to

$$\begin{array}{}\text{(7)}& \langle \mathbf{By}\rangle =\mathbf{B}\mathit{b}.\end{array}$$

Forming the second-order moments results in

$$\begin{array}{ll}{\displaystyle }\mathbf{B}\langle {\mathbf{yy}}^T\rangle {\mathbf{B}}^T& {\displaystyle }-\langle \mathbf{By}\rangle \langle \mathbf{By}{\rangle }^T=\mathbf{BA}\langle {\mathbf{tt}}^T\rangle {\mathbf{A}}^T{\mathbf{B}}^T+\mathbf{B}\langle \mathit{ϵ}{\mathit{ϵ}}^T\rangle {\mathbf{B}}^T\\ \text{(8)}& {\displaystyle }& {\displaystyle }=\mathbf{B}\langle \mathit{ϵ}{\mathit{ϵ}}^T\rangle {\mathbf{B}}^T=:\mathbf{Z},\end{array}$$

where we have a symmetric $(n_o-n_t)\times (n_o-n_t)$ matrix on both sides of the equation. Because of the symmetry, one gets

$$\begin{array}{}\text{(9)}& m={\displaystyle \frac{(n_o-n_t{)}^{\mathrm{2}}+(n_o-n_t)}{\mathrm{2}}}\end{array}$$

equations. The right-hand side Z is of the form

$$\begin{array}{}\text{(10)}& {\displaystyle }{Z}_{ij}& {\displaystyle }=\sum _{q,k=\mathrm{1}}^{n_o}\langle {\mathit{ϵ}}_q{\mathit{ϵ}}_k\rangle B_{iq}{B}_{jk}=\sum _{k=\mathrm{1}}^{n_o}\langle |{\mathit{ϵ}}_k{|}^{\mathrm{2}}\rangle B_{ik}{B}_{jk}\text{(11)}& {\displaystyle }& {\displaystyle }+\sum _{q<k}\langle {\mathit{ϵ}}_q{\mathit{ϵ}}_k\rangle (B_{iq}{B}_{jk}+B_{ik}{B}_{jq}).\end{array}$$

Eq. (8) is therefore a linear system of equations of the form

$$\begin{array}{}\text{(12)}& \mathit{r}=\mathbf{D}\mathit{ϵ},\end{array}$$

where the vector ϵ contains the unknown variances and covariances of ϵ and r contains elements of the matrix on the left-hand side of Eq. (8). If it is possible to limit the number of unknowns to m, or less, using appropriate assumptions about the variance structure (e.g. independence of error components), this system can be solved when the corresponding system matrix D is regular. Table 1 summarises some feasible combinations of n_t, n_o, and the number of error variances n_var and covariances n_covar that can be estimated if D is regular. Possible observation system configurations corresponding to these cases are shown in Fig. 3. Here, Fig. 3a corresponds to the standard TRIPCOL approach, where all data sources within a certain distance are assumed to measure the same truth. Linear approximations in 1-D and 2-D used in the MULTCOL approach, which relate data sources with a larger distance, are depicted in Fig. 3b and c, respectively.

Table 1 The number of data source error variances n_var and covariances n_covar that can be estimated using different dimensions of the truth parameterisation n_t, and number of data sources n_o according to Eq. (9).

Figure 3 Illustration of three considered observation scenarios. (a) Data sources are assumed to provide information on the same quantity (0-D). (b) Measurements are located along a line and a 1-D linear approximation is employed for the measured quantity (1-D). (c) The same as (b), but using a 2-D interpolation method for more general spatial distributions of data sources (2-D).

If there are more equations than unknowns, a standard linear squares approach can be used to find a reasonable estimate for the unknown variance and covariance components of ϵ. It is interesting to note that this approach also works for biased measurements, although it is in general not possible to estimate the bias explicitly. All that is required is an estimate of Bb and this is easy to obtain via averaging of Eq. (7).

For the case of the triple collocation method, the system matrix D is in fact regular and the inverse is given by

$$\begin{array}{}\text{(13)}& {\mathbf{D}}^{-\mathrm{1}}=\mathrm{2}\left(\begin{array}{ccc}\mathrm{0}& \mathrm{0}& \mathrm{1}\\ \mathrm{1}& \mathrm{0}& -\mathrm{1}\\ \mathrm{0}& \mathrm{1}& -\mathrm{1}\end{array}\right).\end{array}$$

For the triple collocation problem this leads to the well known expressions for the stochastic error variances (Janssen et al., 2007).

$$\begin{array}{}\text{(14)}& {\displaystyle }& {\displaystyle }\langle {\mathit{ϵ}}_{\mathrm{1}}^{\mathrm{2}}\rangle =\langle (y_{\mathrm{1}}-y_{\mathrm{2}})(y_{\mathrm{1}}-y_{\mathrm{3}})\rangle ,\text{(15)}& {\displaystyle }& {\displaystyle }\langle {\mathit{ϵ}}_{\mathrm{2}}^{\mathrm{2}}\rangle =\langle (y_{\mathrm{2}}-y_{\mathrm{1}})(y_{\mathrm{2}}-y_{\mathrm{3}})\rangle ,\text{(16)}& {\displaystyle }& {\displaystyle }\langle {\mathit{ϵ}}_{\mathrm{3}}^{\mathrm{2}}\rangle =\langle (y_{\mathrm{3}}-y_{\mathrm{2}})(y_{\mathrm{3}}-y_{\mathrm{1}})\rangle .\end{array}$$

This corresponds to the “0-D” case in Table 1 and the geometry in Fig. 3a.

If the available number of samples n_s is small, the estimated observation errors may be affected by large errors. To quantify these uncertainties at least in an approximate way, the covariance of the covariance estimator

$$\begin{array}{}\text{(17)}& \stackrel{\mathrm{‾}}{\mathrm{COVAR}}(x_i,x_j)={\displaystyle \frac{\mathrm{1}}{n_{\mathrm{s}}}}\sum _{q=\mathrm{1}}^{n_{\mathrm{s}}}{x}_i^q{x}_j^q\end{array}$$

is considered, where the stochastic vector $(x_{\mathrm{1}},x_{\mathrm{2}},\mathrm{\dots })$ is assumed to be Gaussian and zero mean. The covariance of these estimators ${\mathit{\chi }}_{i,j,i^{\prime },j^{\prime }}$ for different pairs of (i,j) and $(i^{\prime },j^{\prime })$ can then be written as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\chi }}_{i,j,i^{\prime },j^{\prime }}:=\mathrm{COVAR}\left(\stackrel{\mathrm{‾}}{\mathrm{COVAR}}\right(x_i,x_j),\stackrel{\mathrm{‾}}{\mathrm{COVAR}}(x_{i^{\prime }}{x}_{j^{\prime }}\left)\right)\\ \text{(18)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{n_s^{\mathrm{2}}}}\sum _{q{q}^{\prime }}\langle x_i^q{x}_j^q{x}_{i^{\prime }}^{q^{\prime }}{x}_{j^{\prime }}^{q^{\prime }}\rangle -{\displaystyle \frac{\mathrm{1}}{n_s^{\mathrm{2}}}}\sum _{q{q}^{\prime }}\langle x_i^q{x}_j^q\rangle \langle x_{i^{\prime }}^{q^{\prime }}{x}_{j^{\prime }}^{q^{\prime }}\rangle .\end{array}$$

Using standard relationships for the higher-order central moments of Gaussian-distributed variables (Triantafyllopoulos, 2003), this can be expressed as

$$\begin{array}{ll}{\displaystyle }{\mathit{\chi }}_{i,j,i^{\prime },j^{\prime }}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{n_{\mathrm{s}}}}\mathrm{COVAR}(x_i,x_{i^{\prime }})\mathrm{COVAR}(x_j,x_{j^{\prime }})\\ \text{(19)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{1}}{n_{\mathrm{s}}}}\mathrm{COVAR}(x_i,x_{j^{\prime }})\mathrm{COVAR}(x_j,x_{i^{\prime }}).\end{array}$$

The latter expression for ${\mathit{\chi }}_{i,j,i^{\prime },j^{\prime }}$ can be used to estimate the variances and covariances in the estimation errors on the left-hand side of Eqs. (8) and (12), respectively. Therefore, the uncertainties in the estimated vector $\stackrel{\mathrm{‾}}{\mathit{ϵ}}$ can be approximated by

$$\begin{array}{}\text{(20)}& \mathrm{covar}\left(\stackrel{\mathrm{‾}}{\mathit{ϵ}}\right)={\mathbf{D}}^{-\mathrm{1}}\mathrm{covar}\left(\mathit{r}\right)({\mathbf{D}}^{-\mathrm{1}}{)}^T.\end{array}$$

From Eqs. (19) and (20) it is evident that observations with large variance and strong positive correlations will tend to lead to stronger estimation errors for $\stackrel{\mathrm{‾}}{\mathit{ϵ}}$. This is in particular the case when the geophysical background statistics already contribute a lot of variance, or when measurements are within the correlation distance of the background fields and the uncorrelated observation errors are relatively small. The usefulness of the approximation Eq. (20) will be considered in Sect. 2.3 based on Monte Carlo simulations.

2.2 Use of reference instruments

In this section a special, but also typical, situation is considered, where for a couple of measurements systematic errors can be neglected. Typically, this assumption is made for standard in situ observations systems, like waverider buoys (Janssen et al., 2007) or wind anemometers (Stoffelen, 1998). In this case, the error model for the different data sources can be formulated as follows:

$$\begin{array}{}\text{(21)}& \left(\begin{array}{c}\mathit{x}\\ \mathit{y}\end{array}\right)=\left(\begin{array}{c}\mathbf{I}\\ \mathit{\lambda }\end{array}\right)\left(\begin{array}{c}{\mathbf{A}}_x\\ {\mathbf{A}}_y\end{array}\right)t+\left(\begin{array}{c}{\mathit{ϵ}}_x\\ {\mathit{ϵ}}_y\end{array}\right)+\left(\begin{array}{c}\mathrm{0}\\ {\mathbf{b}}_y\end{array}\right).\end{array}$$

Here, x represents the vector of reference measurements, and y contains the remaining data sources. In the examples discussed in the following sections, x will contain in situ wave height measurements and y will represent a combination of satellite altimeter and numerical wave model data. The dimensions of x and y are denoted by n_x and n_y in the following. The matrices A_x and A_y translate the truth vector t to the expected reference measurements x and the other data sources y. In addition, it is assumed that the matrix A_x is invertible, i.e. it is possible to obtain an estimate of the truth vector t from the observations x. The matrix I is the identity matrix. Apart from a possible bias, the vector of data sources y may be also affected by systematic calibration errors represented by the diagonal matrix λ.

To obtain expressions for the scaling parameters contained in λ, the first- and second-order moments of the input data x and y are considered. For the first-order moments M_x and M_y of x and y, one gets

$$\begin{array}{}\text{(22)}& {\displaystyle }& {\displaystyle }{M}_x={\mathbf{A}}_x\langle \mathbf{t}\rangle ,\text{(23)}& {\displaystyle }& {\displaystyle }{M}_y=\mathit{\lambda }{\mathbf{A}}_y\langle \mathbf{t}\rangle +{\mathbf{b}}_{\mathbf{y}}=\mathit{\lambda }{\mathbf{A}}_y{\mathbf{A}}_x^{-\mathrm{1}}{M}_x+{\mathit{b}}_y.\end{array}$$

The second-order moments M_xx and M_yy follow as

$$\begin{array}{ll}\text{(24)}& {\displaystyle }{M}_{xx}& {\displaystyle }={\mathbf{A}}_x\langle {\mathbf{tt}}^T\rangle {\mathbf{A}}_x^T+\langle {\mathit{ϵ}}_x{\mathit{ϵ}}_x^T\rangle ,{\displaystyle }{M}_{yy}& {\displaystyle }=\mathit{\lambda }{\mathbf{A}}_y\langle {\mathbf{tt}}^T\rangle {\mathbf{A}}_y^T\mathit{\lambda }+\langle {\mathit{ϵ}}_y{\mathit{ϵ}}_y^T\rangle +{\mathbf{b}}_y{\mathbf{b}}_y^T+\mathit{\lambda }{\mathbf{A}}_y\langle \mathbf{t}\rangle {\mathbf{b}}_y^T\\ \text{(25)}& {\displaystyle }& {\displaystyle }+{\mathbf{b}}_y\langle \mathbf{t}\rangle {\mathbf{A}}_y^T\mathit{\lambda }.\end{array}$$

The covariance functions C_xx and C_yy can then be written as

$$\begin{array}{ll}\text{(26)}& {\displaystyle }{C}_{xx}& {\displaystyle }={\mathbf{A}}_x\langle {\mathbf{tt}}^T\rangle {\mathbf{A}}_x^T+\langle {\mathit{ϵ}}_x{\mathit{ϵ}}_x^T\rangle -{\mathbf{A}}_x\langle \mathbf{t}\rangle \langle \mathbf{t}{\rangle }^T{\mathbf{A}}_x^T,{\displaystyle }{C}_{yy}& {\displaystyle }=\mathit{\lambda }{\mathbf{A}}_y\langle {\mathbf{tt}}^T\rangle {\mathbf{A}}_y^T\mathit{\lambda }+\langle {\mathit{ϵ}}_y{\mathit{ϵ}}_y^T\rangle +{\mathbf{b}}_y{\mathbf{b}}_y^T+\mathit{\lambda }{\mathbf{A}}_y\langle \mathbf{t}\rangle {{\mathbf{b}}_{\mathbf{y}}}^T\\ {\displaystyle }& {\displaystyle }+{\mathbf{b}}_{\mathbf{y}}\langle \mathbf{t}{\rangle }^T{\mathbf{A}}_y^T\mathit{\lambda }-\left(\mathit{\lambda }{\mathbf{A}}_y{\mathbf{A}}_x^{-\mathrm{1}}{\mathbf{M}}_x+{\mathbf{b}}_y\right)\\ \text{(27)}& {\displaystyle }& {\displaystyle }{\left(\mathit{\lambda }{\mathbf{A}}_y{\mathbf{A}}_x^{-\mathrm{1}}{\mathbf{M}}_x+{\mathbf{b}}_y\right)}^T{\displaystyle }& {\displaystyle }=\mathit{\lambda }{\mathbf{A}}_y\langle {\mathbf{tt}}^T\rangle {\mathbf{A}}_y^T\mathit{\lambda }-\mathit{\lambda }{\mathbf{A}}_y{\mathbf{A}}_x^{-\mathrm{1}}{\mathbf{M}}_x{\mathbf{M}}_x^T({\mathbf{A}}_x^{-\mathrm{1}}{)}^T{\mathbf{A}}_y^T\mathit{\lambda }\\ \text{(28)}& {\displaystyle }& {\displaystyle }+\langle {\mathit{ϵ}}_y{\mathit{ϵ}}_y^T\rangle \end{array}$$

and the cross covariance C_xy between x and y as

$$\begin{array}{}\text{(29)}& {\displaystyle }& {\displaystyle }{\mathbf{C}}_{xy}={\mathbf{A}}_x\langle {\mathbf{tt}}^T\rangle {\mathbf{A}}_y^T\mathit{\lambda }+\langle {\mathit{ϵ}}_x{\mathit{ϵ}}_y^T\rangle -{\mathbf{A}}_x\langle \mathbf{t}\rangle \langle \mathbf{t}{\rangle }^T{\mathbf{A}}_y^T\mathit{\lambda }\text{(30)}& {\displaystyle }& {\displaystyle }\mathit{\lambda }{\mathbf{A}}_y{\mathbf{A}}_x^{-\mathrm{1}}{C}_{xy}=\mathit{\lambda }{\mathbf{A}}_y\langle {\mathbf{tt}}^T\rangle {\mathbf{A}}_y^T\mathit{\lambda }-\mathit{\lambda }{\mathbf{A}}_y{\mathbf{A}}_x^{-\mathrm{1}}{M}_x{M}_x^T({\mathbf{A}}_x^{-\mathrm{1}}{)}^T{\mathbf{A}}_y^T\mathit{\lambda }.\end{array}$$

The equation for C_yy then gives

$$\begin{array}{}\text{(31)}& {\displaystyle }{\mathbf{C}}_{yy}=\mathit{\lambda }{\mathbf{A}}_y{\mathbf{A}}_x^{-\mathrm{1}}{\mathbf{C}}_{xy}+\langle {\mathit{ϵ}}_y{\mathit{ϵ}}_y^T\rangle .\end{array}$$

This results in n_y equations for each scaling component according to

$$\begin{array}{}\text{(32)}& {\mathit{\lambda }}_i={\displaystyle \frac{C_{ij}-\langle {\mathit{ϵ}}_i{\mathit{ϵ}}_j\rangle }{{\sum }_q{\mathit{\nu }}_{iq}{C}_{qj}}}:={\displaystyle \frac{{\mathrm{\Omega }}_{\mathrm{1}}}{{\mathrm{\Omega }}_{\mathrm{2}}}}j=\mathrm{1},\mathrm{\dots },n_y,\end{array}$$

where ν_iq are the elements of the matrix ${\mathbf{A}}_y{\mathbf{A}}_x^{-\mathrm{1}}$.

Let us assume for a moment that the scaling parameters λ are available. One can then derive the bias of y from Eq. (23). Furthermore, defining the matrix A in Eq. (4) as

$$\begin{array}{}\text{(33)}& \mathbf{A}=\left(\begin{array}{c}{\mathbf{A}}_x\\ \mathit{\lambda }{\mathbf{A}}_y\end{array}\right),\end{array}$$

the approach in Sect. 2.1 can be applied to estimate the stochastic errors of the different data sources.

There are now two basic approaches to estimate the scaling parameters.

• Direct method: Those terms in Eq. (32) are used for which 〈ϵ_iϵ_j〉 is known, e.g. because the error components are assumed independent. In this case the estimation of the observation errors and the scaling parameters are independent and can be treated separately.

• Iterative method: Terms in Eq. (32) are used for which 〈ϵ_iϵ_j〉 is not known a priori. In this case an iterative method has to be used, where the estimation of the data source errors and the scaling parameters are performed in succession until convergence is achieved. Similar iteration techniques were also discussed for the triple collocation method in Janssen et al. (2007) and Vogelzang and Stoffelen (2012).

In Janssen et al. (2007) an iterative method had to be applied for the triple collocation analysis, because the proposed procedure for the scaling parameter estimation leads to a non-linear expression, which could not be treated in a direct way. The direct method for the standard triple collocation problem leads to the known expressions (Caires and Sterl, 2003):

$$\begin{array}{}\text{(34)}& {\displaystyle }& {\displaystyle }{\mathit{\lambda }}_{y_{\mathrm{1}}}={\displaystyle \frac{C_{y_{\mathrm{1}}{y}_{\mathrm{2}}}}{C_{x_{\mathrm{1}}{y}_{\mathrm{2}}}}},\text{(35)}& {\displaystyle }& {\displaystyle }{\mathit{\lambda }}_{y_{\mathrm{2}}}={\displaystyle \frac{C_{y_{\mathrm{1}}{y}_{\mathrm{2}}}}{C_{x_{\mathrm{1}}{y}_{\mathrm{1}}}}}.\end{array}$$

One can see that for the estimation of ${\mathit{\lambda }}_{y_{\mathrm{1}}}$ no use is made of correlations between y₁ and x₁, which may contain a lot of useful information. This can be overcome by the iterative version with

$$\begin{array}{}\text{(36)}& {\displaystyle }& {\displaystyle }{\mathit{\lambda }}_{y_{\mathrm{1}}}={\displaystyle \frac{C_{y_{\mathrm{1}}{y}_{\mathrm{1}}}-\langle |{\mathit{ϵ}}_{y_{\mathrm{1}}}{|}^{\mathrm{2}}\rangle }{C_{x_{\mathrm{1}}{y}_{\mathrm{1}}}}},\text{(37)}& {\displaystyle }& {\displaystyle }{\mathit{\lambda }}_{y_{\mathrm{2}}}={\displaystyle \frac{C_{y_{\mathrm{2}}{y}_{\mathrm{2}}}-\langle |{\mathit{ϵ}}_{y_{\mathrm{2}}}{|}^{\mathrm{2}}\rangle }{C_{x_{\mathrm{1}}{y}_{\mathrm{2}}}}}.\end{array}$$

In some cases, there may be several equations for one component of λ, and it is then important to have an approximation for the respective estimation errors to pick the estimator with the smallest variance. Quantification of these uncertainties is also of general interest in the statistical analysis of data, in particular if the sample size is small. For the analysis in this study, we only consider the direct method, where the 〈ϵ_iϵ_j〉 in Eq. (32) are known constants. We also do not consider the additional uncertainty, which is caused by estimation errors for these stochastic error variances and covariances. Denoting the nominator and denominator in Eq. (32) by Ω₁ and Ω₂, a Taylor expansion gives

$$\begin{array}{}\text{(38)}& {\mathit{\lambda }}_k\approx {\displaystyle \frac{\mathrm{1}}{\langle {\mathrm{\Omega }}_{\mathrm{2}}\rangle }}({\mathrm{\Omega }}_{\mathrm{1}}-\langle {\mathrm{\Omega }}_{\mathrm{1}}\rangle )-{\displaystyle \frac{\langle {\mathrm{\Omega }}_{\mathrm{1}}\rangle }{\langle {\mathrm{\Omega }}_{\mathrm{2}}{\rangle }^{\mathrm{2}}}}({\mathrm{\Omega }}_{\mathrm{2}}-\langle {\mathrm{\Omega }}_{\mathrm{2}}\rangle )+{\displaystyle \frac{\langle {\mathrm{\Omega }}_{\mathrm{1}}\rangle }{\langle {\mathrm{\Omega }}_{\mathrm{2}}\rangle }}.\end{array}$$

For the variance one gets

$$\begin{array}{ll}{\displaystyle }\mathrm{var}\left({\mathit{\lambda }}_k\right)& {\displaystyle }={\displaystyle \frac{\mathrm{var}\left({\mathrm{\Omega }}_{\mathrm{1}}\right)}{\langle {\mathrm{\Omega }}_{\mathrm{2}}{\rangle }^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{var}\left({\mathrm{\Omega }}_{\mathrm{2}}\right)\langle {\mathrm{\Omega }}_{\mathrm{1}}{\rangle }^{\mathrm{2}}}{\langle {\mathrm{\Omega }}_{\mathrm{2}}{\rangle }^{\mathrm{4}}}}\\ \text{(39)}& {\displaystyle }& {\displaystyle }-\mathrm{2}{\displaystyle \frac{\mathrm{covar}({\mathrm{\Omega }}_{\mathrm{1}},{\mathrm{\Omega }}_{\mathrm{2}})\langle {\mathrm{\Omega }}_{\mathrm{1}}\rangle }{\langle {\mathrm{\Omega }}_{\mathrm{2}}{\rangle }^{\mathrm{3}}}}.\end{array}$$

The variances and covariances in Ω₁ and Ω₂ can be derived by making, again, use of Eq. (19).

2.3 Generation of background statistics

In the following, the techniques presented in Sects. 2.1 and 2.2 will be assessed and validated based on synthetic observations for which the observation errors are known a priori. This requires Monte Carlo simulations for which realistic background statistics are desirable. Here, we use parameters derived from an 11-month time series of two buoys in the German Bight. The buoys “ELB” and “HEL” can be identified in Fig. 1 as the instruments closest to the entrance of the river Elbe. The buoy HEL is near the island Heligoland at about 25 m water depth and about 30 km north-west of the buoy ELB, which is at about 27 m water depth. The wave height distributions of both buoys shown in Fig. 4b and c indicate a shape, which can be very well approximated with a log-normal distribution superimposed as green curves. The joint distribution in Fig. 4a shows a quite good correlation between the data sets, which is expected due to the relative close proximity of the buoys. The histogram of the difference between the Elbe buoy and the Heligoland buoy shown in Fig. 4d indicates that the majority of cases have higher waves at the Heligoland location than the Elbe location. This makes sense because north-westerly winds are predominant in the area. Therefore, situations with waves coming from offshore and being dissipated by wave breaking and bottom friction are most often observed in the German Bight. The fewer cases with higher waves near Heligoland are associated with southerly winds, where waves are actually generated near the coast and the wave height increases with fetch length. The respective parameters for the log-normal distribution including the correlations of both buoy time series are given in Table 2.

Figure 4 Background statistics used for the Monte Carlo Simulations. (a) 2-D histogram of the joint distribution for the Elbe and Heligoland buoys derived from the period June 2016–April 2017 with diagonal given as a dashed green line and the black isolines referring to probability density. (b) 1-D histogram for the Heligoland buoy wave heights with log-normal probability density function (pdf) superimposed in green. (c) The same as (b) for the Elbe buoy. (d) Histogram of the difference between the Elbe and Heligoland wave height with Gaussian pdf superimposed in green, and the dashed red line indicating the zero position.

Table 2 Mean, variance (var), covariance (covar), and correlation (corr) parameters used for the simulation of the background wave height statistics at the locations of the Heligoland and Elbe buoys in the German Bight. These numbers were derived from measurements taken during the period June 2016–April 2017. The respective probability distributions with a log-normal approximation are shown in Fig. 4.

2.4 Impact of coastal gradients and spatial data set resolutions on triple collocations

In this section a brief analysis is presented concerning the impact of coastal gradients on the standard triple collocation approach and the role of spatial data set resolutions. The analysis is illustrated using the background statistics presented in Sect. 2.3.

As explained before, the triple collocation method makes the assumption that all three data sets represent the same truth. We consider the case now, where this assumption is violated, and where we have data sets representing the wave height at three different locations x, x^′, and x^′′. Let us denote the wave heights at the three locations by

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{H}_x={\stackrel{\mathrm{‾}}{H}}_x+{\widehat{H}}_x,\\ \text{(40)}& {\displaystyle }& {\displaystyle }{H}_{x^{\prime }}={\stackrel{\mathrm{‾}}{H}}_{x^{\prime }}+{\widehat{H}}_{x^{\prime }},{\displaystyle }& {\displaystyle }{H}_{x^{\prime \prime }}={\stackrel{\mathrm{‾}}{H}}_{x^{\prime \prime }}+{\widehat{H}}_{x^{\prime \prime }},\end{array}$$

where ${\stackrel{\mathrm{‾}}{H}}_x,{\stackrel{\mathrm{‾}}{H}}_{x^{\prime }}$, and ${\stackrel{\mathrm{‾}}{H}}_{x^{\prime \prime }}$ are the respective mean values, and ${\widehat{H}}_x,{\widehat{H}}_{x^{\prime }}$, and ${\widehat{H}}_{x^{\prime \prime }}$ are the departures from the mean. Furthermore, it is assumed that the three wave height data sets are affected by uncorrelated additive zero mean errors ${\mathit{ϵ}}_x,{\mathit{ϵ}}_{x^{\prime }}$, and ${\mathit{ϵ}}_{x^{\prime \prime }}$. According to Eq. (14), the measurement error of the data source corresponding to location x would be estimated as follows if the standard triple collocation method is applied:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\langle {\mathit{ϵ}}_x^{\mathrm{2}}\rangle \approx \langle \left({\widehat{H}}_x+{\stackrel{\mathrm{‾}}{H}}_x+{\mathit{ϵ}}_x-{\widehat{H}}_{x^{\prime }}-{\stackrel{\mathrm{‾}}{H}}_{x^{\prime }}-{\mathit{ϵ}}_{x^{\prime }}\right)\\ {\displaystyle }& {\displaystyle }\left({\widehat{H}}_x+{\stackrel{\mathrm{‾}}{H}}_x+{\mathit{ϵ}}_x-{\widehat{H}}_{x^{\prime \prime }}-{\stackrel{\mathrm{‾}}{H}}_{x^{\prime \prime }}-{\mathit{ϵ}}_{x^{\prime \prime }}\right)\rangle \\ {\displaystyle }& {\displaystyle }=\langle {\mathit{ϵ}}_x^{\mathrm{2}}\rangle +\langle \left({\widehat{H}}_x-{\widehat{H}}_{x^{\prime }}\right)\left({\widehat{H}}_x-{\widehat{H}}_{x^{\prime \prime }}\right)\rangle \\ \text{(41)}& {\displaystyle }& {\displaystyle }+\left({\stackrel{\mathrm{‾}}{H}}_x-{\stackrel{\mathrm{‾}}{H}}_{x^{\prime }}\right)\left({\stackrel{\mathrm{‾}}{H}}_x-{\stackrel{\mathrm{‾}}{H}}_{x^{\prime \prime }}\right)=:\langle {\mathit{ϵ}}_x^{\mathrm{2}}\rangle +R_x+{\stackrel{\mathrm{‾}}{R}}_x.\end{array}$$

The angle brackets refer to averages over different realisations of the background state and data source errors. As one can see, the estimation of $\langle {\mathit{ϵ}}_x^{\mathrm{2}}\rangle$ is affected by an error which has two components. The term R_x is related to correlations of wave height differences in the background statistics. In situations where all three data sources are in a region with a spatial wave height gradient, typically observed in coastal areas, this term will not vanish, at least as long as the data sources are located along the gradient. The term ${\stackrel{\mathrm{‾}}{R}}_x$ is related to the differences in mean wave heights at the three locations. This term can be expected to contribute to the estimation error in coastal areas as well.

We are now estimating these error contributions for the background statistics derived in Sect. 2.3. Let us assume that the wave height along a straight line between the stations HEL and ELB can be approximated reasonably well with a linear function, i.e.

$$\begin{array}{}\text{(42)}& H_x={\displaystyle \frac{x(H_{\mathrm{Hel}}-H_{\mathrm{Elb}})+d{H}_{\mathrm{Elb}}}{d}},\end{array}$$

where x denotes the distance of some point X from the Elbe station in the direction of Heligoland; H_Hel and H_Elb are the wave heights at the Heligoland and Elbe stations; and d=24 km is the distance between the two stations. Defining ${\widehat{H}}_{\mathrm{Hel}}$ and ${\widehat{H}}_{\mathrm{Elb}}$ analogues to Eq. (40), we get for the wave height covariance of two points $x,x^{\prime }$

$$\begin{array}{ll}{\displaystyle }\langle {\widehat{H}}_x{\widehat{H}}_{x^{\prime }}\rangle & {\displaystyle }={\displaystyle \frac{\mathrm{1}}{d^{\mathrm{2}}}}\langle \left(x\right({\widehat{H}}_{\mathrm{Hel}}-{\widehat{H}}_{\mathrm{Elb}})+d{\widehat{H}}_{\mathrm{Elb}})\\ {\displaystyle }& {\displaystyle }\left(x^{\prime }\right({\widehat{H}}_{\mathrm{Hel}}-{\widehat{H}}_{\mathrm{Elb}})+d{\widehat{H}}_{\mathrm{Elb}})\rangle \\ {\displaystyle }& {\displaystyle }=x{x}^{\prime }\langle {\displaystyle \frac{({\widehat{H}}_{\mathrm{Hel}}-{\widehat{H}}_{\mathrm{Elb}}{)}^{\mathrm{2}}}{d^{\mathrm{2}}}}\rangle +(x+x^{\prime })\\ \text{(43)}& {\displaystyle }& {\displaystyle }\langle {\widehat{H}}_{\mathrm{Elb}}{\displaystyle \frac{({\widehat{H}}_{\mathrm{Hel}}-{\widehat{H}}_{\mathrm{Elb}})}{d}}\rangle +\langle {\widehat{H}}_{\mathrm{Elb}}^{\mathrm{2}}\rangle {\displaystyle }& {\displaystyle }=:x{x}^{\prime }{\mathit{\alpha }}_{\mathrm{1}}+(x+x^{\prime }){\mathit{\alpha }}_{\mathrm{2}}+{\mathit{\alpha }}_{\mathrm{3}},\end{array}$$

where α₁=0.000147, α₂=0.003174 m, and α₃=1.7529 m² for the considered case. With this information Eq. (41) can be evaluated. For the spatial distribution of the three data sources we assume that x=10, $x^{\prime }=\mathrm{5}$, and $x^{\prime \prime }=\mathrm{15}$ km, i.e. all three data sources are within 10 km distance. The resulting errors for the triple collocation method then follow as (see Eq. 41):

$$\begin{array}{ll}{\displaystyle }{R}_x^{\mathrm{10}\mathrm{km}}& {\displaystyle }+{\stackrel{\mathrm{‾}}{R}}_x^{\mathrm{10}\mathrm{km}}=-\mathrm{0.003675}{\mathrm{m}}^{\mathrm{2}}-\mathrm{0.00035}{\mathrm{m}}^{\mathrm{2}}\\ \text{(44)}& {\displaystyle }& {\displaystyle }=-\mathrm{0.0040}{\mathrm{m}}^{\mathrm{2}}.\end{array}$$

If an instrument at location x is considered, which has a truth observation error standard deviation of $\sqrt{\langle {\mathit{ϵ}}_x^{\mathrm{2}}\rangle }=\mathrm{0.1}$ m, these error terms would lead to an estimation error by the triple collocation method of about 40 % in terms of variance. If the collocation distance is increased to 20 km, and one has data sources at x=10, $x^{\prime }=\mathrm{0}$, and $x^{\prime \prime }=\mathrm{20}$ km, the following error is obtained:

$$\begin{array}{ll}{\displaystyle }{R}_x^{\mathrm{20}\mathrm{km}}& {\displaystyle }+{\stackrel{\mathrm{‾}}{R}}_x^{\mathrm{20}\mathrm{km}}=-\mathrm{0.01473}{\mathrm{m}}^{\mathrm{2}}-\mathrm{0.00138}{\mathrm{m}}^{\mathrm{2}}\\ \text{(45)}& {\displaystyle }& {\displaystyle }=-\mathrm{0.0161}{\mathrm{m}}^{\mathrm{2}}\end{array}$$

In this case the collocation error grows to 160 % with respect to the truth observation error. As explained in Sect. 2.1, the multi-collocation method proposed in this study is designed to take spatial gradients, as discussed above, into account; however, this is at the cost of requiring a larger number of data sources. This will be illustrated in Sect. 2.5 using the same background statistics.

The second issue to be discussed in this section is the role of the spatial resolution of the different models and observations. The main point to consider here is that sub-resolution variations in wave height become part of the estimated data set error if the triple or multi-collocation methods are applied. This has two main consequences:

• The estimated data source errors are influenced by the background statistics.

• For two data sources with a common unresolved band of spatial scales, the data source errors are correlated (Vogelzang and Stoffelen, 2012).

In general, the expected unresolved sub-resolution variance in wave height is given by

$$\begin{array}{}\text{(46)}& H_{\mathrm{sub}}^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{A}}\underset{A}{\int }\langle (H_{x^{\prime }}-{\stackrel{\mathrm{‾}}{H}}_x^{\mathrm{data}}{)}^{\mathrm{2}}\rangle \mathrm{d}{x}^{\prime },\end{array}$$

where A is the resolution cell of the assumed data source and the respective data source wave height ${\stackrel{\mathrm{‾}}{H}}_x^{\mathrm{data}}$ is computed as

$$\begin{array}{}\text{(47)}& {\stackrel{\mathrm{‾}}{H}}_x^{\mathrm{data}}={\displaystyle \frac{\mathrm{1}}{A}}\underset{A}{\int }{H}_{x^{\prime }}\mathrm{d}{x}^{\prime },\end{array}$$

where $H_{x^{\prime }}$ is the truth wave height at location x^′ within the resolution cell. We now evaluate these terms, again using the background statistics presented in Sect. 2.3. For simplicity, we assume that the resolution cell is one-dimensional and spans from the Elbe station (x=0) to some point x=a in the direction of Heligoland. The data set then corresponds to averages of the form

$$\begin{array}{}\text{(48)}& {\stackrel{\mathrm{‾}}{H}}_{a/\mathrm{2}}^{\mathrm{data}}={\displaystyle \frac{\frac{a}{\mathrm{2}}(H_{\mathrm{Hel}}-H_{\mathrm{Elb}})+\mathrm{d}{H}_{\mathrm{Elb}}}{d}}.\end{array}$$

The mean unresolved variance within one resolution cell can then be written as

$$\begin{array}{ll}{\displaystyle }{H}_{\mathrm{sub}}^{\mathrm{2}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{a}}\underset{\mathrm{0}}{\overset{a}{\int }}\langle ({\displaystyle \frac{x^{\prime }(H_{\mathrm{Hel}}-H_{\mathrm{Elb}})+d{H}_{\mathrm{Elb}}}{d}}-{\stackrel{\mathrm{‾}}{H}}_{a/\mathrm{2}}^{\mathrm{data}}{)}^{\mathrm{2}}\rangle \mathrm{d}{x}^{\prime }\\ \text{(49)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{a}}\underset{\mathrm{0}}{\overset{a}{\int }}\langle ({\displaystyle \frac{(x^{\prime }-a/\mathrm{2})(H_{\mathrm{Hel}}-H_{\mathrm{Elb}})}{d}}{)}^{\mathrm{2}}\rangle \mathrm{d}{x}^{\prime }{\displaystyle }& {\displaystyle }={\displaystyle \frac{a^{\mathrm{2}}}{\mathrm{12}}}({\displaystyle \frac{\langle ({\widehat{H}}_{\mathrm{Hel}}-{\widehat{H}}_{\mathrm{Elb}}{)}^{\mathrm{2}}\rangle }{d^{\mathrm{2}}}}+{\displaystyle \frac{({\stackrel{\mathrm{‾}}{H}}_{\mathrm{Hel}}-{\stackrel{\mathrm{‾}}{H}}_{\mathrm{Elb}}{)}^{\mathrm{2}}}{d^{\mathrm{2}}}}).\end{array}$$

One can see that the mean sub-resolution variance is depending on the mean gradient, as well as the variance in the gradient within the resolution cell. Using the background statistic values in Sect. 2.3, the variance $H_{\mathrm{sub}}^{\mathrm{2}}$ was computed for different values of the cell size a. For a=5 km one gets $H_{\mathrm{sub}}^{\mathrm{2}}=\mathrm{0.0003}$m², for a=10 km the result is $H_{\mathrm{sub}}^{\mathrm{2}}=\mathrm{0.0012}$m², and a=20 km gives $H_{\mathrm{sub}}^{\mathrm{2}}=\mathrm{0.005}$m². Let us imagine an observation instrument located at $x=a/\mathrm{2}$ with a measurement error standard deviation of 0.1 m. This value is supposed to only describe the instrumental errors, i.e. the errors that one always has, even if the wave height within the resolution cell is constant. Furthermore, assume that the observations represent averages over wave heights in the resolution cell of size a as described by Eq. (48). For a=10 km this averaging process adds about 10 % to the data set error variance, and for a=20 km this increases to 50 %.

The above analysis has shown that both the collocation distance and the spatial data set resolutions are important factors for the quantification and interpretation of the respective data set errors. The separation of instrumental errors and sub-resolution-related errors is a challenge, because it requires knowledge about the truth background wave statistics on a sub-resolution scale. In general, such information can only be obtained if one of the data sources has a significantly higher spatial resolution than the other data sources.

2.5 Monte Carlo simulation for the 1-D case

As an example, we consider the case where we have data sources which are approximately located along a straight line. This corresponds to the scenario depicted in Fig. 3b. We approximate the truth state by a linear model with two parameters. From Eq. (9) it follows that we need at least five data sources to estimate the errors. Let us assume we have two buoys, a satellite altimeter with two measurements close to the buoys, and a numerical model estimate in the middle between the two buoys. Using the wave heights at the buoy positions as the state parameters t one gets

$$\begin{array}{}\text{(50)}& \mathbf{A}=\left(\begin{array}{cc}\mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}\\ \mathrm{1}/\mathrm{7}& \mathrm{6}/\mathrm{7}\\ \mathrm{6}/\mathrm{7}& \mathrm{1}/\mathrm{7}\\ \mathrm{1}/\mathrm{2}& \mathrm{1}/\mathrm{2}\end{array}\right)\end{array}$$

for the matrix A, which relates the truth vector to the observations (see Eq. 4). Here, we have assumed a geometry as depicted in Fig. 3b. The first and second row of A refer to the two buoy measurements, which are assumed to be without systematic errors. The third and fourth row correspond to the two altimeter measurements near the ELB buoys and the HEL buoy, which are assumed to be affected by calibration errors with scaling parameters of 1.2 and 1.3. The last row represents the wave height estimate provided by the wave model in the middle between the two buoys. The model is assumed to have a calibration error with a factor of 0.9.

The Monte Carlo experiments were then performed as follows:

• 120 observation vectors y were created using a random simulator with prescribed variances and covariances for the background statistics and the observation errors.

• The observation errors and their uncertainty were estimated using the approach described in Sect. 2.1.

• These experiments were repeated 1000 times to obtain statistically robust results.

The parameters used for the simulations, as well as the obtained results, are summarised in Table 3. The first three columns refer to the assumed observation error statistics for the buoys, the altimeter, and the numerical model. One can see that a covariance of 0.056 m² was used for the two satellite measurements, which corresponds to an error correlation of 0.5. The last three columns refer to the estimation errors, which were obtained in two different ways:

• The uncertainties are estimated directly by computing the variance of the estimated observation errors over all experiments. This is called the “averaged experiments” approach (AVEXP) in the following.

• The uncertainties were estimated for each experiment from the input data covariance matrices as explained in Sect. 2.1. These estimates were than averaged over all experiments. This is called the “covariance matrix” approach (COMAT) in the following.

For the obtained data source errors averaged over all experiments, the numbers agree with the assumed errors to within three decimals, which illustrates the validity of Eq. (12). The same is also true for the estimated uncertainties for the variances and covariances estimated from Eq. (20). The last three columns in Table 3 show that the covariance matrix method and the numbers from the averaged experiments are in very good agreement. The last column contains the respective comparison for the covariance of the altimeter measurement errors, where the two approaches also give very consistent results. Overall, these results confirm that the estimation of uncertainties in the estimated stochastic errors by Eq. (20) is a reasonable approach.

Table 3 Parameters used for the Monte Carlo simulations in Sect. 2.5. The first two columns refer to the stochastic wave height error standard deviation (SD) and variance (var) assumed for the considered data sources. The third column gives the assumed error cross covariance (covar) values for the two altimeter measurements and the two buoy data sets. The fourth column is the error standard deviation of the estimator for the observation error variances obtained by averaging over 1000 estimation experiments (AVEXP approach). The values in column five refer to the same estimation errors, but derived by application of the method described in Sect. 2.1 (COMAT approach). The last column gives the COMAT and AVEXP standard deviations for the covariance estimation errors.

In a second step the same exercise was done for the estimation of the systematic errors. The first column of Table 4 shows the assumed calibration errors, i.e. scaling parameters used in the generation of the synthetic observations. In this case the estimated calibration factors averaged over all experiments shown in the second column agree with the theoretical values to within two decimals, which seems reasonable. The values for the estimation errors obtained with the COMAT approach (third column) and the AVEXP approach (fourth column) are also in good agreement, considering that several approximations (e.g. Eq. 38) were used.

Table 4 Parameters used for the Monte Carlo simulations to validate the approach described in Sect. 2.2 for the quantification of errors in the calibration factor estimates. The MULTCOL technique is applied to a 1-D configuration with five data sources, of which three (two altimeter (Alt) observations and one model estimate) are affected by calibration errors. The first column gives the assumed truth scaling parameters and the second column gives the respective estimates. The last two columns represent the uncertainty estimates for the derived scaling parameters in terms of standard deviation (SD) based on two different procedures. See text for details.

3 Description of data sets

In this section the observation and numerical model data used for the multi-collocation analysis are introduced. The data sets are from the period April 2016 to August 2017.

3.1 Satellite altimeter data

The satellite data used here were taken by the European satellite Sentinel-3A launched in February 2016. The satellite flies on a sun-synchronous orbit with an exact repeat cycle of 27 days. The spatial accuracy of the revisit is ± 1 km in the longitudinal direction. Among other instruments, the platform hosts a radar altimeter (SRAL) operating at Ku and C bands (Le Roy et al., 2007). The main frequency used for range measurements is in the Ku band (13.575 GHz), while the C-band frequency (5.41 GHz) is used for ionospheric correction. The basic footprint of the altimeter antenna is a disc with approximately 20 km diameter. However, the effective area actually influencing the measurements is more narrowly centred around the nadir point with a diameter of about

$$\begin{array}{}\text{(51)}& A={\displaystyle \frac{\mathit{\pi }{R}_{\mathrm{0}}(c\mathit{\tau }+\mathrm{2}{H}_{\mathrm{s}})}{\mathrm{1}+R_{\mathrm{0}}/R_{\mathrm{e}}}}.\end{array}$$

Here, R₀=814 km is the altitude of the satellite, R_e is the radius of the earth, c is the speed of light, τ is the pulse duration, and H_s is significant wave height (Chelton et al., 1989). For the typical pulse durations of the order of 3 ns, the effective footprint varies between 1 and 10 km with larger footprints at high sea states. In particular in coastal areas, the altimeter data processing is quite involved (Chelton et al., 2001), and a number of instrument and processing parameters can have a strong impact on the characteristics of the wave height estimates.

In this study Sentinel-3A data with 1 Hz sampling are analysed, which correspond to measurements taken every 7 km along the track. The analysed data were acquired in the so-called reduced synthetic aperture radar (RDSAR) mode, which provides data comparable to measurements from a traditional satellite altimeter. A comparison of different Sentinel-3A altimeter modes can be found in Wiese et al. (2018).

Figure 2a shows the distribution of Sentinel-3A tracks over the North Sea. “Ascending” passes are from south-south-east to north-north-west, whereas “descending” passes are from north-north-east to south-south-west.

3.2 In situ measurements

In this study in situ wave height measurements distributed over the Global Telecommunication System (GTS) were used, which are archived at the European Centre for Medium-Range Weather Forecasts (ECMWF; Bidlot and Holt, 2006). Additional wave observation data were gathered by ECMWF as part of the JCOMM Forecast Verification project (Bidlot et al., 2002). These measurements have a quite inhomogeneous geographical distribution as shown in Fig. 2. As one can see, the focus of the observations is on coastal areas and regions with intense offshore activities, like the northern part of the North Sea. Some of the in situ stations shown in Figs. 1 and 2b, which are referenced in the subsequent analysis, are labelled by either 5 digit numbers (e.g. 62168) or three character strings (e.g. ELB). Due to the lack of respective metadata, it was not possible to distinguish between different types of instruments, e.g. waverider buoys or platform-mounted devices. One exception is the station 62170 near the east English Channel entrance, which is identical to the light ship “F3” mentioned in Anderson et al. (2016). In addition to the GTS data, in situ wave measurements taken in the German Bight were obtained from the Bundesamt für Seeschifffahrt und Hydrographie (BSH). The GTS data have a temporal sampling of 1 h, while the BSH buoys provide observations every 30 min. The in situ observations represent raw values and were checked for unrealistic wave heights. Looking at all the in situ stations for the analysed period in summary, the provided significant wave heights were in the range between 0.1 and 7.8 m. These are realistic values for the North Sea (Semedo et al., 2015) and hence all observations were used in the analysis.

3.3 Wave model WAM and meteorological input data used

For this study, data generated with the spectral wave model WAM were used (Komen et al., 1996). The model version cycle 4.6.2 considered here includes depth refraction and wave breaking and is therefore suitable for coastal applications (Staneva et al., 2017). Spatial variations in bathymetry are taken into account; however, temporal variations in water depth due to tides are not included in the simulations. The 2-D-wave spectra are calculated on a polar grid with 30 directional 15^∘ sectors and 30 logarithmically spaced frequencies ranging from 0.042 to 0.66 Hz. A spherical grid is used for the space dimensions with ∼0.06^∘ resolution in zonal and ∼0.03^∘ resolution in meridional direction. The required forcing at the open boundaries of the North Sea model domain are derived from a coarser model simulation for the whole North Atlantic. Model output with 1 h time steps was available for the analysis. ERA-5 data are used as meteorological forcing for the North Sea model runs (Hersbach and Dee, 2016). This data set is a global re-analysis product from ECMWF with a spatial resolution of 31 km. The model results are interpolated to a 0.25^∘ grid, and the time step is 1 h in the final product. A detailed comparison of different model setups with satellite altimeter data can be found in Wiese et al. (2018).

Compared to previous studies (Janssen et al., 2007; Caires and Sterl, 2003), the spatial resolutions of the three analysed data sources are in quite close agreement. The effective resolutions of the altimeter and the in situ instruments both depend on the actual sea state. For the altimeter typical footprint sizes are between 1 and 10 km as explained in Sect. 3.1. For the in situ data, the translation of the typical 20 min averages to spatial averages is determined by the group velocity. For example, the energy will propagate with about 15 km h⁻¹ if the dominant wave length is 50 m long and the water is deep (>50 m). A 20 min temporal average would therefore correspond to a 5 km spatial average in this case, which is in good correspondence to the spatial model resolution of about 3.5 km. We have therefore used the original data for the analysis and not generated super-observations by averaging, as done in Janssen et al. (2007) and Caires and Sterl (2003), who used wave model data with significantly coarser resolution.

4 Triple collocations for the entire North Sea

In this section the triple collocation method, as a special case of the multi-collocation approach, is applied to the Sentinel-3A altimeter wave height measurements introduced in Sect. 3.1 to assess the respective systematic and stochastic errors. The novelty lies in the analysis of a new satellite data set and the provision of error bars for the estimated stochastic and systematic errors.

Traditionally, validations of new data sets are performed by comparing to data from established in situ measurements, which are regarded as a reference. Here, the following assumptions are made:

• Sentinel-3A and the WAM model may be affected by calibration problems represented by the calibrations factors λ_S3 and λ_WAM.

• Sentinel-3A and the WAM model may be affected by biases b_WAM and b_S3.

• Buoys are regarded as reference systems, i.e. they are assumed as bias free and without calibration errors.

Each of the Sentinel-3A tracks shown in Fig. 2a is passed by the satellite about once a month. Figure 2b shows the respective number of collocations found if a maximum distance of 10 km is accepted. The collocation involves some necessary interpolation steps, which were performed as follows (Janssen et al., 2007):

• The model is interpolated to the buoy using linear interpolation.

• The model is interpolated to the closest altimeter point using linear interpolation.

• Both the buoy and the model are interpolated to the satellite overflight time.

• The model value used for the location is taken as the average of the buoy and the satellite interpolation (see Janssen et al., 2007).

The triple collocation technique was applied to each in situ platform for which altimeter data within the acceptable collocation could be found. The direct method as described in Sect. 2.2 was used for this analysis.

As an example, Fig. 5 shows the obtained results for the Elbe buoy ELB located at 54.0^∘ N 8.1^∘ E. The location of this buoy can also be found in Fig. 1. The three scatter plots show the used data sets in different combinations (buoy versus WAM, Fig. 5a; buoy versus Sentinel-3A, Fig. 5b; and WAM versus Sentinel-3A, Fig. 5c). The three data sets were corrected according to the slope and bias parameters estimated in the collocation procedure. The slope parameters for both the model and Sentinel-3A were found to be below 1, and there exists a larger positive bias for the altimeter. The red triangles correspond to ascending satellite passes and green triangles indicate descending satellite heading. A connection between the satellite flight direction and errors is not clearly visible. This is an important result, because the altimeter data processing is known to be more challenging for passes going from land to sea. It is evident that the best agreement is between the buoy and the model. The smallest stochastic error is found for the buoy with 0.04 m standard deviation. For this location, the collocation procedure gives the largest stochastic error of 0.25 m for the altimeter data.

Figure 5 (a, b) Comparison of wave heights measured by the Elbe buoy with numerical model results (a), and Sentinel-3A altimeter measurements (b). The model and satellite data sets were corrected according to the calibration factors λ_WAM and λ_S3 and bias parameters b_WAM and b_S3 estimated in the triple collocation procedure (see Eq. 21). (c) Direct comparison between model and satellite data. Numbers are given for the estimated calibration factors, bias, and stochastic error standard deviations ϵ_WAM, ϵ_S3, and ϵ_Buoy. The red triangles refer to ascending (asc) satellite passes and the green ones to descending (desc) passes.

Figure 6 shows the estimated biases (Fig. 6a and d) for the Sentinel-3A altimeter (Fig. 6a, b, c) and the wave model (Fig. 6d, e, d). One can see that the altimeter seems to be either bias free or slightly biased high for most of the cases (Fig. 6a). Averaging over all buoys gives a bias estimate of

$$\begin{array}{}\text{(52)}& \langle b_{\mathrm{S}\mathrm{3}}{\rangle }_{\mathrm{buoys}}\approx \mathrm{0.07}\pm \mathrm{0.31}\mathrm{m}.\end{array}$$

Concerning the spatial distribution of observation errors, it is hard to draw a conclusion. It is, however, evident that the few cases with low bias are far offshore (Fig. 6a). For the wave model there a few more cases where a small low bias is found (Fig. 6d). Again, averaging over all buoys gives

$$\begin{array}{}\text{(53)}& \langle b_{\mathrm{WAM}}{\rangle }_{\mathrm{buoys}}\approx -\mathrm{0.03}\pm \mathrm{0.26}\mathrm{m}.\end{array}$$

The spatial distribution shows a weak clustering of low-bias cases in the northern part of the North Sea. It is interesting to note that for the location of the light vessel 62170 near the east entrance of the English Channel (see Fig. 1), the satellite and the model show a positive bias of about 0.3 and 0.2 m, respectively. According to Anderson et al. (2016), one can expect a systematic low bias for wave height measurements from light vessels of about 0.3 to 0.4 m. It is thus possible that the estimated high bias for satellite and model is in this case an artefact caused by the violated assumption of bias-free in situ observations.

Figure 6 Colour coded biases (a, d) and calibration factors (b, e) derived by triple collocation for the Sentinel-3A altimeter (a, b, c) and WAM model wave heights (d, e, f). The right column (c, f) gives the uncertainties in the slope estimations derived using the approach in Sect. 2.2 as standard deviation.

The scaling parameter for the satellite altimeter shown in Fig. 6b indicates values above 1 for most of the cases. In fact, averaging over all buoys gives

$$\begin{array}{}\text{(54)}& \langle {\mathit{\lambda }}_{\mathrm{S}\mathrm{3}}{\rangle }_{\mathrm{buoys}}\approx \mathrm{1.11}\pm \mathrm{0.27}.\end{array}$$

The respective scaling parameter estimation errors derived using the approach described in Sect. 2.2 are shown in Fig. 6c. It is evident that quite a few of the cases with exceptionally high scaling values (around 1.2) are affected by large estimation errors. This is a good illustration of the added value provided by the error estimation procedure presented in this study. The corresponding scaling parameters for the WAM model shown in Fig. 6e show values which are closer to unity for the most part. The respective mean value is

$$\begin{array}{}\text{(55)}& \langle {\mathit{\lambda }}_{\mathrm{WAM}}{\rangle }_{\mathrm{buoys}}\approx \mathrm{1.02}\pm \mathrm{0.20}\end{array}$$

with higher values (around 1.1) found in the English Channel. Most of the other cases closer to the coast have slope values slightly below unity. Most of the cases with large estimation errors for the scaling factor (Fig. 6f) are found close to the coast.

Results for the stochastic errors are summarised in Fig. 7. The columns refer to Sentinel-3A (Fig. 7a and d), the WAM model (Fig. 6b and e), and the buoys (Fig. 6c and f). Figure 7a, b, and c show the estimated stochastic error standard deviation and the Fig. 7d, e, and f the respective relative estimation errors ν, defined as

$$\begin{array}{}\text{(56)}& \mathit{\nu }=\mathrm{100}\mathit{\%}{\displaystyle \frac{\mathrm{SD}(\langle {\mathit{ϵ}}^{\mathrm{2}}\rangle )}{\langle {\mathit{ϵ}}^{\mathrm{2}}\rangle }},\end{array}$$

where 〈ϵ²〉 is the error variance, and SD(〈ϵ²〉) is the standard deviation of the respective estimator, derived using the approach described in Sect. 2.1. One can see that, overall, the smallest stochastic errors are found for the buoys, as expected (Fig. 7c). In fact, one gets

$$\begin{array}{}\text{(57)}& \langle \langle {\mathit{ϵ}}_{\mathrm{Buoy}}^{\mathrm{2}}{\rangle }^{\mathrm{1}/\mathrm{2}}{\rangle }_{\mathrm{buoys}}\approx \mathrm{0.12}\pm \mathrm{0.11}\mathrm{m}\end{array}$$

averaging over all buoys. There are two buoys which stand out with errors above 0.25 m in the northern part of the North Sea. In this case, the estimation errors are not exceptionally high and possible reasons for these relatively high error levels should be further investigated. In general, one can see that the estimation errors are quite large, exceeding in most cases 20 % (Fig. 7f).

Figure 7 (a, b, c) Colour coded stochastic error standard deviations of wave heights provided by the Sentinel-3A altimeter (a), the WAM model (b), and the in situ stations (c) estimated by triple collocation. (d, e, f) Relative uncertainties in the stochastic error estimates derived by using the approach in Sect. 2.1.

The stochastic errors of the WAM model (Fig. 7b) and the altimeter (Fig. 7a) are quite similar in their average values

$$\begin{array}{}\text{(58)}& {\displaystyle }& {\displaystyle }\langle \langle {\mathit{ϵ}}_{\mathrm{WAM}}^{\mathrm{2}}{\rangle }^{\mathrm{1}/\mathrm{2}}{\rangle }_{\mathrm{buoys}}\approx \mathrm{0.17}\pm \mathrm{0.07}\mathrm{m},\text{(59)}& {\displaystyle }& {\displaystyle }\langle \langle {\mathit{ϵ}}_{\mathrm{S}\mathrm{3}}^{\mathrm{2}}{\rangle }^{\mathrm{1}/\mathrm{2}}{\rangle }_{\mathrm{buoys}}\approx \mathrm{0.18}\pm \mathrm{0.14}\mathrm{m}.\end{array}$$

It is interesting to see that the two buoys mentioned above also stand out with respect to the corresponding model errors. Theoretically, this could be due to a correlation between the background statistics and both the model and buoy errors. However, because this is observed in a quite homogeneous offshore area, with neighbouring buoys not showing the same effect, this explanation is not very likely. It is more likely that the basic assumptions about zero bias or absent calibration errors are violated for these buoys.

The finding that, on average, the in situ stations have the smallest stochastic errors is at first sight in disagreement with results presented in Janssen et al. (2007). One has to take into account, however that there are a number of significant differences in the analysis. First of all, a global wave model with 55 km resolution was used in the former study, whereas the computational model grid used in our analysis has a resolution more than 15 times higher. It is unlikely, however that the coarser model resolution is the only factor, because Caires and Sterl (2003) also concluded that the in situ stations have the smallest stochastic errors using wave model output with even coarser resolution (1.5^∘) than used by Janssen et al. (2007). Both studies introduced altimeter super-observations (averages over subsequent measurements) to make the altimeter observations more consistent with the model estimates. In the present study this was not considered necessary, because the altimeter and model resolutions are in much closer agreement. The second major difference with respect to previous studies is the geographic locations and the type of altimeter data considered in the analysis. Janssen et al. (2007) investigated global ERS-2 and Envisat altimeter data sets, while Caires and Sterl (2003) concentrated on TOPEX and ERS-1 altimeter data acquired over the Pacific and the US east coast. This means that there are certainly differences both with regard to the background wave statistics and the satellite and in situ observation errors. A third important difference between the studies is the applied collocation criteria. Janssen et al. (2007) required the model, in situ, and satellite estimates to be within 200 km distance and Caires and Sterl (2003) used a smaller collocation distance of 0.75^∘. The allowed distance of 10 km used in the present study is still significantly smaller than that, and the collocation errors are therefore also likely to be smaller. For the above reasons, one cannot conclude that the present study contradicts the results in Janssen et al. (2007). The conclusion is rather that a common set of reference in situ data and collocation criteria are desirable to make different studies more comparable.

It is evident that the observed heterogeneity of in situ measurements is a big complicating factor in the analysis. Wave model computations and satellite altimeter observations have reached a level of accuracy where further improvements require a very careful selection and treatment of validation data sets. This in particular requires more knowledge about the type of in situ instruments and applied data processing techniques (e.g. averaging intervals). This could also be an argument for investments into dedicated validation instruments with more transparent and better documented error characteristics and quality control. The deployment of such instruments should take into account both research aspects and requirements for operational use.

5 Multi-collocations

In this section different examples are presented where more than three observations are combined, i.e. this is beyond the standard triple collocation approach. The two examples discussed in the following are typical situations encountered when analysing in situ data, model data, and satellite measurements in combination.

5.1 1-D example

The geometry of the first example is depicted in Fig. 8a. Here, an ascending Sentinel-3A track passes between the two in situ stations 62150 and 62289. The station on the easterly side is within 10 km distance of the track and was therefore used in the triple collocation study presented in Sect. 4. Station 62150 did not match the criteria and was disregarded for the analysis. Both stations can be found in Figs. 1 and 2b, where they are indicated by triangle symbols.

Figure 8 (a) Example of observation configuration with a descending Sentinel-3A altimeter track passing between two in situ wave observation instruments (62150 and 62289), (b) Observation geometry with Sentinel-3A track passing through a group of three in situ wave measurement devices. The blue and red altimeter measurements are used to estimate error correlations for both the altimeter and the numerical model.

The idea to relate both in situ measurements to the altimeter track is to use a linear interpolation of the truth wave height between the two stations, which makes the use of the instrument with the larger distance more acceptable in the collocation procedure. In principle, this corresponds to the 1-D case depicted in Fig. 3b with the role of altimeter and model interchanged. In the present situation there is one altimeter measurement between the two reference instruments and for simplicity the numerical model wave height estimate is taken at the location of the buoys. Because of the small number of available samples, we have also used altimeter measurements which are slightly above and below the connecting line (red dots in Fig. 8a). This resulted in n_s=14 common data samples that could be used for the statistics.

Using this geometry allows for the estimation of the errors of all data sources, as well as the error correlations between the model wave heights (see Table 1). The calibration factors and their respective standard deviations were estimated with the direct and iterative method and are as follows.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\lambda }}_{\mathrm{62150}}^{\mathrm{WAM}}=\mathrm{0.662}\pm \mathrm{0.147}(\mathrm{0.788}\pm \mathrm{0.161})\\ \text{(60)}& {\displaystyle }& {\displaystyle }{\mathit{\lambda }}_{\mathrm{62289}}^{\mathrm{WAM}}=\mathrm{0.779}\pm \mathrm{0.113}(\mathrm{0.778}\pm \mathrm{0.100}){\displaystyle }& {\displaystyle }{\mathit{\lambda }}_{\mathrm{S}\mathrm{3}}=\mathrm{1.023}\pm \mathrm{0.246}(\mathrm{1.023}\pm \mathrm{0.360})\end{array}$$

The values in brackets were obtained with the iterative method. Significant differences are only found for the first scaling parameter. However, both methods provide consistent results if the error bars are taken into account. It is interesting to note that the scaling value for Sentinel-3A is closer to unity than the smaller value of about 0.8 found by the triple collocation method (Fig. 6b). This value was exceptional among the other buoys for which numbers above 1 were found for the most part. This could be an indication of a problem with station 62289, which also stands out in the stochastic errors shown in Fig. 7c.

The numbers obtained for the stochastic errors are as follows.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathrm{var}\left({\mathit{ϵ}}_{\mathrm{62150}}^{\mathrm{Buoy}}\right)=-\mathrm{0.0890}\pm \mathrm{0.0914}{\mathrm{m}}^{\mathrm{2}}(-\mathrm{0.0889}\pm \mathrm{0.0781}{\mathrm{m}}^{\mathrm{2}})\\ {\displaystyle }& {\displaystyle }\mathrm{var}\left({\mathit{ϵ}}_{\mathrm{62289}}^{\mathrm{Buoy}}\right)=-\mathrm{0.0072}\pm \mathrm{0.0234}{\mathrm{m}}^{\mathrm{2}}(-\mathrm{0.0072}\pm \mathrm{0.0235}{\mathrm{m}}^{\mathrm{2}})\\ {\displaystyle }& {\displaystyle }\mathrm{var}\left({\mathit{ϵ}}_{\mathrm{62150}}^{\mathrm{WAM}}\right)=\mathrm{0.0749}\pm \mathrm{0.0467}{\mathrm{m}}^{\mathrm{2}}(\mathrm{0.0913}\pm \mathrm{0.0557}{\mathrm{m}}^{\mathrm{2}})\\ {\displaystyle }& {\displaystyle }\mathrm{var}\left({\mathit{ϵ}}_{\mathrm{62289}}^{\mathrm{WAM}}\right)=\mathrm{0.0234}\pm \mathrm{0.0167}{\mathrm{m}}^{\mathrm{2}}(\mathrm{0.0234}\pm \mathrm{0.0167}{\mathrm{m}}^{\mathrm{2}})\\ {\displaystyle }& {\displaystyle }\mathrm{covar}({\mathit{ϵ}}_{\mathrm{62150}}^{\mathrm{WAM}},{\mathit{ϵ}}_{\mathrm{62289}}^{\mathrm{WAM}})=\mathrm{0.0242}\pm \mathrm{0.0095}{\mathrm{m}}^{\mathrm{2}}(\mathrm{0.0241}\pm \mathrm{0.0095}{\mathrm{m}}^{\mathrm{2}})\\ \text{(61)}& {\displaystyle }& {\displaystyle }\mathrm{var}\left({\mathit{ϵ}}_{S\mathrm{3}}\right)=\mathrm{0.1372}\pm \mathrm{0.0555}{\mathrm{m}}^{\mathrm{2}}(\mathrm{0.1372}\pm \mathrm{0.0550}{\mathrm{m}}^{\mathrm{2}})\end{array}$$

It can be seen that the estimates for the buoys are slightly negative, which is not meaningful for a variance. This can in fact happen for small sample sizes, since the estimators do not guarantee positive values. In this case it is helpful to look at the respective error bars, which are given as standard deviations. For a Gaussian-distributed variable the interval given by ± SD gives a 68 % confidence interval, i.e. more than 30 % of the cases are outside of this value range. This means that the estimated values for buoys are consistent with small positive error variances. The largest value is found for the Sentinel-3A altimeter with a relatively small error bar. This is consistent with the finding already made with the triple collocation method (see Fig. 7a). For the WAM model at the location of the 62289 station, the triple collocation method gave a similarly high value, but with almost 100 % error margin. The estimate obtained with the multi-collocation is significantly lower, but again with a large relative estimation error of about 100 %. Because of the smaller mean value, the latter estimates still point towards a smaller model error than indicated by the triple collocation method.

The covariances estimated for the WAM wave height errors at the two buoy locations correspond to a correlation value of 0.58. If we assume that the error autocorrelation function is Gaussian shaped, i.e.

$$\begin{array}{}\text{(62)}& \mathrm{ACF}\left(\mathrm{\Delta }x\right)=\exp (-{\displaystyle \frac{\mathrm{\Delta }x}{{\mathit{\lambda }}_{\mathrm{C}}}})\end{array}$$

with correlation length λ_C and spatial distance Δx, the above value results in λ_C=55 km.

Because of

$$\begin{array}{ll}{\displaystyle }\langle ({\mathit{ϵ}}_{\mathrm{62150}}^{\mathrm{WAM}}-{\mathit{ϵ}}_{\mathrm{62289}}^{\mathrm{WAM}}{)}^{\mathrm{2}}\rangle & {\displaystyle }=\mathrm{var}\left({\mathit{ϵ}}_{\mathrm{62289}}^{\mathrm{WAM}}\right)+\mathrm{var}\left({\mathit{ϵ}}_{\mathrm{62150}}^{\mathrm{WAM}}\right)\\ \text{(63)}& {\displaystyle }& {\displaystyle }-\mathrm{2}\mathrm{covar}({\mathit{ϵ}}_{\mathrm{62289}}^{\mathrm{WAM}},{\mathit{ϵ}}_{\mathrm{62150}}^{\mathrm{WAM}}),\end{array}$$

knowledge about the variances and covariances also allows an estimate of the uncertainties in the gradient estimates. In this case an error standard deviation of 0.31 m was obtained for the difference of the WAM model wave heights at the two buoy locations.

5.2 2-D example

The geometry of the second example is depicted in Fig. 8b. This is an area in the northern part of the North Sea around 58^∘ latitude between England and Norway. In this case an ascending Sentinel-3A track is passing through a group of three in situ wave observation platforms, which are shown in Figs. 1 and 2b. Here, we concentrate on two locations covered by the satellite, which appear as two clusters in Fig. 8b. The “North” group of satellite observations is shown in blue and the “South” group in red. Including the numerical model estimates at those locations, the situation is then as described by the last row in Table 1. One has 7 wave height estimates in total, and a 2-D plane approximation is used for the observed area. The multi-collocation method then allows an estimation of the errors of all components, as well as three covariances. As the buoy measurements can be assumed as independent, only two covariances are required in this example; this is the covariance between the model errors at the two locations and the same for the altimeter measurements. With this configuration the number of available data sets was n_s=11.

The scaling values and their standard deviations obtained with the direct method are as follows.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathit{\lambda }}_{\mathrm{WAM}}^{\mathrm{South}}=\mathrm{1.130}\pm \mathrm{0.006}\\ {\displaystyle }& {\displaystyle }{\mathit{\lambda }}_{\mathrm{WAM}}^{\mathrm{North}}=\mathrm{1.104}\pm \mathrm{0.004}\\ \text{(64)}& {\displaystyle }& {\displaystyle }{\mathit{\lambda }}_{\mathrm{S}\mathrm{3}}^{\mathrm{South}}=\mathrm{1.270}\pm \mathrm{0.002}{\displaystyle }& {\displaystyle }{\mathit{\lambda }}_{\mathrm{S}\mathrm{3}}^{\mathrm{North}}=\mathrm{1.272}\pm \mathrm{0.003}\end{array}$$

Here, the values labelled “North” refer to the northern cluster of Sentinel-3A measurements (blue points in Fig. 8b) and the values labelled with “South” refer to the southern group of observations (red points in Fig. 8b). These estimates seem to be quite robust, because of the small error bars and the fact that the errors in both areas are very similar. The results also confirm the overall finding of the triple collocation analysis which indicated a wave height overestimation by the Sentinel-3A altimeter.

The respective values for the stochastic errors and their standard deviations with the same naming convention and obtained with the direct method are as follows:

$$\begin{array}{ll}{\displaystyle }\mathrm{var}\left({\mathit{ϵ}}_{\mathrm{62168}}\right)& {\displaystyle }=\mathrm{0.003}\pm \mathrm{0.007}{\mathrm{m}}^{\mathrm{2}}\\ {\displaystyle }\mathrm{var}\left({\mathit{ϵ}}_{\mathrm{62161}}\right)& {\displaystyle }=\mathrm{0.010}\pm \mathrm{0.006}{\mathrm{m}}^{\mathrm{2}}\\ {\displaystyle }\mathrm{var}\left({\mathit{ϵ}}_{\mathrm{62134}}\right)& {\displaystyle }=\mathrm{0.014}\pm \mathrm{0.007}{\mathrm{m}}^{\mathrm{2}}\\ {\displaystyle }\mathrm{var}\left({\mathit{ϵ}}_{\mathrm{WAM}}^{\mathrm{South}}\right)& {\displaystyle }=\mathrm{0.016}\pm \mathrm{0.008}{\mathrm{m}}^{\mathrm{2}}\\ \text{(65)}& {\displaystyle }\mathrm{covar}({\mathit{ϵ}}_{\mathrm{WAM}}^{\mathrm{North}},{\mathit{ϵ}}_{\mathrm{WAM}}^{\mathrm{South}})& {\displaystyle }=\mathrm{0.009}\pm \mathrm{0.005}{\mathrm{m}}^{\mathrm{2}}{\displaystyle }\mathrm{var}\left({\mathit{ϵ}}_{\mathrm{WAM}}^{\mathrm{North}}\right)& {\displaystyle }=\mathrm{0.005}\pm \mathrm{0.003}{\mathrm{m}}^{\mathrm{2}}\\ {\displaystyle }\mathrm{var}\left({\mathit{ϵ}}_{\mathrm{S}\mathrm{3}}^{\mathrm{South}}\right)& {\displaystyle }=\mathrm{0.011}\pm \mathrm{0.007}{\mathrm{m}}^{\mathrm{2}}\\ {\displaystyle }\mathrm{covar}({\mathit{ϵ}}_{\mathrm{S}\mathrm{3}}^{\mathrm{North}},{\mathit{ϵ}}_{\mathrm{S}\mathrm{3}}^{\mathrm{South}})& {\displaystyle }=\mathrm{0.005}\pm \mathrm{0.005}{\mathrm{m}}^{\mathrm{2}}\\ {\displaystyle }\mathrm{var}\left({\mathit{ϵ}}_{\mathrm{S}\mathrm{3}}^{\mathrm{North}}\right)& {\displaystyle }=\mathrm{0.012}\pm \mathrm{0.007}{\mathrm{m}}^{\mathrm{2}}\end{array}$$

Due to the significant estimation errors it is hard to tell which data source has the smallest errors. The obtained numbers are consistent with an error standard deviation of around 0.1 m for all data sets. The error estimates for the altimeter at the two locations agree very well and are also consistent with the values found with the triple collocation method (Fig. 7a). The difference in the error variances for the WAM model at the two locations appear to be quite big considering the distance of about 30 km. But again, the error bars show that there is a significant probability that the errors are actually in closer agreement. In principle, it would be possible to force the WAM error variances at the two locations to be the same, using a respective formulation of the linear system Eq. (12). However, looking at the spatial variations in the bathymetry in Fig. 1, this is hard to justify.

For the correlation, a value of 0.39 was found for the altimeter and a value close to 1 for the WAM model. This corresponds to a correlation length of about 30 km for the satellite data. It makes sense that the correlation length for the WAM model is longer in this case compared to the configuration discussed in the previous section, because the analysed area is in deeper water quite far offshore, and can therefore be assumed as more homogeneous with respect to model errors.

The examples show that the multi-collocation method is in fact applicable to real data source configurations. In particular, the matrix D in Eq. (12) is regular for the considered geometries, and estimates for error correlations can be obtained. It is also evident, of course, that the limited number of samples results in significant estimation errors. According to Eq. (19), the variance of the error variance estimation scales with 1∕n_s, i.e. in order to reduce the error bars given in Eqs. (61) and (65) by a factor of 2, the number of samples has to be increased by a factor of 4.

6 Conclusions

The presented study provides an extension of the known triple collocation method, which can be useful in areas with stronger gradients, like coastal regions, where nearest-neighbour approximations maybe critical. The proposed method is very flexible in the way that various parameterisations can be used to describe the spatial variability of the measured quantities. In this study we considered only linear models, but this is not a restriction of the method, since more sophisticated functional forms (e.g. bilinear functions) can be easily integrated. Such higher-order approaches are certainly desirable for coastal areas with strong spatial variations; however, they require a larger number of data sources (see Eq. 9).

An approach was proposed to estimate the uncertainties in estimated calibration and stochastic errors, which is also useful in the context of the standard triple collocation method, which is a special case of the multi-collocation technique. The technique uses the covariance matrices of the input data and the number of samples as input, i.e. bootstrapping is not required. These uncertainty estimates are seen as very valuable, in particular in the context of new instruments for which only a limited data set is available for the assessment.

The proposed techniques were validated using Monte Carlo simulations with realistic background statistics. It was shown that the obtained error estimates and their respective uncertainties are in good agreement with the expected values, although a couple of approximations had to be used in the derivation.

The method was applied to a data set of in situ wave measurements, Sentinel-3A altimeter observations, and numerical wave model data. The number of available samples was relatively small and estimation errors had therefore to be taken into account. The usefulness of the derived error bars for the interpretation of the data could be demonstrated. For the analysed 16 months data set presented here, the estimation errors are significant, in particular if individual geographic locations are analysed. It would therefore be interesting to continue some parts of the analysis at a later stage of the Sentinel-3A mission, when a larger data set will be available. More robust results are obtained if the systematic and stochastic data set errors estimated for different in situ instrument locations are averaged. The results obtained for the North Sea indicate the smallest stochastic errors for the in situ measurements, as expected. The stochastic errors of the model and the altimeter seem comparable if averaged over all in situ locations. The analysis indicates that on average the altimeter is overestimating wave heights by about 10 % for above-mean wave conditions. Two examples of multi-collocations were analysed, which included a group of two and three in situ platforms. In both cases a Sentinel-3A track passed nearby, and model data were used in addition. The use of 1-D and 2-D parameterisations for the first and second example, respectively, resulted in estimates for the spatial decorrelation of model and altimeter errors.

The proposed method can be used for many other applications not discussed in this study. For example, it is straightforward to extend the analysis of error correlations to the time domain. The method can also be applied in situations where different instruments do not measure exactly the same quantity, but different components of a truth vector, for example HF radar providing 2-D current vectors and satellite SAR providing one current component (e.g. Hansen et al., 2011).

This study is supposed to make a contribution to the optimal use of the growing number of observations, in particular in coastal areas. For applications, like data assimilation, knowledge about the errors of different data sources is essential. Analysis of observation errors is also a critical component in the design and extension of observatories used for various applications. This subject will be of growing concern, for example, in the context of the European marine core service (CMEMS), where in situ data are required to optimise forecasts for all European seas.