Spatial Interpolation of Tidal Data Using a Multiple-Order Harmonic Equation for Unstructured Grids

A general spatial interpolation method for tidal properties has been developed by solving a partial differential equation with a combination of different orders of harmonic operators using a mixed finite element method. Numerically, the equation is solved implicitly without iteration on an unstructured triangular mesh grid. The paper demonstrates the performance of the method for tidal property fields with different characteristics, boundary complexity, number of input data points, and data point distribution. The method has been successfully applied under several different tidal environments, including an idealized distribution in a square basin, coamplitude and cophase lines in the Taylor semi-infiite rotating channel, and tide coamplitude and cophase lines in the Bohai Sea and Chesapeake Bay. Compared to Laplace’s equation that NOAA/NOS currently uses for interpolation in hydrographic and oceanographic applications, the multiple-order harmonic equation method eliminates the problem of singularities at data points, and produces interpolation results with better accuracy and precision.


Introduction
Coastal ocean water level observations are only available at a limited number of locations, usually at tide stations located primarily along the coastline.Therefore, any tidal characteristics derived from data analysis of water level measurements, like tidal constituents, residual water levels, and tidal datums, are limited to these tide station locations.However, the spatial distribution of tidal properties is always needed in the open waters far from coasts for navigation and hydrographic survey purposes.Numerous approaches have been developed to interpolate scattered geo-spatial data [1][2][3], including self-correcting [4] and thin plate spline methods [5].These methods have been developed for rectangular, simply-connected regions and thus would have difficulty in interpolating tidal properties in irregular regions with complex waterways interspersed with land features such as islands and peninsulas.To overcome this difficulty, NOAA's National Ocean Service (NOS) has used a discrete tidal zoning process in which the survey area is divided into a number of polygon-shaped zones [6].Each zone is assigned to a range ratio and time difference; these values are applied to tidal data at the appropriate reference tide station to obtain the range and phase in the zone.Since the value within each zone is a constant, there is a discontinuity in the interpolated tidal characteristics across the boundaries between adjacent zoning polygons.The development of zoning is also subjective and requires extensive cross checking and quality control.
NOS more recently adapted a solution of Laplace's equation for spatial interpolation of tidal constituents, residual water levels, and tidal datums in irregular regions [7][8][9].Laplace's equation is   where Δ is the Laplace or harmonic operator x y     , f is the property to be interpolated, i F is the forcing value at point  , and N is the number of data points.The equation is solved iteratively using finite difference methods in a regular rectangular mesh grid [9].Laplace's equation has been proven to be conservative in that the minimum and maximum of an interpolated field are always at an input point or at the boundary.If the normal derivative at the boundary is zero, the range of the interpolated field is within the minimum and maximum of the data values.This property will eliminate the risk of overshooting (i.e., producing interpolated values that are higher/lower than the maximum/minimum data values) of the interpolation/ extrapolation, although it should be remembered that the maximum value of certain properties, such as tide range, is undetermined, since the spatial distribution of the range is usually not known unless numerous tide stations are installed.But Laplace's equation has some limitations: it may be too conservative, and most of all, singularities (i.e., discontinuities in the slope) exist at all data points.Another option is to use a higher-order harmonic differential equation.The pioneer in this approach was Briggs [10], who used the biharmonic equation to interpolate gravity and aeromagnetic fields but was limited to regular grids and rectangular, simply-connected regions.The equation is derived by minimizing total domain-integrated squared curvature, , where   is the domain on which f is defined.Smith and Wessel [11] further extend Briggs' work by applying a linear combination of a biharmonic operator and Laplace's operator: where T is a tension parameter.Their applications were also restricted to regular grids and rectangular, simplyconnected regions.Equation ( 3), which we call minimum curvature with tension (MCT), can be thought of as the analog of a thin elastic plate, bent to fit the data points with a tension applied at the boundary.Tension helps to reduce the problem of overshooting.Higher order harmonic interpolation functions such as and are also widely applied in graphic processing [12], but their use rarely has been reported in hydrographic and oceanographic applications.We include them here for completeness.
In this paper, we present a new, generalized interpola-tion method by solving a partial differential equation (PDE) with a linear combination of different orders of harmonic operators (the 1 st , 2 nd , 3 rd and 4 th ) on an unstructured, triangular-element grid.Our approach has three main advantages over previous solution methods.First, the use of multiple orders of differentials allows for more possible harmonic functions to be incorporated into the solution.With proper parameter selection, we can eliminate the apparent flaw in the occurrence of singularities at the data points when using Laplace's equation alone.Second, the use of unstructured grids to represent coastal water avoids issues related to interpolation over intervening land features, the necessity of over-water distance calculations, and the problem of multiple connected regions.The third is the streamlining of the multiple-order harmonic PDE interpolation into one simple implementation and solving the equation implicitly without the need for iterations.We note that this approach incorporates the (regular grid) equations of Hess [7], Briggs [10], Smith and Wessel [11], and Botsch and Kobbelt [12] through parameter selection.
We explore the application of the multiple-order harmonic equation to oceanographic properties, namely tidal fields.The outline of the paper is as follows: first, we introduce the mathematics of the governing equations, then we give accounts of the numerical method for the solution, the quantitative measure of error, and the jackknife method of testing.After that, we present four cases of applications: an idealized function in a square domain; Kelvin wave reflection in a semi-infinite rotating channel; M 2 , S 2 , K 1 and O 1 tidal constituents in the Bohai Sea; and the M 2 tidal constituent amplitude and phase in Chesapeake Bay.Finally, we offer conclusions based on the applications of the generalized multiple-order harmonic equation interpolation method.

