Radiational Tides: their double-counting in storm surge forecasts and contribution to the Highest Astronomical Tide.

. Tide predictions based on tide-gauge observations are not just the astronomical tides, they also contain radiational tides - periodic sea level changes due to atmospheric conditions and solar forcing. This poses a problem of double-counting for operational forecasts of total water level during 5 storm surges. In some surge forecasting, a regional model is run in two modes: tide-only, with astronomic forcing alone; and tide-and-surge, forced additionally by surface winds and pressure. The surge residual is deﬁned to be the difference between these conﬁgurations and is added to the local har- 10 monic predictions from gauges. Here we use the Global Tide and Surge Model based on Delft-FM to investigate this in the UK and elsewhere, quantifying the weather-related tides that may be double-counted in operational forecasts. We show that the global S 2 atmospheric tide is captured by the 15 tide-and-surge model, and observe changes in other major constituents, including M 2 . The Lowest and Highest Astronomical Tides levels, used in navigation datums and design heights, are derived from tide predictions based on observations. We use our ﬁndings on radiational tides to quantify 20 the extent to which these levels may contain weather-related components.


Introduction
The operational forecast in several countries of storm surge still-water levels is based on a combination of a harmonic 25 tidal prediction and a model-derived forecast of the meteorologically induced storm surge component. The forecast is based on the "non-tidal residual", the difference of two model runs with and without weather effects. This is linearly added to the "astronomical prediction" derived from lo-30 cal tide-gauge harmonics (Flowerdew et al., 2010). This approach is taken in the UK because the complexity and large range of the tides is such that it has historically been difficult to model them to sufficient accuracy. The same method was applied in the Netherlands until 2015 when improvements 35 to the local surge model DCSM-v6 made it unnecessary (Zijl et al., 2013). It is still in use operationally in the extra-tropical US, where results of the SLOSH surge model are added to local tidal predictions (National Weather Service, 2018); similarly in Germany, using the BSHsmod model (BSH, 2018); 40 and is also used in the new Aggregate Sea-level Forecasting under evaluation in Australia, which also incorporates sealevel anomalies from a global baroclinic model (Taylor and Brassington, 2017).
There are several possible sources of error in this proce-45 dure. The purpose of the combined tide-and-surge model is to capture the well-documented non-linear interactions of the tide and surge. (e.g. Proudman, 1955). Yet the forecasting procedure assumes that the non-tidal residual may be added linearly to a gauge-based tide prediction. There is also an as-50 sumption that the tide-only model and the harmonic prediction from the gauge are equivalent. In fact, the harmonics at the gauge will also be affected by the weather, so there is the potential for double-counting of radiational (weather-related) tidal constituents. 55 In section 2, we show that the double counting of radiational tides has a potential contribution to forecasting errror not just on long time scales (through S a , S sa ) but also on a fortnightly cycle due to variations in S 2 and in the phase of M 2 . We also show that the assumption of non-linearity may 60 introduce errors if phase predictions disagree between model and observations. Specific radiational tides have been studied using response analysis, for example the solar-diurnal S 1 by Ray and Egbert (2004), and semi-diurnal S 2 by Dobslaw and Thomas (2005). In section 3 we look at more constituents, and demonstrate that the atmospheric tide at S 2 may be observed in the GTSM 5 model.
Highest and Lowest Astronomical Tide (HAT and LAT) are important datums used for navigation, and are calculated from tidal predictions. In section 4 we use the model predictions to quantify to what extent HAT and LAT are influenced 10 by weather-related tides, and show that in many places several cm of what is reported as HAT is attributable to periodic weather patterns.
There are other contributors to water level, including steric effects and river flow, that will also create differences be-15 tween the tide gauge and the forecast water levels, particularly seasonally and which may be out of phase with the atmospheric contribution. The problem of double-counting of periodic changes does not arise if they are omitted from the surge model entirely, but they may contribute to HAT and 20 LAT calculations. These effects are not included in this study.

Surge forecasting
The current procedure for forecasting total water level in the UK is as follows: 1. Run a barotropic shelf model (CS3X, currently transi-25 tioning to NEMO Surge ) in tide-and-surge mode, forced by an ensemble of wind and pressure from the current weather forecast to give timeseries M s (x, t) at each location x. Also run the shelf model in tide-only mode, to get M t (x, t). Get the 30 residual from these models, M r = M s − M t .
2. At individual tide-gauge locations, derive a tide harmonic predictionG(x g , t) based on past records. This is assumed to be more accurate locally than the model tide.

