## Abstract

As part of the U.K. Hydrographic Office (UKHO)-sponsored Vertical Offshore Reference Frames (VORF) project, a high-resolution model of lowest astronomical tide (LAT) with respect to mean sea level has been developed for U.K.–Irish waters. In offshore areas the model relies on data from satellite altimetry, while in coastal areas data from a 3.5-km-resolution hydrodynamic tide-surge model and tide gauges have been included. To provide for a smooth surface and predict tidal levels in unobserved areas, the data have been merged and interpolated using the thin plate spline method, which has been appropriately tuned by an empirical prediction test whereby observed values at tide gauges were removed from the solution space and surrounding data used to predict its behavior. To allow for the complex coastal morphology, a sea distance function has been implemented within the data weighting, which is shown to significantly enhance the solution. The tuning process allows for independent validation giving a standard error of the resulting surface of 0.2 m for areas with no tidal observations.

## 1. Introduction

The United Kingdom Hydrographic Office (UKHO) commissioned the Vertical Offshore Reference Frames (VORF) project in 2005, with a final completion date of early 2008. The project was led by University College London (UCL), with contributions from the Danish National Space Centre (DNSC) and the Proudman Oceanographic Laboratory (POL). The objective was the research and development of a set of marine reference surface models for U.K.–Irish waters, which would permit the transformation of coastal and hydrographic data between land, sea, and global vertical datums. Of particular import was the modeling of Chart Datum in ETRF89, which would permit the integration of GPS height observations into the hydrographic data reduction process and lessen the reliance on coastal and offshore tide gauges and tidal models (Ziebart et al. 2007; Iliffe et al. 2006). Chart Datum is an arbitrary level used as the reference plane for both tidal predictions and the depths in marine charts. In the U.K. and Ireland Chart Datum is determined to be identical or very close to lowest astronomical tide (LAT). LAT is defined as the level below which water is not expected to fall given normal meteorological conditions, and may be derived by the analysis of a number of years of tidal data/predictions (normally 18.6 yr to account for the full nodal cycle).

In our method an important step in deriving a model of Chart Datum is to first determine a model of lowest astronomical tide with respect to the ETRF89 ellipsoid, and this in turn is separated into two successive steps of deriving a model of mean sea level (MSL) with respect to ETRF89 and then deriving the height below this of LAT through tidal modeling (Iliffe et al. 2007a). This paper is concerned with the derivation of a model of LAT with respect to MSL. A previous paper was concerned with the derivation of a mean sea level surface with respect to ETRF89 (Iliffe et al. 2007b). Further tidal surfaces mean low water springs (MLWS), mean high water springs (MHWS), and highest astronomical tide (HAT) are similarly required; however, this paper will concentrate on the more significant LAT. The area of interest is bounded by 46.80°–63.90°N, 25.00°W–3.50°E. For the entire area the required resolution is 0.008°.

The available data sources are discrete measurements of tidal levels at tide gauges both on the coast and offshore, and continuous predictions of tidal levels from hydrodynamic tide models and satellite altimetry–derived global ocean tide models. The desired solution is to produce an LAT surface from an optimally merged combination of all data sources, with the coastal tide gauges being considered as “truth” data. Away from the coast the tide gauges are envisaged as having little or no effect, while in the coastal zone it is desired that the resulting LAT surface sufficiently represents the directly observed regimes. A further prerequisite is a smooth surface with no discernable transitions between zones.

Previous studies have highlighted the difficulties inherent to predicting tidal levels in areas that have had no in situ measurements (Le Provost 2001). Hydrodynamic models such as those used in operational surge forecasting (e.g., Flather 2000) have the advantage of spatial continuity; however, they are sensitive to input bathymetry and are not always sufficiently well resolved to recreate the higher tidal harmonics responsible for tidal curves in very shallow water. For these reasons they should not be the sole source of datum modeling in the coastal zone. Hess (2002) used a numerical solution of Laplace’s equation to interpolate tidal constituents and posited the use of his work for the determination of ellipsoid-referenced water levels; while Martin and Broadbent (2004) interpolated tidal data using a simple distance-based method. In addition to these, the least squares collocation method in Iliffe et al. (2007b) has been used as an interpolation method on irregularly spaced sea level data in the coastal zone.