Governing Equation and Boundary Conditions
The generalized governing equation for spatial interpolation, which incorporates multiple orders of harmonics, is where k order harmonic operator, and with the constraint We note that, with appropriate selection of coefficients β, this equation includes Laplace's equation (Equation (1)), minimum curvature (Equation ( 2)), and MCT (Equation (3)).Boundary conditions are specified as: is an adjustable dimensionless parameter and f n   is the spatial average of the normal derivative in a small region adjacent to the boundary.Equation (5a) allows some flexibility in simulating tidal properties, and has been used successfully in computing tidal constituent distributions in Galveston Bay and San Francisco Bay [8].
The left hand side of Equation ( 4) is a linear combination of different order harmonic operators.By using different combinations of k 4) provides a class of interpolation methods (Table 1).In our applications, we avoid any order of harmonic operator higher than 4 due to the high computational cost, high computer memory requirements, and the natural tendency towards overshooting or oscillation, which can be especially pronounced in a domain interior in the case of interpolating from boundary data without constraints from internal data points.
Equation (5a) can have different physical explanations, depending on the order and the value of α.For Laplace's equation, in analogy to a stretched membrane, 0   corresponds to a zero slope at the boundary, and 1

 
to a constant, non-zero slope at the boundary.For 2 nd order or MCT, the equation describes an elastic thin plate.In this case, 0   1 forces the solution to flatten at the edge and   corresponds to a free-edge condition without bending stress.In general, the α value indicates the degree of balance between the 1 st and 2 nd normal derivative at the boundary, and it provides a tool to adjust the boundary condition and internal field.In our current application, we take 0.9

 
, which has been shown to yield realistic solutions to tidal constituent distributions [7].In addition, many of the equations represented by the different combinations listed in Table 1 can be derived directly from optimization of an energy function through variational calculus [12,13].For example, for k = 1, Laplace's equation describes a surface which minimizes the area.For k = 2, a surface that minimizes the curvature is represented, and k = 3 represents a surface that minimizes the variation of linearized curvature.In the onedimensional case, the 1 st , 2 nd , 3 rd , and 4 th order harmonic interpolation equations represent piecewise linear, cubic, quintic, and septic polynomial interpolation, respectively.
As will be shown in the next section, the PDE (Equations ( 4), ( 5)) can be readily solved using mixed finite element methods with linear triangular elements.

