Non-assimilated tidal modeling of the South China

The tides in the South China Sea were simulated using an established tidal model, with the purpose to evaluate if non-assimilated modeling of the area is feasible. Simulations were done for the locally dominating diurnal (K1) and semi-diurnal (M2) tidal constituents, and the model was shown to provide reasonably accurate results in terms of both elevations and levels of dissipation. However, this was only the case when a realistic tidal conversion parameterization was included in the model, and it is suggested that tidal conversion is a missing process in other model efforts of the area. Compared to observations, the modeled dissipation levels were slightly overestimated when integrated over the entire domain, and far larger in the model at topography with a slope which is supercritical for the baroclinic tidal waves. A crude, empirical correction of the tidal conversion rates at supercritical topography is suggested and implemented in the model and shown to improve the model results in terms of both elevations and dissipation rates. It is concluded that the presented model set up is suitable for investigations of how perturbations, e.g., future sealevel rise, will affect the tidal dynamics in the South China Sea. & 2013 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/3.0/).


Introduction
Our ability to model tidal dynamics, both regionally and globally, with high resolution and relatively high accuracy has improved significantly over the last decade (e.g., Egbert et al., 2004;Green, 2010;Müller et al., 2011;Pelling et al., 2013;Green and Nycander, 2013). One key model component to achieve good accuracy is to implement a parameterization for the conversion of tidal energy from barotropic to baroclinic modes (Egbert et al., 2004;Green and Nycander, 2013), especially in global simulations or simulations of areas where the conversion is significant. Recently, Zu et al. (2008) described the tidal dynamic in the South China Sea (Fig. 1a) in great detail, and the reader is referred to their paper for a detailed dynamical description. They used the Oregon State University Tidal Inversion Software (OTIS) to assimilate satellite altimetry into a numerical solution for the tides, but they did not specifically include tidal conversion. They concluded that the solution without assimilation seriously underestimates the dissipation in the area, especially in the Luzon Strait. As a consequence, the tidal amplitudes in the South China Sea itself are overestimated in the non-assimilated runs (see Figs. 2 and 5 in Zu et al., 2008). With data assimilation included they achieved a root mean square difference between the model and nonassimilated observations of some 4-5 cm for the M 2 amplitudes. However, for future forecasting and perturbation response analyses, using assimilated models is dubious as we do not know a priori how the data will change. It is thus necessary to be able to represent an area with a forward solution only, and to include those processes which may or may not become important. Here, we therefore ask if OTIS can accurately reproduce the tides and tidal conversion in the South China Sea without assimilation if a tidal conversion parameterization is included.
The South China Sea (SCS in the following; see Fig. 1) is made up of a central deep basin and shallow shelf seas in the North and Southwest. It links with the Java Sea in the South, with the Sulu Sea through several narrow channels between the Philippine Islands, and directly with the Pacific through the very energetic Luzon Strait south of Taiwan. The dynamics in the area is complicated and influenced by the seasonal monsoon via modifications of the location of the Kuroshio current (e.g., Hu et al., 2000;Jan et al., 2012), and complex topography, especially in the Luzon Strait, leading to generation and scattering of internal waves with a variety of frequencies (see, e.g., Egbert and Ray, 2000;Zu et al., 2008;Alford et al., 2011;Li and Farmer, 2011;Jan et al., 2012, and references therein). The tides in the SCS enter mainly through the Luzon Strait (Zu et al., 2008), leading to a tidal system with several amphidromic systems, especially on the shelfs.
It is estimated that some 50-75 GW of the M 2 energy is dissipated in the SCS (Egbert and Ray, 2000), with a significant fraction being lost in the Luzon Strait (Zu et al., 2008;StLaurent et al., 2011). The conversion rate in the strait has been estimated to be around 1 W m −2 , with the baroclinic energy flux into the SCS exceeding 60 kW m −1 (Alford et al., 2011). The conversion of barotropic energy into baroclinic waves acts as a drag on the barotropic flow, and the drag force can be parameterized as an additional, linear, friction term in tidal models (see Green and Nycander, 2013, for an evaluation of three conversion schemes). Tidal conversion parameterizations focus on propagating internal (tidal) waves generated at topography with subcritical slopes. At supercritical topography, i.e., where tan ðθÞ 4 ðN 2 −ω 2 Þ=ðω 2 −f 2 Þ (θ is the slope of the bathymetry, N is the buoyancy frequency, ω the tidal frequency, and f the Coriolis parameter), the actual levels of energy conversion are less clear, and a few different approaches can be found in the literature. Nycander (2005) simply sets his conversion rate to 0 at supercritical topography, Green and Nycander (2013) somewhat arbitrarily cap the conversion rate there at 1 W m −2 , whereas specific conversion scheme for supercritical topography is presented by Balmforth and Peacock (2009). Here, we will investigate how much energy is lost in the model at supercritical topography if we use the scheme given by Nycander (2005), and compare that to computed dissipation levels from the heavily assimilated solution for the area (see http://volkov. oce.orst.edu/tides/YS.html for details). We refer to this as the ATLAS in the following as the solution is part of the OTIS ATLAS for the Western Pacific. More details can be found in http://volkov.oce.orst. edu/tides/atlas.html. The objectives in this paper are thus to (i) investigate if we can reproduce the tides and the tidal dissipation in SCS without assimilation, and (ii) estimate the tidal conversion rates at supercritical topography in the model and in the ATLAS. We continue in the next section with a description of the tidal model, its set-up and the conversion parameterizations used. Sections 3 and 4 contain the results, and we finish the investigation with a discussion in Section 5.

Model set-up
OTIS has previously been used in a number of investigations with and without data-assimilation (e.g. Egbert et al., 2004;Zu et al., 2008;Green, 2010, to mention but a few). Here, we only use the forward stepping part of the model, i.e., we do not do any assimilation of data, and thus solve the (linearized) shallow water equations with rotation where U ¼ uH is the depth-integrated volume transport given by the velocity u multiplied by the water depth H, f is again the Coriolis parameter, k the vertical unit vector, η the tidal elevation, η SAL the self-attraction and loading elevation, η eq is the equilibrium tidal forcing, C d ∼1 Â 10 −2 is a dimensionless bottom drag coefficient, u a is the tidal velocity from all constituents combined, and Cðx; yÞ is the internal wave drag tensor. The model was set up with a horizontal resolution of 1/301 in both longitude and latitude using the bathymetric database in the ATLAS for the area, albeit at a reduced domain size here. The bathymetry initially came from GEBCO2008 (http://www.gebco. net/) and was averaged to 1/301 resolution (see Fig. 1). The ATLAS itself is a tidal solution which has assimilated satellite altimetry, and is a further development of the TPXO database (Egbert and Erofeeva, 2002). The version here is an updated version of the simulations presented by Zu et al. (2008) and is accurate within a few centimeter of data from tide gauges (see Egbert and Erofeeva, 2002;Egbert et al., 2010, for details about the methodology).
Simulations were done using forcing (and subsequent harmonic analysis) for the M 2 and K 1 constituents. The forcing consisted of the tidal potential for each of these constituents applied over the entire domain, and of prescribed elevation amplitudes and phases at the open boundaries taken from the ATLAS. Our model output consists of surface amplitudes and phases, and the tidal transport amplitudes and their phases for each of the constituents under evaluation. The ATLAS solution was also used to verify the present model simulations, along with a comparison of the coastal tide gauge (TG) data presented in Zu et al. (2008), to which a reader is referred for more details.

Conversion parameterizations
Following Green and Nycander (2013), the parameterization by Nycander (2005) is given by where The Green's function G is defined by where I 0 is a modified Bessel function of the first kind, and the cutoff length a is given by where β ¼ 1:455 is a numerical coefficient, N B is the buoyancy frequency at the seabed, and H is the bottom topography (increasing upward). The generation of internal waves by topographic features with a larger length scale than the horizontal wavelength of the first internal wave mode is very weak (Smith and Young, 2002), and using the "filtered" Green's function given by Eq. (6) is an approximate way of taking this into account. Note that due to the frequency dependency of the conversion parameterisation separate runs will have to be done for each of the different tidal constituents. Also note that C is a spatially varying 2 Â 2 tensor (Green and Nycander, 2013) whereas C d is spatially homogenous.

Simulations and computations
The stratification used to calculate the buoyancy frequency in the conversion scheme came from the WOCE database, which has a resolution 11 Â 11 (Gouretski and Koltermann, 2004). It was interpolated linearly to the model grid, but fails to represent the stratification in some shallow regions. These are simply left without any conversion in the following, but are generally shallow and far away from large topographic features and should thus provide only a minor contribution to the conversion.
By taking (1)Áu þ ηg (2), introducing the energy density E ¼ 0:5½Hu 2 þ gη 2 , assuming that ∂E=∂t ¼ 0, and taking the timeaverage, we arrive at the well known expression for the tidal dissipation D (in W m −2 , see, e.g., Egbert and Ray, 2001): Here the work rate W and energy flux P are defined as in which 〈 〉 denote the time-averages. Using ATLAS data or the output from the model it is possible to calculate the dissipation rate for each constituent. Note, however, that using Eq. (8) gives the total dissipation rate, e.g., due to both conversion and bed friction.

Tidal amplitudes
The model simulations were tuned to give the lowest possible difference between the M 2 amplitudes from the model and the ATLAS database (see Fig. 2a for details). The bottom drag coefficient, C d , was tuned to minimize the root-mean square (RMS) difference and maximize the variance capture between the model and the ATLAS. We arrived at an optimum valued for C d of 0.01, i.e. about three times higher than what is usually used. A simulation using the optimized values is referred to as "control" in the following, whereas runs labeled "SCn" use a tidal conversion around supercritical topography modified by the factor given in "n" (see Section 4 for details).
The RMS difference between the model and the ATLAS for the control is about 9 cm, and the model captures some 93% of the variance in the ATLAS M 2 amplitude. The corresponding values for K 1 are a 10 cm RMS and 91.2% variance capture when compared to the ATLAS (see Figs. 2 and 3).
When comparing the M 2 amplitudes to the TG data, the agreement is less good, with an RMS of some 17 cm. However, this large error is down to 5-10 stations (depending on constituent) where the model and TG data disagree with more than one standard deviation (see Fig. 2a and c). These TG stations may be in areas where the local topography is not accurately captured in the bathymetric database, and if they are removed from the comparisons the RMS difference between the model and the TG data is 7 cm for M 2 . The same comparison for the K 1 data gives a 17 cm RMS, whereas we get a 6 cm RMS difference for K 1 data points within one standard deviation. As a comparison, Zu et al. (2008) achieve an RMS difference of some 4-5 cm using all TG data, but they do assimilate altimetry into their solution so a better fit is expected. We argue that the model is doing a decent job of recreating both the diurnal and semidiurnal tides in the SCS.
There are obvious discrepancies between the modeled and ATLAS M 2 -amplitudes at the entrance to the Gulf of Thailand, along the coast of Borneo, and in the Strait of Taiwan (Fig. 2). This is due to an offset in the location of the virtual M 2 amphidromic point along the coast of Vietnam (Fig. 3). The model also underestimates the M 2 tide in the Sulu Sea, and the overall averaged difference between modeled M 2 amplitudes and those in the ATLAS is some −2 cm (i.e., the model underestimates the tidal amplitudes). All of these effects are most likely down to bathymetric effects which are compensated for by the assimilation in the ATLAS solution.
For K 1 there is a large underestimate of the tidal amplitudes northeast of Vietnam, and again an offset amphidrome in the Gulf of Thailand. The averaged difference between the model and the ATLAS is −4 cm for K 1 amplitudes, i.e., again an underestimate. The modeled amplitudes are too large in the Sulu Sea, which again suggests issues with the topography in that area. The central SCS sees a general underestimate of the tidal amplitude, suggesting too large K 1 energy losses in the the Luzon Strait (which is supported by the dissipation rates, see Table 1).
In the control run the discrepancy between the modeled K 1 amplitudes and those found in the ATLAS is relatively small but covers nearly the whole domain, whereas for the M 2 results the differences are larger but confined to embayments. This is most likely due to the M 2 constituent being more energetic than K 1 for the same tidal amplitude, and thus more sensitive to shallower water. This constituent dependency of the tuning points towards the missing process being tidal conversion of a form which is not captured fully by the current scheme. Consequently, since we optimized our runs for the semi-diurnal tide, we obtain a bed friction which is too high for K 1 . In fact, running the model with a friction coefficient of 0.003 improves the K 1 amplitudes but distorts M 2 , suggesting that we underestimate the M 2 dissipation Fig. 3. (a-b) The M 2 and K 1 tidal amplitudes (color; in meters) and relative phases (white contours separated by 601) from the ATLAS. Panels (c-d) are the same as (a) and (b) but show model results.

Fig. 2. (a) M 2 amplitude difference in meters between the ATLAS data and the model shown in color, and the TG stations marked and numbered. Plus-signs mark stations with exact interpolation, circles mark best fit interpolation locations (see text for details). (b) as in (a) but for K 1 . (c)
The M 2 (with pluses and circles for exact and interpolated locations, respectively) and K 1 (marked with stars and diamonds) amplitude difference between the model and the TG data for each station. The dashed (dash-dotted) lines mark 7 1 standard deviation of the differences for M 2 (K 1 ). rates in the path of the tidal wave (e.g., in Luzon Strait) despite the conversion scheme in runs with C d ¼0.003.

SCS
The dissipation of M 2 tidal energy in our control simulation amounts to 71.6 GW when integrated over the whole domain and some 20.9 GW of this is lost in waters deeper than 100 m (see Table 1 and Fig. 4). The ATLAS points towards 69.9 GW in total. At supercritical topography, the ATLAS implies a dissipation of 11.6 GW, whereas the model comes up with almost twice that: 20.1 GW. Neither shown, nor discussed further, is the amount of energy lost in water shallower than 100 m, which for M 2 amounts to some with 44.8 GW in the ATLAS. This loss should be dominated by bed friction, and gives a handle on the ratio between the losses of energy due to friction and tidal conversion.
The model overestimates the conversion in the domain, especially near supercritical topography where the model dissipates nearly twice as much energy as is suggested by the ATLAS (cf. Fig. 4). However, running without any conversion at supercritical topography (this is referred to as SC0 in the following), implemented by setting C ¼ 0 there, leads to an underestimated conversion rate of 58.0 GW in total. Furthermore, 3.4 GW is still lost at supercritical topography due to bed friction in this case. This indicates that without the drag induced by the tidal conversion, bed friction compensates (on a basin scale) for the lack of internal wave generation, but it significantly underestimates the local dissipation rates. These results suggest that we could obtain a more realistic modeled conversion if the conversion coefficient was reduced, but not equal to 0, at supercritical topography.
The dissipation values presented suggest that we should use C crit ¼ 0:57C at supercritical topography only, as this is the ratio between the ATLAS and control simulation dissipation rates at supercritical topography. Running the model with this modification does indeed improve the simulations, albeit very slightly: the RMS is still 9 cm and the variance capture has increased to 94.0% for M 2 . Furthermore, the dissipation at supercritical topography has decreased to 14.9 GW, which is significantly better when compared to the ATLAS' 11.6 GW, especially when taking into account that there is 3.4 GW lost at supercritical topography due to bed friction in the model (see the discussion above regarding simulation SC0).
For the dominating diurnal constituent (K 1 ), the model underestimates both the overall conversion (some 10% larger in the ATLAS) and the losses in the Luzon Strait (62% larger), but overestimates it at supercritical topography (the ATLAS dissipation is some 75% of the modeled). The total dissipation rates are comparable to the M 2 losses, albeit with a larger fraction being lost in shallow water for M 2 due to the larger M 2 velocities in shallow water (not shown).

Table 1
A summary of the horizontally integrated dissipation rates in GW from the different simulations. "SC0" refers to computations or model simulations without any conversion at supercritical topography, whereas "SC0.57" and "SC0.74" are results from when the internal wave drag coefficient C was set to 0:57C control for M 2 and 0:74C control for K 1 at supercritical topography (C control is the conversion coefficient used in the control run).  The subcritical to supercritical K 1 dissipation ratio between the ATLAS and the model results is 0.74, indicating that the impact of supercritical topography is weaker for diurnal than semidiurnal tides. Indeed, simulation SC0 gives an enhanced K 1 -dissipation compared to the control, although the conversion rates at supercritical topography are obviously reduced in the SC0 run. This enhanced dissipation over the domain is a surprising result, but can be explained by too much K 1 -energy entering the SCS itself because of significant areas in the Luzon Strait with supercritical topography where the energy loss in the barotropic wave is reduced. As seen from Table 1, the dissipation in the Luzon Strait is about 1/3 in simulation SC0 compared to that in the control. Repeating the simulations, but with the conversion coefficient for K 1 set to 0.74 (simulation SC0.74) of its control value at supercritical topography gave 56.1 GW in total, of which 25.1 GW is lost at supercritical topography. SC0.74 thus underestimates the total dissipation rate but provides a more accurate supercritical dissipation rate. It also improves the model accuracy, giving a 9 cm RMS for K 1 with a 92.1% variance capture.

The Luzon Strait
The M 2 constituent loses 9.7 GW in the Luzon Strait in the model, whereas the ATLAS points towards 13.4 GW. The K 1 dissipation rate is 17.6 GW in the Luzon Strait according to the ATLAS, but only 10.7 GW in the model. Both the modeled M 2 and K 1 dissipation rates are obviously too small in the strait. There are possible explanations for this, with the most likely being erroneous topography in the GEBCO database for the area or an incorrect stratification for the strait itself in the global climatology. Correcting for supercritical topography does nothing to improve this because we generally overestimate dissipation at supercritical topography and we are already on the low side for the Luzon Strait. This also highlights the importance of having accurate topography in complex regions in tidal model simulations.

Discussion
An established tidal model was used to simulate the dominating diurnal (K 1 ) and semidiurnal (M 2 ) tidal constituents in the South China Sea, without data assimilation but with a realistic tidal conversion scheme. The model accuracy is not as high as that in Zu et al. (2008), but they used heavy data assimilation, which makes it dubious to use their model formulation to look at any effects of perturbations as we do not know how the assimilated data will change. The present model uses the linearized equations without horizontal diffusion. Simulations were initially done with these terms included but these runs did not provide better results than those presented here, and for computational reasons we opted for the linearized set (see Egbert et al., 2004, for a discussion).
The stratification database used is a relatively low-resolution global climatology, but we feel that it is accurate enough for this type of study (see Nycander, 2005, for a further discussion). The main lack of data in the WOCE database can be found around shallow waters, where the conversion should be relatively small anyway because of the limited topographic gradients. Furthermore, in sensitivity runs with the conversion multiplied by factors 1.5 or 0.5 (chosen arbitrarily) the accuracy of the runs did not improve and led to large deviations in terms of the dissipation field. Of more significance for the model accuracy is the bathymetric database used. Initial tests with the widely used Smith and Sandwell global bathymetry (see Smith and Sandwell, 1997, the version tried here was v14.1) showed that it was very difficult to get accuracies near those we obtained with GEBCO and described earlier. This could be down to an overall difference in water depth in the two databases: Smith and Sandwell's bathymetry is on an average 42 m, or some 2%, deeper than GEBCO (see Fig. 5). The RMS difference between the two datasets, however, is 210 man extraordinary difference for such a shallow domain. Without arguing that one is more accurate than the other, we opted for the database which provided us with the best possible results, which fortuitously was the GEBCO-based bathymetry used by Zu et al. (2008).
There is a justification for our higher than normal value for the bottom drag coefficient, C d . The SCS can be described as a highly dissipative Helmholtz resonator because of the increase in tidal elevation for K 1 over the Luzon Strait, whereas the M 2 elevation drops (Zu et al., 2008). Sutherland et al. (2005), investigating the Juan de Fuca Strait, justify the necessary use of large bed friction coefficients with their area being highly dissipative. Another motivation is that there may be large areas of complex topography not fully resolved in the bathymetric database. The use of an increased drag coefficient would then compensate for the lack of drag-both frictional and tidal conversion-there. Also, the semi-diurnal internal tide is most likely enhanced (compared to the diurnal ride) in Luzon Strait, because the spacing of the two ridges in the Luzon Strait is comparable to the wavelength of the internal M 2 tide (Li and Farmer, 2011). The present conversion scheme may thus underestimate the dissipation rate of the semi-diurnal tide there (but not of the diurnal). The default parameterization in the model is that presented by Zaron and Egbert (2006), which when tried here gave results in which the dissipation rates were strongly overestimated, especially at supercritical topography (hence the capping done a posteriori in Green and Nycander, 2013). This is again an argument to use Nycander's (2005) conversion scheme in tidal models, albeit with the caveat that it does not allow for tidal conversion poleward of the critical latitude. Zu et al. (2008) accurately simulate the tidal dynamics in the SCS using assimilation but without any tidal conversion. Here, we set out to investigate if we can reproduce the dominating tides in the SCS without assimilation by including tidal conversion, and we argue that we canwithin limits. Our accuracy is not as good as that in Zu et al. (2008), but still good enough for us to argue that using GEBCO in a combination of OTIS with Nycander's (2005) conversion scheme is a valid model system which can be used for perturbation simulations, e.g., the impact of climate change or future sea-level rises on the tides in the SCS. There are growing demands to assess the impact of the future sea-level rise on coastal environments, especially in low-lying countries and many areas in southeast Asia. Furthermore, there is quite a large interest in the tidal community regarding changes in dissipation over geological time scales (e.g. Egbert et al., 2004;Green and Huber, 2013), and we hope that the present investigation improves the quality of models used for both types of investigations.