35
3. Forecast the total water level F at each location as model residual plus gauge harmonic prediction, 4. Finally, it has been proposed (Hibbert et al., 2015) that the forecast could apply various "empirical corrections" 40 to nudge the forecast towards the observed level G based on the mismatch of the peak tide over the last few days. However no formal correction schemes have been implemented.

45
Similar procedures are implemented elsewhere in the world, so in this paper we replace the regional models with GTSM. This is the forward Global Tide and Surge Model developed at Deltares, on a base of Delft-FM (Flexible Mesh) (Verlaan et al., 2015.;Irazoqui Apecechea et al., 2018 (In review)). 50 The version used in this paper has resolution from around 50 km in the open ocean to around 5 km at the coast. We ran the model in two modes, tide-only M t and tide-and-surge M s . The atmospheric forcing used was the ECMWF (European Centre for Medium-Range Weather Forecasting) ERA-55 Interim 6-hourly reanalysis (Dee et al., 2011), downloaded at 0.25 • resolution but from a spherical harmonic equivalent to ∼ 0.75 • . Validation of the major tidal coefficients has been favourable, and although the model under-predicts the effect of tropical cyclones, due to coarse temporal and spatial res-60 olution in the weather reanalysis, most surge events are captured. We make the assumption that tropical cyclones at any given location are sufficiently rare that the tidal coefficients fitted over a year should not be very different if those surges are underestimated. Due to limitations of data storage and 65 post-processing the output from the model was only saved at high frequency at all grid points for one month (Jan 2012) and a subset of coastal points for the year 2013. All runs were preceded by 11 days spin-up. constituents Harmonic analysis (Pugh and Woodworth, 2014, Chapter 4) gives a tidal predictionG as:

Harmonic analysis and selection of tidal
where Z 0 is the mean of the gauge data and the amplitudes 75 A n and phases g n are associated with the tidal constituents with astronomically-determined frequencies σ n . f n (t) and u n (t) are nodal modulations to amplitude and phase, applied in order to allow for the 18.61 year nodal cycle and 8.85 year longitude of lunar perigee cycle. V n are the phases 80 of the equilibrium tide, which we take as for Greenwich, using UTC for all times. Throughout this paper an overhead tilde indicates "time series derived from harmonics", as the shape is reminiscent of a sine wave. The choice and number of tidal constituents determined 85 by harmonic analysis are typically chosen according to the length and frequency of data available. In this paper we use 62 harmonics where there is one year's data, as listed in table B1. To derive harmonics from the global model from only 1 month's data, we use 26 independent primary constituents, 90 and a further 8 related constituents. We will useM s andM t to indicate harmonic prediction time series from the tide-andsurge model and tide model respectively.