To perform the merging and interpolation of the data sources a modified implementation of the thin plate spline (TPS) method is introduced, which has the capability to perform interpolation while accounting for complex coastal morphology. The TPS method is a 2D interpolation process whereby the behavior of a thin metal sheet is simulated [method introduced by Duchon (1977)]. This theoretical metal sheet may be constrained to pass through a number of control points (*x _{i}*,

*y*,

_{i}*z*) and, subsequently, the bending energy required for passing through these points is minimized, thus resulting in a smooth, efficient surface. The method is used for a number of purposes, for example, the construction of digital elevation models (Mitas and Mitasova 1999), coordinate transforms in computer graphics (Donato and Belongie 2002), and extraction of tidal flow patterns (Vennell and Beatson 2006).

_{i}Coastal morphology is modeled by the use of a function that models “sea distance,” which is the shortest distance that it would be necessary to travel between two points by sea. This influences the weighting given to individual data points. To tune the interpolation procedure so that the sea distance mechanism has optimal influence, a test is introduced whereby known tide gauge data points are removed, in turn, from the solution space and their behavior is subsequently predicted by the use of the interpolation procedure and the surrounding information. Thus, the set of weighting parameters that best predict the known data points may be used for the full surface computation.

The TPS method models a real, constrained, physical system and therefore is likely to be of higher applicability to the modeling of tidal surfaces compared to least squares collocation, which is more suited to the prediction of data in an equipotential system. Both methods are implemented for this paper in order to afford comparisons.

The United Kingdom and Ireland have a very tide gauge–rich environment (UKHO 2009; POL 2008), having over the years many tidal observations with excellent spatial distribution. Consequently, in a derivation of a tidal surface such as LAT, optimal use of this data is desired. Further data sources are satellite altimetry–derived tidal models and hydrodynamic tidal models. A number of these wide-area tidal models are utilized and assessed by comparisons with in situ data to determine suitability for offshore areas.

## 2. Data sources

Data for this study may be generally said to fall into one of three definitions: point data from tide gauge observations, wide-area hydrodynamic tide models, and global satellite altimetry–derived tide models. The *onshore* tide gauge data points are sourced from the Admiralty Tide Tables (ATT) of which there are approximately 700 in the area of interest (UKHO 2009), in addition data from 180 *offshore* gauges supplied by the UKHO were used. The duration and epochs of observation at individual gauges vary considerably; for example, the ATT tide gauge at Plymouth, England, has 19 yr of observations, while the Lizard Point tide gauge has 1 month. This impacts on the availability of tidal level observations; Plymouth has sufficient data to enable a direct prediction of LAT while Lizard Point does not. However, what is available at all gauges is the level known as Chart Datum (CD), which has approximated LAT by scaling the available lowest level from the shorter time series through a comparison with identical time series at nearby primary tide stations where LAT had been derived from long-term observations. In practice we have found that where CD and LAT are both known there is an RMS difference of 0.12 m, and therefore to increase the number of observations available on the coast we have used CD as a proxy for a direct observation of LAT, with the proviso that the appropriate weight to be used for these points should be tested. This has been done and is discussed in section 5. Consequently, values of LAT at tide gauge locations are derived from

where is the height of LAT with respect to mean sea level (thus a negative quantity), *Z*_{0} is the observation of mean sea level above local CD, and , where available, is the height of lowest astronomical tide above local CD. In total there are 391 gauges tide gauges with LAT observations from this method. A further 97 ATT gauges had tidal levels derived using

where CD is chart datum expressed in ETRF89 and MSL is mean sea level in ETRF89 from the previously computed mean sea level model. Note this is an approximation to LAT only used for gauges where there exists a chart datum connection to a recognized datum but no observations of LAT or in situ *Z*_{0} are available. Figure 1 plots the locations of all gauges used in the study.

The wide-area models are as follows:

North and Irish Sea and English Channel (NISE10) hydrodynamic tide-surge model (Flather and Williams 2004)

Center for Space Research version 4.0 (CSR4.0) global ocean tide model (Eanes and Bettadpur 1995)

Global Tide Finite Element Solution (FES2004) global ocean tide model (Lefèvre et al. 2002)

TOPEX/Poseidon crossover solution (TPXO 7.0) global ocean tide model (Egbert and Erofeeva 2002)

Goddard Ocean Tide (GOT00.2) global ocean tide model (Ray 1999)

The NISE10 model is a two-dimensional, depth-averaged numerical model with a horizontal resolution of 3.5 km. It is part of the U.K. operational surge forecasting suite, where it is forced by surface pressure and 10-m wind fields from the mesoscale atmospheric model. In this project it is used with tidal input only, specifically the largest 26 constituents applied to the lateral boundaries through a radiation condition. The remaining models form a subset of a family of global ocean tide models derived from data from various satellite altimetry missions.

