On the combination of altimetry and tide gauge observations in the Norwegian coastal zone

1 Introduction

Lately, the Norwegian Mapping Authority (NMA) has seen a need (Kartverket, 2019, p. 6) to clarify and improve connections between height reference frames in the coastal zone. Consequently, the NMA can now offer a unified framework for planning (https://www.kartverket.no/api-og-data/separasjonsmodellar) on land and at sea within the territorial border of Norway. Essential in the development of this administrative tool is knowledge of the mean sea level (MSL).

Today’s official description of MSL by the NMA is determined by sea level measurements at tide gauges along the Norwegian coast. One can therefore expect a decline in the model’s performance with increasing distance from the mainland. Further from the coast, however, the sea level is observed by satellites. We therefore want to investigate how information from altimetry missions can be included in the services of the NMA. For our work, we choose the ESA satellites Sentinel-3A/B since we are familiar with their use.

Hopefully, our work will shed some light on the challenges of a seamless tide gauge connected, altimetry-based MSL model covering the entire Norwegian Economic Zone. On the open sea, a reliable model can be determined by altimetry, but must be consistent with datum defining tide gauges as these are approached. If combined properly, the area to which the NMA models can be applied with confidence increases manyfold, from within the outer black line segment (territorial border) of Figure 1 to within the green border (Norwegian Economic Zone).

Figure 1

Permanent tide gauges are marked in blue, while red shows the locations of the temporary tide gauges. Special areas of interest are also highlighted: 1. Finnmark, 2. Trøndelag/Nordland, 3. Vestlandet, 4. Sør-Norge. The baseline and territorial borders are given by the inner and outer black line segments, respectively, while the perimeter of the Norwegian Economic Zone is drawn in green.

Previous works that utilize both tide gauge and altimetry observations include the studies by Iliffe et al. (2007), Slobbe et al. (2014), Robin et al. (2016). Specifically, Slobbe et al. (2014) point out that nontrivial spatial interpolation is needed to close the data gap between altimetry and tide gauges along the coastline. In the following, we take a closer look at this interpolation problem using robust geostatistical methods applied in a moving window approach. When appropriate, it has several benefits.

First, the models that need to be considered are very simple compared to global models, which must describe not only the region around the point of interest but also an area covering all data points. Second, predictions can be made using only observations in the vicinity, which is very convenient if altimetry is to be included. As a result, matrices that must be inverted are small. Third, the high speed and simplicity of the method allow for effective tests of what are our ultimate goals, to make good predictions at locations where observations are not available, and to say something about how uncertain these predictions are.

This cross validation can be accomplished by removing one or more observations in the region of focus, and comparing these with predictions based on other nearby observations. Notably, with our approach, the computed uncertainties will be representative on a local scale. This cannot be expected from many global models, like the current official one by the NMA, due to how their characteristics are derived from all available observations.

To be more specific, a fundamental thought behind our study is to predict MSL at the locations of the tide gauges, with the use of nearby altimetry-derived MSL. The deviation between predicted and actual values at the gauges are then to be compared with estimates of the deviation’s uncertainty, constituting a type of validation that is applied to evaluate consistency (e.g., Lark 2000, Mysen 2014, 2020). If the deviations correspond statistically to the computed uncertainties, a representative model for the two different data types of the same quantity is found. As a result, we can have faith in what the model tells us about the relative weight of tide gauge measurements and altimetry observations in regions between data points.

2 Data types

The determination of MSL will, in this work, be partly based on observations from permanent and temporary tide gauges along the Norwegian coast, see Figure 1. The quantity MSL is defined as the average sea level over the period 1996–2014^[1] and has been observed directly at the 26 permanent gauges with an uncertainty of about 5 mm (Kartverket, 2020, p. 33).

At 151 locations, temporary gauges have been placed to yield observations which, when compared sea levels at nearby permanent gauges, can be used to find MSL there. The accuracies lie between 1 and 4 cm, but their specific values are not so important in our work focused on MSL outside the Norwegian baseline, encompassing the outer islets. To facilitate a comparison with altimetry, we also correct the tide gauge observations for atmospheric pressure^[2].

Estimates of MSL can also be derived from observations made by the ESA satellites Sentinel-3A (2016–2022) and -3B (2018–2022) and are available along (crossing) ground tracks separated by 15–25 km at our latitudes. The altimetry-derived sea levels are grouped in zones with radius 3.5 km and averaged, after removing known temporal and general annual and semiannual variations, to one point for each zone, see Figure 2. These points will sometimes be referred to as altimetry observations (of MSL) in the following. To avoid land reflected radar pulses, we do not use sea levels computed this way which are located closer to land than 10 km or in the fjords. Finally, we add a global constant to these estimates, so that they are aligned on average with the datum of the permanent tide gauges.

Figure 2

The prediction point of interest is here the tide gauge at Kjøllefjord, Finnmark (circled red box). Locations associated with altimetry-based MDT are given by the black points, with the point nearest to the tide gauge marked by a red circle. The values at locations (small circles) within 75 km of the red circle are used to compute the variogram, in this case consisting of 157 samples.

Three other MSL surfaces are also considered in our work, namely, NMBU18 (Ophaug et al. 2021) based on recent generations of altimetry satellites, DTU21 (Andersen et al. 2023) with emphasis on older missions, and IC2-MSS (Tomic and Andersen 2023) from ICESat-2. The IC2-MSS based on laser altimetry will be brought up when model performance in the deep Sognefjorden is evaluated.

While MSL relative to the ellipsoid decreases by more than 25 m from south to north along the Norwegian coast, the height of the MSL above the geoid changes far less, and is therefore easier to handle for some applications. In our case, we subtract the gravimetric NKG2015 geoid (Ågren et al. 2016) from our MSL values, both converted to the mean-tide system (e.g., Ekman 1989), to obtain estimates of the mean dynamic topography (MDT). The accuracy of this MDT is dependent not only on the quality of the MSL estimates but also on the subtracted geoid. From comparisons with precise GNSS/levelling on the mainland, we find the error of the NKG2015 geoid to increase from 0 to 4 cm over spatial separations of about 300 km (Kartverket 2021).

The NMA already offers an official model (https://www.kartverket.no/api-og-data/separasjonsmodellar) for heights in the coastal zone, based on tide gauge observations. The previous work (Breili 2023) shows that the official model corresponds fairly well with Sentinel-3 derived MDT at the edge of the model’s validity, along the territorial border of Norway, 12 nautical miles ( ∼ 22 km ) outside Norway’s baseline. Here, we take the comparison a step further by focusing on how to combine the two different information sources on MDT.

3 Nearest neighbour

To exploit the available information, a problem is often sought solved using as few unknown parameters as possible. We can, for instance, subtract other MSL surfaces, say DTU21, from our Sentinel-3 based MSL estimates, and try to describe this new and far less variable quantity. The success of this proposal will be, for interpolation purposes, highly dependent on the stability of the systematic errors of DTU21 as land is approached.

We will, instead, focus on the MDT, i.e., the MSL above the geoid, as target quantity, which changes by about 20 cm only (e.g., Ophaug et al. 2021) along the Norwegian coast. Our results therefore obtain a more direct physical interpretation, for future reference.

Accounting for observation noise n i , we have samples z i of MDT along the Sentinel-3 tracks and at the tide gauges

where the index i ≥ 1 indicates the sample’s georeference. The index i ≡ 0 we reserve for prediction points, which may or may not coincide with the location of an observation.

In the following, we set

for both data types. In the case of our Sentinel-3 based MDT, we see little indication of spatially uncorrelated patterns. As for the tide gauge observations, equation (2) is not descriptive. Within the framework we are to introduce, however, a nonzero observation noise at an isolated gauge determines the uncertainty in its surroundings only. Away from the gauges and outside the Norwegian baseline, which we are mainly interested in here, observation noise will be of less importance merely because the MDT has changed significantly from what has been measured at the gauge.

At terms with our previous discussion, we will attempt to describe the MDT in a small region by a local constant x 1

with s i as the signal, often viewed as a realization of a spatial stochastic model. Independent of how suitable equation (3) is, but certainly motivated by it, we start by letting the nearest value z i define the prediction

This simple and transparent estimator is robust in the sense that it is bounded by a nearby (more or less) directly observed value of MDT in the vicinity of the prediction point. Eventually, we will allow for more than one observation to define the prediction. But, for our preliminary studies close to the coast, the inclusion of more observations will lead to technicalities that are beyond the scope of our investigations. From least squares collocation (Moritz 1980), the nearest neighbour estimate equation (4) follows from equation (3) if only one observation is available.

The uncertainty σ 0 of the nearest neighbour estimator is given by

where E ( ⋅ ) is the statistical expectancy, which we parameterize by the spatial separation Δ between z i and z 0 . This does not mean that we assume isotropy, but that the uncertainty estimates should be considered as averages with respect to direction. Now, we will take a closer look at the variograms, equation (5).

4 Variograms

To evaluate the local degree of variation in MDT, we will make use of the empirical variogram 2 γ ˆ derived solely from altimetry-based MDT values. That is, we first find the location of the altimetry-based MDT sample that is nearest to the point of prediction. From this central sample, and altimetry-based MDT samples within 75 km of it (Figure 2), we compute all unique differences z j − z k . These are then grouped according to their spatial separation Δ j k to get^[3]

with n Δ as the number of unique pairs separated by the distance Δ ≠ 0 . The variogram itself (Matheron 1963) is given by

based on equations (1)–(3), providing the nearest neighbour uncertainty, equation (5).

We determine variogram’s defining parameters α by minimizing

for a chosen function 2 γ ( α , Δ ) , up to a maximum spatial separation Δ t = 0.5 Δ max ∼ 75 km . The weights

are from Cressie (1985), related to the variogram’s dispersion. The minimization is performed by the downhill simplex method^[4] (e.g., Press et al. 2006, p. 402).

Due to the variogram’s definition as an expectation (or average) of squared increments, our choices are restricted to conditionally negative definite (e.g., Cressie 2015, p. 60, Chiles and Delfiner 2012, p. 63) function types. Since we compute Δ in plane UTM33 coordinates, we have that

is valid (e.g., Cressie 2015, p. 61, Chiles and Delfiner 2012, pp. 31 and 88). The exponential form is motivated by empirical variograms along the Norwegian coast, see Figure 3 based on altimetry derived MDT samples outside Kjøllefjord.

Figure 3

The solid line represents equation (10) fit to the altimetry-based empirical variogram (points) outside Kjøllefjord, Finnmark. The assumed dispersion derived from equation (9) is shown by the dashed lines.

We see, however, that a linear increase with distance, as shown in Figure 4 for Rong on the west coast of Southern Norway, is more common.

Figure 4

The altimetry-based variogram for the tide gauge near Rong, in the Øygarden municipality off the west coast of Southern Norway, is given together with a fit (solid line).

But, as the figure demonstrates, the proposal equation (10) can easily adapt by fixing α 2 to a value that is considerably larger than typical interpolation distances Δ ≲ Δ t . For the minimization of equation (8) to converge in general, it is therefore important that the correlation length α 2 is not allowed to attain values above a defined maximum, here 300 km. In the applied downhill simplex algorithm, such constraints are trivial to implement since it is solely based on evaluations of the function to be minimized, and not its derivatives. That is, we can let 2 γ return a very high value if α 2 > 300 km .

It must be mentioned that it is necessary to add a small fixed quantity, its specific value not important, to equation (10) for the fits to converge. Typical values for the linear variograms, which we find along the whole coast of mainland Norway, lie in the range α 1 ∕ α 2 ∼ 0.3 – 2 cm 2 ∕ km .

A strong linear deterministic trend in our observations will be visible as a quadratic increase (Cressie 2015, p. 72, Bliznyuk et al. 2012) in the variogram, which we sometimes see indications of, especially for the larger spatial separations Δ ≳ 75 km considered. We do not rule out that our linear variograms can be caused by couplings between spatially correlated signals and weaker deterministic trends.

5 Validation

5.2 Observations: tide gauges

Up until now, we have only used altimetry derived MDT to make predictions and estimates of uncertainty at the location of tide gauges. Table 1, however, demonstrates the performance of our approach mostly in the direction towards the coast, and our impression so far may be too optimistic. To check this, we will now make predictions more along the coast, by using tide gauge observations as the nearest neighbour equation (4). As always in this work, the variogram is computed using altimetry derived MDT nearest to the point of prediction, see Figure 2.

Close to the coast, it is more important that Δ of equation (10) is computed as the sea distance. Representing an extreme case, we have the temporary tide gauge at Lestoneset shown in Figure 5. Although the nearest tide gauge is only 3.6 km away in UTM33 coordinates, their separation by sea is 67 km. Instead, the nearest sea connected tide gauge is at Osmundsvåg 15 km away. To obtain these sea paths, we use the Dijkstra algorithm (Dijkstra 1959)^[5], where we found it necessary with a sea grid spacing as small as 50 m to connect all 177 tide gauges to the open sea. Smoother paths can be obtained by increasing the resolution further, and allowing each node in the grid to connect directly with more distant nodes, but this comes at a great computational burden.

Figure 5

The temporary tide gauge at Lestoneset is marked by a circle. Its nearest other tide gauge in sea distance is at Osmundsvåg (solid red line), while the distance to Kjøde is shortest measured in a Euclidean metric.

By using sea distances, we obtain the results of Table 2, by predicting (comparing) the MDT at tide gauges from the MDT at the nearest tide gauge. Accounting for uncorrelated observation noise at the gauges, which we have not done in our work, will always lower RMS( τ ), and especially so if the absolute deviations are small. This could explain the too optimistic uncertainties we see for the shortest interpolation range from 0 to 10 km. The, in other cases, smaller RMS( τ ) indicate that our estimates of uncertainty can be somewhat conservative, but acceptable also for this type of interpolation.

Table 2

The RMS of δ and τ of equations (11) and (12) are given, obtained by comparing the MDT at each tide gauge with the MDT at a different nearby tide gauge (nearest neighbour prediction). In the first columns, the results are listed according to the sea distance Δ between the two gauges, with the number of samples in each range category in parenthesis. The last four columns show the corresponding results for the regions marked in Figure 1. Tide gauges farther away from altimetry-based MDT than 50 km are not included among the samples

Δ ( # ) ∣ ∣ zone ( # )	0–10 km (45)	10–20 km (37)	20 km + (46)	1. (14)	2. (9)	3. (39)	4. (30)
RMS ( δ ) [ cm ]	2.8	2.9	4.2	4.2	3.7	2.4	4.0
RMS ( τ )	1.14	0.73	0.87	0.76	0.79	0.95	1.15

6 Kriging

The nearest neighbour estimator equation (4) and its accuracy equation (5) can, now that we have a handle on the underlying model, be strengthened with more observations. The estimates then become smoother and more robust, and improvements can also be expected in some areas.

A framework that is often applied for this purpose is least squares collocation (LSC) (e.g., Moritz 1980, Ligas 2022). LSC, though, requires an a priori signal covariance which, at the numerous locations along the coast having approximately linear variograms, is not obviously given by the exponentially correlated form implied by equation (10). Indeed, in such cases the nonstationary Brownian field (e.g., Cressie 2015, p. 68, Stein 2002) is a more appropriate choice. To avoid the supposed need to switch between forms at different locations, we turn to a type of kriging that only requires the variogram of the signal to be well defined. For other alternatives, we refer to Ouassou et al. (2015) and references therein.

The kriging and LSC estimators are derived in very much the same way, for instance, as a linear combination of n (nearby) observations

where the weights λ i are to be found. Let the local constant x 1 of equation (3) be defined by E ( s i ) = 0 , where E ( ⋅ ) is the statistical expectancy, as mentioned earlier. From the unbiasedness requirement E ( z ˜ 0 − z 0 ) = 0 , we obtain the constraint

The accuracy of the estimator is given by

with γ i j ≡ γ ( Δ i j ) , see equation (10), where we have used equation (15) and the identity (e.g. Chilès and Delfiner 2012, p. 64)

The maximization of accuracy under the constraint equation (15) is accomplished by the minimization of

with respect to λ i , where μ is an unknown Lagrange multiplier. The solution is given by the system (e.g., Cressie 2015, p. 153–156, Chiles and Delfiner 2012, p. 171)

where γ is an n × n matrix determined by the variograms γ i j between the observations. The elements of the column vector γ 0 is given by the variogram γ 0 i between the prediction point and the observations, while 1 n is a column vector of length n that only contains 1’s. The superscript T signifies that the matrix or vector is transposed.

If only one observation z i is included, equations (14) and (17) yield the nearest neighbour estimator equation (4) and its uncertainty σ 0 2 = 2 γ i 0 , see equations (3), (5), and (7). In any case, we note that λ i = 1 and λ j ≠ i = 0 at the location of an observation z i .

As for questions of long range stability, we multiply the solution of equation (17) with the row vector ( z T , 0 ) to obtain the prediction on the form

where the vector b and scalar e are given by the system (e.g., Dubrule 1983)

Clearly, if the solution of equation (17) exists and is unique, the same applies to the solution of equation (19). To proceed, we define Δ γ i c ≡ γ i 0 − γ c 0 with c indexing some (central) location in the data region. From the last row of equation (19), we then obtain

which remains bounded as the distance between prediction and data points increases, provided Δ γ i c converges. The latter is the case for exponential and linear variograms we have mentioned in this work. A discussion of stability given two observations is provided by Chiles and Delfiner (2012, pp. 165–167).

Equation (17) can, for instance, be applied to increase the number of nearby altimetry-derived MDT samples from 1 to 10 in the prediction of MDT at a tide gauge, which yields the results of Table 3. Compared to Table 1, computed on the basis of the nearest neighbour only, we observe that the absolute performance is, for the most part, somewhat improved. Also, the statistical consistency for the Finnmark region (1) is much better, at some cost for the Trøndelag/Nordland (2) result. We have previously reported that our predictions are particularly biased in this latter case, which indicates suboptimal modelling, and could therefore explain why our uncertainties are too optimistic here.

The impact of optimistic uncertainties in the computation of an MDT model is difficult to anticipate and will be dependent on, for instance, the geometry of the observations and the nature of the adopted and true variograms. But, if this optimism is due to an estimated correlation length that is too long, our interpolation between data points would possibly be too smooth (see Chilès and Delfiner, 2012, pp. 165–167).

As previously mentioned, if predictions that also include the effect of tide gauge observations are sought close to the coast, we should use sea paths to compute distances. Equation (10) with Δ as the spatial separation by sea, however, may then no longer be conditionally negative definite, which could result in serious loss of precision and negative uncertainty estimates. Solutions to this problem (Løland and Høst 2003, Davis and Curriero 2019) will not be explored here since we are mainly interested in the region outside the Norwegian baseline, where sea paths to a greater extent are given by our adopted Euclidean metric.

We include the kriged MDT and its uncertainty off the west coast as Figures 6 and 7, respectively, based on the 10 tide gauge and the 10 altimetry derived MDT samples nearest to each prediction point. The figures illustrate significant variability, but we stress that our work is concerned with interpolation only. That is, Figure 6 may contain errors due to, for instance, imperfect corrections of the altimetry observations in the computation of MSL. Such errors in the computed MDT surface are not properly accounted for in Figure 7. Indeed, the physicality of the strong and distinct ocean currents implied by the rapid change in the MDT of Figure 6 can be questioned, see the study by Idzanovic et al. (2017).

Figure 6

MDT estimates off the coast of Norway are illustrated, based on tide gauge and altimetry observations. Values are provided up to 80 km outside the territorial border (outer set of black line segments), which is speculative since we have not studied the variograms so far out, and up to 10 km inside the baseline (inner set of black line segments). The locations of tide gauges are marked by boxes (blue: permanent, red: temporary).

Figure 7

The interpolation uncertainty of the MDT estimates of Figure 6 is shown.

The locations of altimetry derived MDT is clearly discernible in Figure 7, by the corresponding reduction in uncertainty. As for the change of this reduction with spatial separation, it will be dependent on the local variability of the MDT. For instance, the uncertainty increases rapidly near Bjørnafjorden in the lower part of Figure 7, while variograms at the entrance of Sognefjorden (upper part) tells us that the MDT changes less there. Still, as can be seen from the figure, the MDT between baseline and territorial border is mostly^[6] determined by altimetry data. The potential of consistently adding other missions with different ground tracks or higher resolution should also be evident.

7 Discussion

Overall, we see that the local isotropic variogram for the MDT increases linearly, or more slowly, with spatial separation. As a result, the optimal^[7] predictor remains bounded outside the data region and improves somewhat upon the nearest neighbour estimator.

Cross validation over a range of distances shows, for the most part, that our estimated uncertainties are reliable. This type of consistency gives us some confidence in how tide gauge and altimetry observations are weighted in a final prediction. As an example, we have mapped the interpolation uncertainties off the west coast of Southern Norway. From the map, we see that altimetry can potentially provide much improved information on MDT in the zone between the Norwegian baseline and territorial border, which is the main focus of our work.

We consider it straightforward to include the entire Norwegian Economic Zone in our computations, extending up to 200 nautical miles ( ∼ 370 km) beyond the baseline, but this generalization of our approach will be associated with considerable uncertainty. For instance, we need to establish that the local variograms far from the coast are similar to those we have encountered in this work.

Indeed, even close to the coast, we see a discrepancy between tide gauge and altimetry-derived MDT in the Trøndelag/Nordland region. Although the bias may have several causes, it indicates that our modelling is suboptimal there, which is at terms with the centred predictions of the trend accommodating NMBU18 MDT. Alternatively, the difference in performance may be related to the different ground tracks of Sentinel-3, which are relevant in our work, and those of CryoSat-2 and SARAL/AltiKa, also included in the computation of NMBU18.