Quantifying the effect on forecast of
double-counting radiational tides 95 A significant source of error for this method is that a tidegauge is measuring the total water level, and hence the harmonic predictionG includes all wave, steric, river levels and surge effects. This is not therefore a prediction of the astronomical tide alone. Steric, wave and river effects are omitted by the barotropic model, but M s does include periodic radiational effects, which may be double-counted. We can test a minimum effect of this double-counting purely within the 5 model by usingM s , the harmonic prediction of the model including surge, as a proxy for the harmonics of the observations at gauges. Then the forecast procedure can be estimated as M r +M s .
To estimate ∆, the error in this model forecast, we can once again use the model, assuming M s ≈ G. Hence ∆ = M s − (M r +M s ) = M t −M s . That is, the minimum error from the current forecast procedure is equal to the error in the harmonic prediction from the model including surge at esti-5 mating the tide-only model, figure 1(top). There are several striking features here, annual cycles peaking around March in the Arctic, January in South-East Asia, and June in Europe. Fortnightly cycles occur almost everywhere, with amplitudes of several cm. We will examine the causes of these 10 below.
If it were possible to avoid the double-counting, and provide astronomical tidal harmonics for the observations, the prediction would instead be equivalent to M r +M t , and the error would become as shown in figure 1(bottom). Since we are using the model as proxy for observations, if the harmonic prediction were an exact reproduction of the tide-only model then ∆ = 0. In practice ∆ < 5 cm at most UK sites and the monthly cycle has gone, but in the Bristol Channel there is still an error of 20 around 50 cm, indicating that the 62 harmonic constituents are not capturing all of the model tide, and further shallowwater constituents may be required. This is consistent with the conclusions of Flowerdew et al. (2010), who found an "average [across UK ports] rms error [in harmonic prediction 25 of a tide-only run] of 7 cm with a maximum value of 29 cm at Newport, in the Bristol Channel", using 50 constituents on the CS3X model. M 2 has a period of 12.42 hours and S 2 exactly 12 hours. They move in and out of phase with each other twice in a lunar month, resulting in the spring-neap cycle. A small change in phase to S 2 harmonic would result in a change of which days it is in phase with M 2 , and hence a substantial change 35 in total tidal amplitude at a given date. For example, near Avonmouth in the Bristol channel S 2 derived from M s has an amplitude 3.5 cm greater than S 2 derived from M t , however there is a phase change of around 3.5 • , so the tide arrives 7 minutes later. Figure 2 shows how this and smaller changes 40 in M 2 account for differences betweenM t andM s of up to 5-8 cm, on a fortnightly cycle between these limits. This can account for about half the error in forecasted high-water at Avonmouth, which varies between 5 and 20cm on a fortnightly cycle (Byrne et al., 2017, figure 4). Similar variation 45 in error of the forecast was seen by Flowerdew et al. (2010).