Mixed Finite Element Method
Laplace's and Poisson's equations can be readily solved using a standard finite element method with linear triangular elements [14].To directly solve a higher order harmonic equation, we employ a mixed finite element method [15][16][17] to transform Equation ( 4) into a 1 st order linear Poisson's equation system.First, we assume We then transform the high order linear harmonic Equation (4) into a low order linear PDE system, with the boundary conditions, Equation ( 6b) is a screened Poisson equation.Applying the mixed finite element method to the above mixed PDE system, i.e., Equation ( 6) with boundary conditions Equation ( 7), and using a discrete Laplace operator [18,19], we have a linear system the unstructured triangular mesh nodes.This sparse linear equation AG = B is solvable by straightforward array arithmetic (we used the MATLAB© built-in linear system function).The equation not only solves the interpolation field, but also produces the k th (k = 1, 2, 3,•••, K − 1) order Laplacian field as a byproduct.The method can also be used for solving the Poisson equation, the Helmholtz equation, and the L. SHI ET AL. 1428 bi-Helmholtz equation with the addition of appropriate terms to Equation (4).

Experiments, Error Estimation and the Jackknife Method
To test the generalized multiple-order harmonic spatial interpolation method, we use two different types of data sets.The first type of data set has a known reference field such as an array of observations or an analytic solution.
The second type of data set does not have a known reference field, as only the values at a limited number of locations are known.Each type of data set requires a unique experimental approach.

Experiments with a Known Reference Field
Assuming there is a two-dimensional distributed reference data set A containing N points, we randomly subsample A to create a data set B, a subset of A with M points, where M N  .The spatial interpolation method uses B to create an interpolated field C, which has the same size as A and contains values at all locations in the domain.The number of points selected, M, starts at a small value (for example, we use 5), then increases gradually by adding 5 randomly selected points to the previously selected points.The points are selected from the entire domain.Every time 5 additional points are added, interpolation is performed to create a new field C, and the associated maximum absolute error (MAXE), mean absolute error (MAE), and root mean square error (RMSE) are calculated by comparing the interpolated field C with the reference field A. This experiment is designed to quantify the performance of different methods under different data density scenarios.
Since a tidal property is known primarily at locations along the coastal boundary of the ocean domain, we also conduct one parallel experiment with a realistic constraint: select subset B from boundary points only, again starting with 5 randomly-selected boundary points, then adding another 5 boundary points, and so on.The error statistics are calculated with this data series in the same way as for the random selection from all points.

Experiments with an Unknown Reference Field
In many coastal applications, there is no known background reference field, and only scattered station data located mostly along the land-sea boundary are available.Because of a limited number of stations, subsampling may encounter the problem of small total sample size.To reduce the uncertainty due to sample size, we will fix the subsample size and repeat the subsampling for any possible combination.In our application, the delete-1 jackknife method will be used [20].Specifically, a subsample size of N − 1 out of a total of N samples will be used, and a total of N repeat samplings will be performed by re-moving one station at a time.At the end we have N error values, one at each station, and the same MAXE, MAE, and RMSE will be calculated.This experiment is designed to quantify the performance of different methods using limited data points at fixed locations without a known reference field.

Applications
The performance of the new interpolation method using the high order harmonic PDE is evaluated in four test cases, each with different numbers and distributions of data points and varying complexity of the domain.In each case we first tested pure 1 st , 2 nd , 3 rd and 4 th order interpolation.If 3 rd and 4 th order interpolation is not as good (i.e., has higher error measures) as both 1 st and 2 nd order interpolation, a few MCT cases are then tested, and the 3 rd and 4 th order results are not presented.Usually we start by testing the use of all internal and boundary data points, then proceed to testing the use of boundary points only.For the cases in which an analytic solution is known, we test the solution using various harmonics with progressively larger numbers of data points.

Analytic Function
Here we apply the PDE to an idealized square basin with a known analytic function describing the distribution.Since the reference field is a known analytic function, we may compare the interpolated solution and reference field at any desired point, especially at randomly-selected internal points.The analytic function, Peaks, included in the MATLAB© software, takes the form of: Peaks is a combination of three two-dimensional Gaussian functions, each multiplied by a polynomial function, which ensures that the function approaches zero at locations far from the origin (in our idealized basin, x and y each varies between −3 and 3).The function ranges in value from a minimum of −6 to a maximum of 8.In our experiments, we represent the square basin with an unstructured grid mesh having 161 equally-spaced nodes in the x-and y-direction.The reference field was generated by evaluating Equation (8) at the nodes of the mesh.

 , h x y