To acquire the full tidal range within the 18.6-yr nodal cycle the models were run over extended periods (usually 20 yr) across the desired model domain; the four properties of interest derived from each run were modeled lowest tide, modeled highest tide, mean high water springs, and mean low water springs. To assess the wide-area models comparisons were made with the tide gauge observations, using LAT as the chosen metric. Tables 1 and 2 summarize the comparisons between all wide-area models and the onshore and offshore gauges, respectively. These have been calculated from the differences between the values of LAT at the tide gauge (TG) and the value of LAT derived from the wide-area model (WAM) at the tide gauge location, as follows:

A further assessment is performed by deriving the RMS agreement between each possible pair of the wide-area models. This is not an independent validation of model accuracy; however, it is a necessary precondition that two models must agree with one another if both are to be said to accurately represent tidal behavior. The RMS differences between the models (at the offshore gauge locations) are given in Table 3. As may be expected, because of their similar reliance on Ocean Topography Experiment (TOPEX)/Poseidon data, the models derived from satellite altimetry have the smaller differences.

Of particular regard is the expected degradation of the altimetric-derived models close to land and that it is this area in which the hydrodynamic modeling performs better. Offshore, however, the altimetric-derived models better represent the directly observed tidal behavior. Consequently, the LAT surface derived from the satellite altimetry ocean tide model CSR4.0 was selected to be the single data source for all areas at least 30 km from the coast because of its good agreement with the offshore tide gauges, good agreement with the FES2004 model, and superior computational efficiency. As an additional benefit it may be said to be a “safer” option for hydrographic application because of it predicting, on average, a greater tidal range. Within 30 km of the coast the CSR4.0 data was suitably decimated and points rejected if either within ∼14 km of land (Deng et al. 2002) or seen to insufficiently reflect the local tidal regime in comparison with in situ gauges (such as in the Bristol Channel).

## 3. Review of the TPS method

The TPS function has the following form:

where *z* is the value interpolated at point (*x*, *y*) using a dataset of *p* control points

Here [*a*_{1}, *a _{x}*,

*a*]

_{y}^{T}and

*w*are, respectively, coefficients and weights that have been derived from

_{i}and 𝗞 is populated from the *p* control points as follows:

The *i*th row of 𝗣 is (1, *x _{i}*,

*y*), 𝗢 is a 3 × 3 matrix of zeros,

_{i}**o**is a 3 × 1 column vector of zeros, and

**is a**

*υ**p*× 1 column vector consisting of the surface values

*z*. Lambda (

_{i}*λ*) is termed the

_{i}*regularization*parameter and controls the degree of exactness of the interpolation. If all

*λ*s are zero, then the interpolation is exact; that is, the spline surface passes exactly through the control points, and as all

_{i}*λ*→ ∞, then the resulting interpolation is equivalent to a least squares fit plane. By analogy,

_{i}*λ*defines the rigidity of the plate, and the cases (

_{i}*λ*= ∞) and (

_{i}*λ*= 0) run from the extreme weakness to the extreme strength of the control points.

_{i}## 4. Application of the TPS to an LAT surface computation

The various data sources have differing areas of relevance and applicability. Away from the coast the altimetry-derived modeling is of utmost importance; consequently, outside a zone defined by a buffer ∼30 km from the coast the LAT surface points are derived purely from the bicubic interpolation of the CSR LAT surface. In the coastal zone TPS functions are used to interpolate to the surface point locations. The quantity of source data is such that it is impractical to use all for each interpolation point; therefore, an approach is used whereby a specific TPS function is created for each point. The control points are selected as follows:

The nearest eight valid LAT values derived from the CSR4.0 model. Even if the CSR4.0 data are distant, the inclusion of these points ensures a smooth transition from the coastal to the offshore zone.

LAT surface value at the point in question derived from the hydrodynamic model NISE10.

All tide gauge data within a radius of 80 km of the point in question.

To ensure a smoothly continuous transition from the inshore to the offshore zone all non-CSR data were removed from the analysis if situated at least 14 km from the coast. To enable true distance calculations for the radial basis function in Eq. (4) all station and point coordinates are converted to Eastings/Northings (kilometers) in OSGB36. This moving window approach for the control points, using a subset of all data and individual TPS functions, is necessary not only for computational reasons but also to account for coastal morphology. Each estimation point will have a unique relationship with the control points due to location relative to the coast. The coastal morphology is modeled by the introduction of the concept of sea distance (Iliffe et al. 2007b). This represents the shortest distance that it would be necessary to travel between two points by sea (which is different to the distance along the coast as the path traveled will “jump” from one headland to another). For each individual control point the associated value of the regularization parameter *λ _{i}* in each TPS function is derived from