Quantifying surge-forecasting error due to disregarding non-linearity
The forecasting approach of linear addition of a non-linear model residual to a harmonic prediction, F = M r +G, can 50 also cause errors. Disagreements in phase between the model tide M t and harmonic prediction from the gaugeG affect the forecast of an individual surge event.
Consider a simplified example, where the tide can be modelled by a single constituent, M t = A cos(σt). Suppose there 55 is a storm-surge in which there is an uniform additional water level A s and an advancement of the tide of t = δ, so the tideand-surge model is M s = A s + A cos(σ(t + δ)). As before, the model residual is given by Suppose the harmonic prediction at the gauge agrees in 60 amplitude to the tide-only model, but has slightly different phase:G = A cos(σ(t + )).
The skew surge is defined as the difference between the maximum water level, here max(M s ), and max(G). The error in the skew surge forecast is E = max(M r +G) − 65 max(M s ). Substituting in and assuming phase changes are small, we find A s cancels out and can show analytically that E ≈ A (cos(σ ) − cos(σ(δ + )) + cos(σδ) − 1) . This is illustrated in figure 3, with A = 3 m, σ = 2π/12.42hr −1 (M 2 ), and the surge advancing the tide by 70 δ = 30 min. The residual M r is decreasing during High Water due to the advanced tide. So if the observed harmonics have High Water later than the model ( =5 min), the forecast skew surge is underestimated by 3 cm. If the observed harmonics predict High Water earlier than the model ( =-5 min), 75 the forecast skew surge is overestimated by 3 cm.
Although in practice there are more constituents, a similar relationship will still hold in a small window about each high tide. Where there are frequent surges with consistent effect on the tidal phase we would expect to have the same sign 80 as δ, as the gauge registers water levels more like the tideand-surge model than the tide-only model and the harmonic predictions would follow suit. 3 The difference of specific harmonics Figure 4 shows the vector difference in individual constituents between tide-and-surge and tide-only models run for 2013, along the coast globally. With some exceptions in the Arctic and Antarctic, the effect on S a is around 5-20 cm, 5 with around half that effect on S sa , although in the Indian Ocean there is a change to S a only. Since the model was only run for one year, S a may not be representative of all years, but figure 4 indicates typical changes. In the Baltic, the seasonal change is wind-forced, but elsewhere it is con-10 sistent with the annual and semi-annual cycles in sea-level atmospheric pressure (Chen et al., 2012). MS f is affected by the surge component, as a side effect of the interaction between M 2 and S 2 . This is because MS f is the fortnightly constituent which arises from the combination 15 of M 2 and S 2 , with a speed equal to the difference of their speeds. MS 4 is the counterpart to this, with a speed equal to the sum of the speeds of M 2 and S 2 (Pugh and Woodworth, 2014). Less explicable is the effect on M m and M f , but it may be due to insufficient separation with MS f over 20 a relatively short record. Another possibility is that non-tidal power in the tide-and-surge model is leaking into M m and M f estimates. Eliminating this would require a many-year model run.
The diurnal constituents K 1 and O 1 are affected by less 25 than 5 cm, and are only changed regionally in the Antarctic. S 1 however is everywhere less than 0.1 cm in the tide-only model, but with the surge model peaks at 0.5 cm in northern Australia, the broadest regional effect being has 0.2-0.3 cm, in South-East Asia, consistent with the findings of (Ray and It may come as a surprise that constituents such as M 2 , which has a purely lunar frequency, could possibly be affected by the weather. There is a very small atmospheric tide at M 2 , peaking at the equator at about 0.1 mbar (Schin-35 delegger and Dobslaw, 2016). But more significant is the non-linear interaction of surge and tide. The surge may consistently advance the phase of the tide during low pressure events and certain wind configurations. A high pressure system could delay the phase of the tide, but there is asymmetry 40 between these events, so there is a net bias on the phase when the weather is included.
The effect on higher order constituents is everywhere less than 5 cm. The maximum difference in the UK and globally for each constituent is given in Appendix B. In the UK, the 45 constituents affected the most by including the surge are S 2 , S sa , M 2 , S a , M m , MS 4 , MS f and M f , with a maximum change of > 2 cm, and a further 19 constituents change 1-2 cm. Globally, S a and S sa are far more significant, but S 2 , , and MS 4 all 50 change more than 4 cm (somewhere on the global coast). A vector difference of 13 cm in S 2 is seen in north-west Australia.
We tested the stability of these results to the number of constituents fitted, using the list of 115 harmonics usually 55 associated with 18.6 years of data (see supplemtneary information) and found the changes remain to within 0.2 cm.

S 2 atmospheric tide
Some of the difference between harmonics of surge and tideonly models is directly attributable to the atmospheric tides. 60 The global atmospheric pressure field contains S 2 variations with amplitude of about 1.25 cos 3 φ millibars, for latitude φ (Pugh and Woodworth, 2014). GSTM air pressure and wind forcing is taken from the ERA-Interim data set (Appendix A), and the ocean response to that forcing at S 2 is contained 65 in the difference between harmonic predictions of the M s and M t model runss (figure 5). It is consistent with response analysis based on the S 2 tides seen in ECMWF reanalysis data (figure 2, Dobslaw and Thomas, 2005), and in a 2-layer model forced by 8 constituents (figure 1b, Arbic, 2005). The 70 6-hour sampling prevents ERA-Interim forcing from capturing the S 2 atmospheric tide correctly (Dobslaw and Thomas, 2005), but the analysis in this paper is self-consistent with the forcing used.

Highest Astronomical Tide and Lowest Astronomical 75 Tide
The Highest Astronomical Tide (HAT) is used internationally for flood-forecasting references levels and in navigation for clearance under bridges. HAT can be used in structural design alongside skew surge as an independent variable for 80 determining return period water levels. Lowest Astronomi- cal Tide (LAT) is also an important parameter, recommended for use as the datum on navigation charts (IHO, 2017). Once the phases and amplitudes A n and g n are known,G(t) is fully determined for all time by equation (1), and the future HAT and LAT are given by max(G(t)) and min(G(t)). But 5 because of the overlap in phase of the forcing between the constituents, and the f n and u n nodal modulations, it is not trivial to write HAT or LAT algebraically. They are therefore determined by inspection of the predicted tides, preferably over a 18.6 year nodal cycle. Figure 6a shows the range, HAT 10 minus LAT, when we do this by synthesising a predicted tide at 15 minute intervals over 18.6 years, globally. Radiational effects are omitted from this figure, which is based on a tideonly run. Since the GTSM data was limited to 1 month, it uses only 34 constituents, therefore omitting S 1 and the long-5 period contributions to HAT and LAT. An approximate calculation of Range = 2(M 2 +S 2 +O 1 + K 1 ) is occasionally used (e.g. Yotsukuri et al., 2017), but the error due to this can be over 1 metre (figure 6b). N 2 is a significant contributor, at about 20% of M 2 in many sites 10 worldwide. A few tens of centimetres are accounted for by the omission of the nodal modulations, and there are also the shallow water constituents at the coast. Figure 6c shows the effect on HAT and LAT of including surge in the GTSM. Coastal locations are shown and 15 62 constituuents used. In many places round the world the HAT is higher when the tide-and-surge model is used. So the observation-based HAT has been raised by some radiational component. But in most of the UK, the HAT goes down when the tide-and-surge model is used to generate the tidal predic-20 tions. This is because the peak of the weather-related components does not coincide with the maximum astronomical effects alone.This implies that since the tide-gauge predictions include surge, the observation-based HAT in the UK is actually about 10 cm lower than true astronomical-only tidal 25 height.
LAT tends to move the opposite way, so in most places the maximum tidal range is increased by using the tide-andsurge model. That is, the true astronomical-only tidal range is slightly less than that quoted from harmonics based on pre-30 dictions. In Scotland, (just above Liverpool in figure 6c) both LAT and HAT go down when the surge model is used to generate the tidal predictions, so the quoted LAT and HAT are actually about 10 cm lower than astronomical only.
The most extreme changes shown in Figure 6c are in the 35 Arctic and Antarctic, and should be interpreted with some caution as these areas are the least well represented in the model.
In places with small tide, seasonal signals may be dominant and they may be important to include for practical pur-40 poses. For example along the French/Italian coast from Mallorca to Sicily there is about a 7 cm increase in HAT and 3 cm decrease in LAT using the surge rather than tide-only model, so a Highest "Astronomical" Tide based on predicted tide from observations actually contains about 7 cm due to 45 seasonal winds.

Conclusions
There are substantial changes in tidal constituents fitted to tide-only and tide-and-surge model results. Even constituents with purely lunar frequencies, including M 2 , may be af-50 fected by the surge, perhaps owing to asymmetry in phase changes of the tide under high and low pressure weather systems.
Some effects of the weather on tides are double-counted in the forecast procedure used in the UK, where model resid-55 uals are added to gauge-based tide predictions. Even if the model were perfect, the minimum error from the current forecast procedure would be at least the error in the harmonic prediction including surge at estimating the tide-only model. If 62 constituents are fitted, this has a standard deviation of 60 20 cm at Avonmouth and 4-10 cm at most other UK gauges. 5-8 cm of the error at Avonmouth is due simply to a small change in phase of the S 2 harmonic. Further errors in total water level and skew surge arise directly from the linear addition of the harmonic prediction to the non-linear residual, 65 particularly where there is a phase difference between model and gauge tidal harmonics.
Understanding and quantifying these errors is extremely important for forecasters, who will often need to advise or intervene on the expected surge risk, often based on a direct 70 comparison between observed residuals and the forecast nontidal residual. Where, for example, such a comparison may lead to the observed residual falling outside the bounds of an ensemble of forecast non-tidal residuals, the forecaster may significantly (and potentially incorrectly) reduce their con-75 fidence in the model's estimate of surge if they are unaware of the additional errors associated with the harmonic tide and whether or not they have been addressed within the ensemble forecast's post processing system. For comparison, across the UK tide-gauge network, short range ensemble forecast RMS spread is of order 5-10 cm (Flowerdew et al., 2013). It is noted that, in the UK, the majority of coastal flood events occur around peak spring tides (Haigh et al., 2015), where the 5 sensitivity to any errors in the M 2 -S 2 phase relationship is arguably at its highest. The atmospheric tide, at S 2 is present in the ERA Interim forcing, and the ocean response to it, with amplitude about 1-5 cm, can be seen in the difference between the model re-10 sults with and without surge. There is hence an argument for including an atmospheric tide forcing in a "tide-only" model, and this is being explored by Irazoqui Apecechea et al. (2018 (In review)). In this case, care would need to be taken to omit the direct atmospheric tide forcing in the tide-and-surge version, to avoid a different form of double-counting.
The estimates of Highest and Lowest Astronomical Tide are influenced by radiational tides. HAT and LAT are most readily calculated by inspecting long time series of predicted tides, and if observation-based, these predictions will include weather-related components. In most places globally this re-10 sults in HAT being calculated as higher than the strictly astronomical component, and LAT being lower, however the opposite is true in the UK. The effects are of the order of ∼ 10 cm.
For many practical purposes it is correct to include pre-15 dictable seasonal and daily weather-related cycles in the HAT and LAT. However the separate effects should be understood, as the radiational constituents may be subject to changing weather patterns due to climate change. It is also important not to double-count weather effects, if HAT or LAT are used 20 in combinations with surge for estimating return-period water levels.
These considerations about HAT would also apply (proportionally less) to other key metrics such as mean high water. 25 Data availability. Model data is available from the authors on request.

Appendix A: Ordering of model sites around the coast
The coastal points in the model output are spaced roughly every 80 km, and also wherever a tide gauge is situated, ac-30 cording to the GESLA data set (Woodworth et al., 2017). Due to automatic procedures to select output sites, a few may be incorrectly sited at model dry sites -these are clearly seen in plots as lacking sufficient high-frequency variability. The along-coast plots are ordered approximately from 35 west to east around the world coastline, starting and ending at Alaska. The order is indicated in figure A1.
The algorithm for coastal order is as follows: 1) Define a single global coastline polygon. This is done using the GSHHG (Global Self-consistent, 40 Hierarchical, High-resolution Geography) data set (Wessel and Smith, 1996), version gshhg2.3.6 (August 19, 2016, downloaded from www.ngdc.noaa.gov) . We use the coarse resolution, with only Level 1 (coastline) and Level 6 (Antarctic Ice Shelf), although consistent results for this technique 45 can be obtained including enclosed lakes. To merge the separate landmasses and islands into a weakly simple polygon, topologically equivalent to a disc, we start with a single landmass and add others in turn using pairs of identical edges as "bridges". We start with the main landmass of Eurasia 50 L 1 , and find the closest vertex l to a vertex p from any of the remaining polygons [P 2 , ...P N ]. Suppose p belongs to polygon P j . Then we add P j to L 1 using two new edges − → lp and − → pl , to give a new merged polygon L 2 . The vertices of L 2 are then [L 1 (1 : l), P j (p : end, 1 : p − 1), L 1 (l : end)]. 55 Now repeat, searching for the nearest point in L 2 to any vertex in the remaining polygons [P 2 , ..., P j−1 , P j+1 , ...P N ]. It is necessary for all initial polygons to be defined in the same sense (anticlockwise). If inland seas (Level 2) are included, they should be defined clockwise. The GSHHG data is con-60 sistent with this definition. Distance for nearest points are defined as arc-length on a sphere. This technique has the benefit of tending to group island chains together in a consistent order. It cannot produce crossing edges. Because polygons are added in distance order, is-65 lands near continents are added to their neighbouring coast, and remote mid-ocean islands tend to be clustered, and attached to the nearest continent. The coasts of the Pacific, Atlantic and Indian and Arctic Oceans are all treated clockwise. Antarctica is attached across Drake Passage and ordered 70 Westward. Nearby locations across narrow islands (particularly Sumatra), isthmuses (Panama), and straits (Gibraltar) may be widely separated in the order. But neighbouring points in the order can be expected to have fairly smoothly varying oceanography, with the "bridges" often, although not 75 necessarily, approximating shoals.
As a final step we adjust the starting point of L 2 to be in Alaska, for convenience of mapping.
2) Rank the coastal points according to the nearest point on the global polygon.

80
Having defined this coastal order, we can apply it to any coastal data set, for example tide gauges. We number the vertices [1, ..., K]. For each of the gauge locations T we find the nearest vertex k, and then rank the gauges according to T k . In the event of gauges being much closer than the resolution 85 of the vertices, a quick method for refinement is to linearly interpolate with extra vertices along polygon edges. Some problems may also occur with islands not in the coarse resolution data, which will tend to jump to the nearest coast.
A further advantage here is that having defined the coastal 90 polygon, the same order can be applied to different data sets and models, leading to closely comparable along-coast plots. Table B1 lists the constituents used in this paper. For the 1 month model run, related constituents are used, and we fit 34 95 constituents with only 26 independent terms. 62 consituents are used for the 1 year run. The list of 115 usually applied to 18.6 year data is used only as a check on stability of the result in section 3 and is provided as a supplementary spreadsheet. sults in this paper first appeared as an internal National Oceanography Centre report Williams et al. (2018). We thank Martin Verlaan of Deltares and Clare O'Neill for model development work, Phil Woodworth, Richard Ray and two anonymous reviewers for helpful suggestions during the final preparation of the paper.