A flaw in the Laplace's equation solution is very obvious from a plot of the test with 90 random data points, in which singularities occur at the data points (Figure 1(a)), while the higher order harmonic interpolation generates a very smooth surface (Figure 1(b)).The surface contour plots (Figure 2) indicate that the pure 3 rd order Open Access IJG L. SHI ET AL. harmonic interpolation gives the best results, compared to either the 1 st or 2 nd order interpolation with 90 random data points (this conclusion is also supported by the plots in Figures 3a-c).The MAE for Laplace's interpolation decreases steadily with an increase in the number of data points.At low data point numbers (M < 25), Laplace's interpolation performs better than all of the higher order harmonic interpolations.But at the same time, at low data point numbers, the overall error is relatively high, i.e., on the same order of magnitude as the interpolated values.With an increase in the number of data points, the higher order harmonic interpolations, especially the 3 rd and 4 th order harmonic equations, improve much faster than the Laplace's interpolation.When the number of observed points increases to 100, the 4 th order outperforms the 3 rd order.The trends from all three error measures are similar in that they decrease with increasing number of data points, although the errors decrease faster with the higher order.As for the use of boundary points only, we note that since the values along the boundary are all very close to zero, interpolation using only boundary points will not yield much useful information about the interior of the basin.We therefore do not present any results from boundary only data.
In general, high order interpolation gives a smooth solution with a continuous derivative within the domain, including at the data points.For all interpolation methods, errors decrease when more data points are used.The higher-order PDEs outperform lower-order PDEs, except when the number of data points is small (M < 25).

The Taylor Problem
The original Taylor problem, which was first presented and solved analytically by Taylor [21], is the Kelvin wave reflection in a semi-infinite rotating channel without friction.In our experiments, to be consistent with Taylor [21], we simulated flow in a channel that is 500.4km wide (north to south), 74 m deep, and 1000 km long, east to west (Figure 4).The Coriolis coefficient is 0.000119 s −1 (corresponding to latitude 54.46˚N) and the period of oscillation was set to 12.4 hours to approximate the M 2 period.This example was selected because it has an analytic solution for tidal amplitude and phase distribution, thus providing a reference field, and unlike the previous case it has non-zero boundary values.The dominant feature in the analytic solution is a chain of amphidromes along the center line of the channel.Rienecker and Teubner [22] introduced friction (with a coefficient of 0.00005 s −1 ), and with friction the amphidromes shift laterally from the center line of the channel to the right side of the channel if facing the inbound wave (Figure 4a).There are two amphidromes within the first 1000 km of the channel.The first amphidrome is close to the center line.Further away from the closed western end, the wave is attenuated and the second amphidrome is closer to the southern boundary.Depending on the bottom friction, the amphidrome may be degenerated, and the center may move inland, producing a virtual amphidrome.
To create the interpolation we use an unstructured triangular grid having a uniform 5 km resolution with a   .This complex number interpolation is also used for the interpolation of the M 2 tidal constituent in the following Bohai Sea and Chesapeake Bay cases.
Two tests were conducted for the Taylor case, one interpolated from random internal and boundary data points and the other from only random boundary points.For both tests, along with the usual statistics, another statistic was computed: the average distance between amphidromes in the reference field (here the analytic solution) and corresponding amphidromes in the interpolated field.
In the internal points case, locations of the amphidromes are able to converge quickly (M > 50) to the analytic positions (Figure 5d).For both amplitude and phase, the 2 nd order interpolation dramatically reduces the error in MAXE, MAE, RMSE, and the distance between the amphridromes (Figure 5) in comparison with the 1 st order interpolation.The 1 st order interpolation error is about 10 times that of the 2 nd order interpolation error in all three measures.The 3 rd order interpolation is consistently better than the 2 nd order, but only by a small mar-gin (i.e., only a few percent).The 4 th order is comparable to the 3 rd order, or slightly better, especially for large data numbers (M > 600).However, the 4 th order interpolation performance can deteriorate (most likely from overshooting) as demonstrated by a sudden increase in RMSE at M = 250 (Figure 5f), although it is still better than the 1 st order and performance improves with further increase of data points (M > 600).
In the boundary points only case, besides the pure 1 st and 2 nd order interpolations, we tested three interpolations using MCT, with pairs of coefficients   1 2 ,   = (0.5, 1), (0.2, 1), and (0.1, 1).We do not present results for the 3 rd and 4 th order for the boundary case because we found that the 2 nd order interpolation consistently outperforms both the 3 rd and 4 th order in all error measures.The errors decrease with increasing number of boundary data points for both amplitude and phase for the 1 st order and the three MCT cases.The order of their performance in decreasing magnitude of error corresponds to   values of (1, 0), (0.5, 1), (0.2, 1), and (0.1, 1).By contrast, the pure 2 nd order performance is different for amplitude and phase.For amplitude, there is relatively high accuracy at a low data point number (M < 35), but performance deteriorates very quickly with the increase of boundary point number, and the 2 nd order interpolation is outperformed by all other methods (Figures 6a-c).For phase, the 2 nd order performs the best (Figures 6e and f).With the increase of boundary data points, the phase and amphidrome position improve quickly initially, and after a certain number, the improvement from the additional points levels off.That is reflected in the flat line from M = 150 points and beyond (Figures 6d-f).The interpolated amphidrome is further offshore as compared with the analytic results (Figure 6b), especially for the rightmost amphidrome.Unlike the internal point case, information from the boundary, no matter how detailed, cannot indefinitely improve the results or provide all necessary information for the interior, even when there is still a room to improve.Thus we can conclude that for the internal points case  in a semi-infinite rotating channel, the 2 nd order interpolation provides a solution that dramatically reduces the error to 10% of that for the 1 st order interpolation of both M 2 amplitude and phase.The use of the 3 rd and 4 th order terms subsequently further improves the solution, but with a much smaller margin.For the boundary points only case, the phase is best simulated using the 2 nd order or 3 rd order solution (3 rd order results are not presented).For amplitude, though, MCT gives better solutions.Overall, of the limited pairs of coefficients tested, the solution for  1 2 ,    = (0.1, 1.0) provided the best results for both amplitude and phase.Finally, in the absence of internal data, an increase in the number of boundary data points does not improve the accuracy of the position of the amphidrome after a certain number of data points (M > 150).