where *r*_{sea} is the sea distance between the estimated point and the control point and *r*_{straight} is the straight distance between the points. The *λ*_{0} is the base regularization for control point *i* and *β* is a scale factor determined to optimize the deweighting of points as a result of sea distance. Thus, if there is no intervening land between the points, then there is zero adjustment to the original weighting; however, the farther around a headland a point is located, the less influence it will have on the final solution. The scale factor *β* determines the rate at which points around a headland still have influence. As *β* approaches zero the *i*th data point will have little or no influence on the interpolation, and as *β* approaches infinity the presence of intervening coastal morphology will have no impact on the solution. If there is no intervening land there is no morphology penalty, *r*_{sea} is identical to *r*_{straight} and the sea distance term is eliminated. The values for *λ*_{0} and *β* depend on the tidal behavior in the area being modeled and may be estimated by a tuning process.

True sea distance may be computed by the implementation of Dijkstra’s shortest path algorithm (Dijkstra 1959) with a set of polygons modeling the coastline (Fig. 2a gives an example true sea distance path between two points). The recursive nature of the algorithm combined with the high resolution of the coastline polygons resulted in a computational overhead too high to be of practical use. A second, approximating, sea distance algorithm used a method whereby the polygons intersecting the straight path between the start and end points are split into three equal length sections. The section bounds are determined by the projection of the polygon points onto the straight path as can be seen in Fig. 2b. The polygon vertex within each section at greatest perpendicular distance from the line of intersection is stored. Two candidate distances would be derived for each polygon by the summation of the distances between the stored points for both the clockwise and anticlockwise routes, with the shortest distance then deemed as the return value.

Systematic validation of the approximating method shows that, in comparison to true sea distance, it overestimates by an average of 20%. This is a consistent bias; therefore, the subsequent optimization of the sea distance scale factor *β* will compensate in the TPS function.

## 5. Tuning of the TPS function

The regularization parameter, *λ _{i}*, controls the degree of smoothing in the TPS approximation. If all

*λ*s are zero, then the interpolation is exact; however, this would make no account for issues such as data noise, and would produce a somewhat unrealistic solution. Furthermore, this would take no account of any coastal morphology and the weighting between the various data sources. However, if the

_{i}*λ*s are too great the solution approximates a best-fit plane. Therefore, an empirical method may be implemented in order to determine the optimum values of the base regularization

_{i}*λ*

_{0}for the control points. Similarly, the effect of coastal morphology, as modeled using sea distance, may be optimized by the estimation of

*β*from Eq. (6).

To tune the TPS interpolation, and therefore derive optimal values of *β* and *λ*_{0}, a “take one out” test was implemented. For this, ATT tide gauge data points were removed in turn from the solution and the remaining data points used, along with the TPS methodology, to predict tidal levels at this deleted point. A range of values for the *β* and *λ*_{0} parameters were tested and comparison statistics derived to determine the degree of overall agreement with the independent truth data. In addition further tests were performed with a range of values for *λ*_{nise10}. This is the regularization parameter for data from the NISE10 hydrodynamic model that is designed to have dominance given the absence of availability of tide gauge and CSR data sources.

The RMS agreement between the observed and predicted values is plotted in Fig. 3 for ranges of values of both *β* and *λ*_{0}, with *λ*_{nise10} constant at its optimum value of 10^{3}. The combination of parameters giving the best agreement (RMS 20.20 cm) with the in situ observations is seen to be 1.00 km and 10^{2} for *β* and *λ*_{0}, respectively. Evident from the plot is the parabolic nature of the sea distance weighting; too little or too much weighting results in poorer results.

Included within the overall RMS statistics are individual values up to 1.04 m, which occur at locations of sparsely available data in rapidly changing tidal regimes (e.g., predicting tidal levels at the head of a loch from data available at the open coast).

A further test repeated this work, but eliminated the incorporation of the NISE10 hydrodynamic modeling; the comparable agreement, under identical parameterization, between the predicted and computed tide gauge values was 20.94 cm.