Tides in the Bohai Sea
The Bohai Sea is a semi-enclosed embayment on the coast of China in the northwest Pacific Ocean, and its tidal system has been well documented from observations and hydrodynamic model simulations [23][24][25].The test data set for this case consists of M 2 , K 1 , S 2 , and O 1 tidal constituents from 31 observed tidal stations in the Bohai Sea [23], of which 29 are near the coast, and two are from islands in the mouth of the Bohai Sea.To carry out the interpolations, we generated an unstructured grid with 101,217 nodes, having on average a 1 km resolution.
We tested five pairs of   ,   to estimate the performance of different interpolation schemes of the complex solution using the jackknife method (Table 2).For all tidal constituents and for both amplitude and phase, except for the K 1 phase, there are more than three combinations of   ,   that perform better than the 1 st order.The 2 nd order is better than the 1 st order in all amplitudes of the four tidal constituents, and better than the 1 st order in MAE, and RMSE for the phase of all tidal constituents.The 2 nd order interpolation has worse phase MAXE for the K 1 , S 2 , and O 1 constituents.For the K 1 phase, the 1 st order MAXE is the smallest of all methods.That is the only case that there is no better alternative  1 2 ,    than the 1 st order in all three measures.The 2 nd order gives the best results in K 1 phase, as measured by the MAE and RMSE (Table 2).
The dominant feature of the results using 2 nd order harmonic interpolation is the amphidromic system (Figure 7).For example, the M 2 tide system indicates an amphidrome in the north, and possibly a degenerated amphidrome in the west.The location of the amphidromes for all constituents is very close to the position estimated by Fang et al. [24] from altimetry and observed station data (Table 3).As in the Taylor boundary case, the interpolated position of the semi-diurnal amphidrome is slightly more offshore than the estimation, while the diurnal amphidrome's position is almost identical.

Tides in the Chesapeake Bay
The Chesapeake Bay, on the U.S. East Coast, is representative of an embayment with a long, complex coast whose shoreline length is extremely large compared with its area.In addition, most data locations are situated along the shore and are thus not representative of internal conditions.In this case, the interpolation of the M 2 tide,  from which the coamplitude and cophase lines are extracted, will be performed and compared using only the 1 st order, 2 nd order, and MCT interpolation.The unstructured mesh grid we used has 318,860 nodes, and a spatial resolution varying between 0.02 and 20 km [26].
For comparison purposes, the coamplitude (Figure 8a) and cophase (Figure 9a) lines for the Chesapeake Bay are hand-drawn contour plots from Browne and Fisher [27], using data from 121 present and historical tide stations.There are two features in the tidal amplitude field that are the most apparent.First, as is the situation in the Taylor case, the amplitude is usually higher in the Bay's east side, or the left side when facing the incoming tide wave.This is due to the existence of a virtual amphidrome (i.e., an amphidrome whose center is located on land) just south of latitude 38˚N [27].Second, the tidal amplitude can either increase or decrease in the direction of wave propagation, while the tidal phase is al- ways monotonically increasing.The behavior of the phase indicates that the tide in the Bay is a progressive wave.
In our Chesapeake Bay application, the M 2 tidal constants are obtained from 50 NOS tide stations that had long term tide observations and tidal constituent analysis.The error analysis using jackknifing (Table 4) indicates the results would be the best when using different pairs of β values for amplitude and phase.If we look individually,  1 2 ,    = (0.01, 1) is the best for M 2 amplitude.
However for M 2 phase, the 2 nd order interpolation is the best.
The resulting interpolated fields of M 2 coamplitude and cophase for the optimal pair of  1 2 ,    are presented in Figures 8 and 9, respectively.The major features of the amplitude (i.e., the virtual amphidrome) and phase are reproduced.The MAE of the amplitude interpolation is about 5.90 cm.The phase distribution is more regular and monotonic with a MAE of 13.45 degree using the 2 nd order interpolation (Figure 9).

Summary and Conclusions
In this paper, we developed a PDE containing multiple, high order harmonics for spatial interpolation of tidal properties in the coastal ocean.The equation is solved on an unstructured triangular grid with a mixed finite element method.The use of an unstructured grid allows for the representation of geometrically complex regions, including islands, and eliminates the need for over-water distance calculations or interpolating across land.In its numerical implementation, the method avoids iteration, which may not converge, and instead uses a simple implicit procedure.Four test cases have been examined, including an idealized function in a square basin, tidal properties in a semi-infinite rotating channel, tides in the Bohai Sea, and tides in the Chesapeake Bay.The results demonstrate that the multiple-order harmonic interpolation method eliminates the singularities at the observed data points that occur in Laplace's equation interpolation (Figure 1b).It also reduces the computed error and improves the accuracy and precision of the interpolated field when appropriate values of the adjustable parameter (i.e., the values of β, Table 1).Computationally, the method is comparable in speed to the solution of Laplace's equation by iteration.In all our test cases, the optimal combination of coefficients k  (k = 1, 2, 3,•••, K) is highly dependent on the spatial characteristics of the interpolated property, complexity of the boundary, and the number and spatial distribution of the data points.In practice, different combinations of coefficients have to be tested to find an optimal combination.It is advisable to test pure 1 st , 2 nd , 3 rd and 4 th order interpolations first.If 3 rd and 4 th order interpolations do not reduce errors, they should be ignored and then MCT cases should be tested.The optimal combination of coefficients can be determined by evaluating the error measures from these tests.
To interpolate the tidal constituents' amplitude A and phase φ, it is good practice to separately interpolate the real part R and imaginary part E of the complex number representation, i Ae R E  i   .In our experience, the interpolated R and E fields are usually smoother than amplitude and phase, with maxima and minima at the domain boundary.
Interpolation from boundary data alone can produce good results for tidal properties (Figures 4, 6e and f), but has its limits.No matter how dense the data points are, the missing information from the interior (i.e., offshore) will never be fully recovered from the boundary data alone.This is demonstrated by the persistent error in the distance between amphidromes between the reference field and the interpolated field (Figure 6d), and the flat line of error measures when increasing the number of boundary points (Figures 6a-f).
In practice, tidal amplitude and phase can be interpolated separately, each using different values of β coefficients to achieve an optimal result, as demonstrated in the Taylor basin and the Chesapeake Bay M 2 tide cases.For amplitude in those cases, MCT is better than either the 1 st order or the 2 nd order interpolation alone (Figures 6a-c; Table 4).For phase, the 2 nd order interpolation is the best in all three error measures in Chesapeake Bay (Table 4), as well as for M 2 phase in the Taylor boundary point case (Figures 6d-6f).The fact that different values of β produce optimal results may be due to differences in the underlying properties of tidal amplitude and tidal phase.In other words, the tidal phase field is L. SHI ET AL. 1436 much smoother than the tidal amplitude field, and therefore the tidal phase is more suitable for a higher order harmonic interpolation.
In summary, our development of the high order harmonic interpolation provides a full set of options for the user to choose from for their specific application.Besides being a better tool to achieve more accurate interpolation results, there are also many potential benefits.For example, in hydrographic and oceanographic applications, decisions often have to be made on new locations for data collection.Optimizing the installation of data observation stations (or temporary water level stations for storm surge) may be achieved by differentiating multiple potential station sites by their improvement in overall error measures using this interpolation method.