The base regularization parameter *λ*_{0} was estimated as a single value for use on all tide gauge control points. However, as stated in the section on data sources, of the onshore gauges used in this study there are a number where LAT values are inferred using local values of CD. Therefore, to assess the reliance on tide gauge data with relatively small observation spans, tidal levels may be predicted at the primary ATT gauges (i.e., those that have sufficient data for LAT to be determined) under a number of different weighting scenarios. It was found that deweighting data from these secondary gauges by modifying their value of *λ*_{0} had no significant effect on the solution.

For comparison, the collocation procedure similar to that described in Iliffe et al. (2007b) was also implemented, using a variable covariance model that accounted for regionally varying tidal regimes, and a similar tuning of the data weights performed. The best-case comparison achieved between the predicted and observed values was 26.22 cm.

## 6. Results

A LAT surface was derived, using the set of parameterizations thus derived, for all points in the VORF domain combining the tide gauge data, hydrodynamic modeling, and altimetric-derived modeling. Because of the computationally intensive nature of the approximating sea distance function, the UCL parallel computing facility CONDOR was used to partition the job into manageable segments. The resulting LAT surface is seen in Fig. 4. Final statistics show that for the onshore ATT tide gauges the RMS agreement is 5.81 cm. However, this is purely a measure of precision as all possible data have been used for the surface computation.

Figures 5a–c illustrate the behaviors of the interpolated LAT tidal surface under differing parameterizations of the TPS regularization in an area proximate to Portland Bill off the south coast of England. This is a particularly interesting area as it has a rapidly varying tidal regime (being close to an amphidromic point) over a relatively short distance separated by a narrow spit of land. The three plots depict the respective results when (i) sea distance has no input in the interpolation process, (ii) sea distance is used to exclude all points when even slightly around a headland, and (iii) sea distance has optimal influence.

## 7. Conclusions

A LAT surface (with respect to mean sea level) for the U.K.–Irish waters, encompassing both the coastal and offshore zones, has been derived using an interpolated combination of tide gauge data, satellite altimetry, and numerical tidal modeling. This interpolation was effected by the use of a thin plate spline model. Generally, where altimetry or tide gauge data were available, the thin plate spline interpolation algorithm has given this priority. Toward areas of sparser data the levels derived from the NISE10 hydrodynamic model have gradually been given precedence. The dominant features are the extreme tidal minima in the Bristol Channel and Baie du Mont St Michel, resulting from the large tidal ranges there. The plot in Fig. 4 also reflects the well-known amphidromic system of the European continental shelf (see, e.g., Pugh 1987).

The thin plate spline method has been shown to be more suited to the task of tidal level interpolation compared to the least squares collocation method; the improvement was equivalent to a 23% reduction in modeling error. By virtue of its mathematical basis, the collocation approach is more readily adapted to the interpolation of equipotential fields.

An empirical procedure has been undertaken to assess the extent to which tidal observations can be expected to extend their influence around headlands, the images in Fig. 5 clearly emphasizing the need for the optimum incorporation of coastal morphology when performing near-shore oceanographic interpolations. This method may also be applied to the interpolation of further oceanographic phenomena as well as tidal predictions for specific epochs. Scope for further work includes examining regional variations in the effect of tides around headlands by performing localized tuning and the application of the TPS method and tuning to other regions.

Independent analysis of the surface would involve the introduction of new observations or—equivalently—the exclusion of data points from the solution and then comparison of these with the derived values. However, the process of successive elimination that has been performed as part of the tuning process is in fact the direct equivalent of this and has demonstrated a 20-cm uncertainty in areas in which there is no in situ tidal data. Additionally, contractors working on the U.K.’s Civil Hydrography Programme have been given access to the VORF surfaces since late 2008 and may use the model provided that spot checks are made against new determinations of Chart Datum from cotidal techniques or offshore gauges. For the survey to be accepted, agreement is required at the VORF 2*σ* level of 20 cm, and this has not been exceeded at the locations so far surveyed (offshore Orkney Islands, Northern Ireland, Weymouth, and the sandwave fields off East Anglia).

## Acknowledgments

The authors acknowledge the contribution to this study by the United Kingdom Hydrographic Office, through funding and provision of data. Parallel computing facilities were provided by the UCL Computing Resource Allocation Group (CRAG). Further tide gauge data were provided by the National Tidal and Sea Level Facility (NTSLF) and Permanent Service for Mean Sea Level (PSMSL), both hosted by POL.

## REFERENCES

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

**,**

## Footnotes

*Corresponding author address:* J. C. Iliffe, CEGE, University College London, Chadwick Building, Gower Street, London WC1E 6BT, United Kingdom. Email: jiliffe@cege.ucl.ac.uk