Figure 1 .
Figure 1.Idealized basin and the Peaks function (Equation (8)).(a) Laplace's equation interpolation with 90 random data points (blue circles).(b) Triharmonic equation interpolation with same data points.Numbers along the axes refer to grid points in the 161 × 161 array.

Figure 2 .
Figure 2. Idealized basin and the Peaks function (Equation (8)).Inter-comparison of two-dimensional contour plots between the original data and the 1 st , 2 nd , and 3 rd order interpolation with 90 random data points.Numbers along the axes refer to grid points in the 161 × 161 array.

Figure 3 .
Figure 3.For the idealized basin, the maximum absolute error (MAXE), mean absolute error (MAE), and root mean square error (RMSE) of 4 different order harmonic equation interpolations, and their variation with the increase of internal data points of Peaks: a) MAXE; b) MAE; c) RMSE.

Figure 4 .
Figure 4. Taylor basin solution (the M 2 tide in a semi-infinite rotating channel).The coamplitude (dashed) and cophase (solid) lines for a) analytic results [22] and b) interpolated from 40 random boundary points using 2 nd order interpolation.The green dots are the amphidrome locations from analytic solution.totalnumber of 20,301 nodes.A tidal constituent consists of amplitude A and phase φ, which can be expressed as a complex number,i Ae R E  i  , where R is the real part and E is the imaginary part of the number.We interpolate R and E separately, and after interpolation, we reverse the calculation, 2 2A sqrt R E   , and

Figure 5 .
Figure 5.For the Taylor basin (the M2 tide in a semi-infinite rotating channel), the MAXE, MAE, RMSE and amphidrome distance of four different order harmonic equation interpolations, and their variation with incremental internal data points.a) MAXE, amplitude; b) MAE, amplitude; c) RMSE, amplitude; d) amphidrome distance; e) MAE, phase; f) RMSE, phase.

Figure 6 .
Figure 6.For the Taylor basin, the MAXE, MAE, RMSE and amphidrome distance of five different interpolations, and their variation with incremental boundary data points.M 2 tide in a semi-infinite rotating channel.a) MAXE, amplitude; b) MAE, amplitude; c) RMSE, amplitude; d) amphidrome distance; e) MAE, phase; f) RMSE, phase.

Figure 7 .
Figure 7. Bohai Sea M 2 , S 2 , K 1 , and O 1 cophase lines in degrees (solid line) and coamplitude lines in meters (dashed line).Also displayed are 29 data stations along the coast, and two island stations near the entrance (red dots).The green dots are the amphidrome locations from Fang et al. [24].Table 3.For the Bohai Sea tidal constituents, the location of amphidromes from the Fang et al. [24] model (reference position), our interpolation results, and the approximate spatial difference between the two.Results are for 2 nd order harmonic interpolation only.