Two-dimensional numerical simulations of shoaling internal solitary waves at the ASIAEX site in the South China Sea

The interaction of barotropic tides with Luzon Strait topography generates some of the world’s largest internal solitary waves which eventually shoal and dissipate on the western side of the northern South China Sea. Twodimensional numerical simulations of the shoaling of a single internal solitary wave at the site of the Asian Seas International Acoustic Experiment (ASIAEX) have been undertaken in order to investigate the sensitivity of the shoaling process to the stratification and the underlying bathymetry and to explore the influence of rotation. The bulk of the simulations are inviscid; however, exploratory simulations using a vertical eddy-viscosity confined to a near bottom layer, along with a no-slip boundary condition, suggest that viscous effects may become important in water shallower than about 200 m. A shoaling solitary wave fissions into several waves. At depths of 200–300 m the front of the leading waves become nearly parallel to the bottom and develop a very steep back as has been observed. The leading waves are followed by waves of elevation (pedestals) that are conjugate to the waves of depression ahead and behind them. Horizontal resolutions of at least 50 m are required to simulate these well. Wave breaking was found to occur behind the second or third of the leading solitary waves, never at the back of the leading wave. Comparisons of the shoaling of waves started at depths of 1000 and 3000 m show significant differences and the shoaling waves can be significantly non-adiabatic even at depths greater than 2000 m. When waves reach a depth of 200 m, their amplitudes can be more than 50 % larger than the largest possible solitary wave at that depth. The shoaling behaviour is sensitive to the presence of small-scale features in the bathymetry: a 200 m high bump at 700 m depth can result in the generation of many mode-two waves and of higher mode waves. Sensitivity to the stratification is considered by using three stratifications based on summer observations. They primarily differ in the depth of the thermocline. The generation of mode-two waves and the behaviour of the waves in shallow water is sensitive to this depth. Rotation affects the shoaling waves by reducing the amplitude of the leading waves via the radiation of long trailing inertiagravity waves. The nonlinear-dispersive evolution of these inertia-gravity waves results in the formation of secondary mode-one wave packets.


Introduction
In the South China Sea (SCS) large internal solitary-like waves (ISWs) are frequently observed remotely via synthetic aperture radar (SAR) and via in situ observations (Liu et al., 1998;Ramp et al., 2004;Klymak et al., 2006;Farmer et al., 2009;Li and Farmer, 2011).Analysis of these observations indicate that the ISWs are generated by barotropic tidal motion over the large sills in Luzon Strait and on the slope region of the northern boundary of the SCS (Liu et al., 1998;Hsu and Liu, 2000;Alford et al., 2010).Depressions generated in Luzon Strait propagate westwards across the deep South China Sea basin, where depths can exceed 4000 m, to the Asian Seas International Acoustic Experiment (ASI-AEX) experimental site on the Chinese continental shelf.Nonlinear effects steepen these depressions until they disintegrate into ISWs of large amplitude, through frequency and amplitude dispersion.Zhao and Alford (2006) related the ISWs observed near Dongsha Island to tidal currents passing through Luzon Strait, finding that westward flow through the strait was responsible for the generation of the ISW packets; however, more recent studies have shown that eastward currents are responsible for their generation (Buijsman et al., 2010).Several numerical studies of ISW generation by tidal flow over the ridges in Luzon Strait have been conducted (Niwa andHibiya, 2001, 2004;Warn-Varnas et al., 2010;Buijsman et al., 2010).Li and Farmer (2011) reported deep basin ISW amplitudes as large as 150 m with typical values being around 50 m.In situ measurements in the deep basin by Klymak et al. (2006) show ISWs with amplitudes as large as 170 m and phase speeds of 2.9 m s −1 .During the joint Variations Around the Northern South China Sea (VANS) and Windy Island Soliton Experiment (WISE) (Yang et al., 2009), ISWs with amplitudes of over 190 m and phase speeds over 3.0 m s −1 were observed (Ramp et al., 2010).These measured amplitudes and phase speeds are among the largest observed anywhere, surpassing observed amplitudes of 90 m and phase speeds of 1.8 m s −1 in the Andaman and Sulu seas (Osborne and Burch, 1980;Apel et al., 1985).Their large amplitudes have made these waves the subject of many observational and numerical investigations.The dynamics of ISWs in the northern SCS was recently reviewed by Guo and Chen (2014).
The structure of ISW wave trains in the SCS is variable with large-amplitude wave trains separated by smallamplitude wave trains (Ramp et al., 2004;Warn-Varnas et al., 2010;Buijsman et al., 2010;Vlasenko et al., 2012).Ramp et al. (2004) separated the wave packets observed at the ASI-AEX site into a waves and b waves.The a waves arrived at regular 24 h intervals and were rank ordered, with the largest wave leading the wave packet.The weaker b waves were separated by approximately 25 h and were more irregular, with the largest wave usually in the middle of the packet.The formation of these wave packets is tied to the generation of the waves in Luzon Strait (Alford et al., 2010;Buijsman et al., 2010).
During the ASIAEX experiment, intensive measurements of shoaling ISW trains over the continental shelf off China were undertaken (Lynch et al., 2004;Duda et al., 2004;Orr and Mignerey, 2003;Ramp et al., 2004).Ramp et al. (2004) reported on observations from a series of moorings spanning depths of 350-72 m illustrating the evolution of these waves as they shoaled.At 350 m depth the waves were fairly symmetric with amplitudes ranging from 29 to 142 m (based on the displacement of the 24 • isotherm).By the time they reached a depth of 200 m they had usually deformed significantly with a gently sloping front and a much steeper rear.Orr and Mignerey (2003) tracked a shoaling ISW train with a towed conductivity, temperature, and depth sensor (CTD), observing the evolution from waves of depression to waves of elevation.Mode-two solitary-like waves were observed during the ASIAEX experiment (Orr and Mignerey, 2003) but not reported on in detail.Observations of a mode-two ISW made during a pilot study for the ASIAEX experiment were reported by Yang et al. (2004).Many mode-two solitary-like waves were observed at 350 m depth in a nearby region during the joint VANS/WISE experiments (Yang et al., 2009(Yang et al., , 2010)).Mode-two waves observed in the summer appear to have been generated by shoaling mode-one ISWs.Liu et al. (1998) demonstrated, with a two-layer Gardner equation model that includes cubic nonlinearity, that ISW depressions are transformed into a train of ISWs of elevation as they propagate over sloping topography into shallower water where the upper layer is thinner than the lower layer, as did Orr and Mignerey (2003), who based their simulations on observations from the ASIAEX site.Recently, Grimshaw et al. (2014) used the Ostrovsky equation to investigate the effects of rotation on the evolution of shoaling solitary waves.They too included a cubic nonlinear term.Their simulations were based on bathymetry and continuous stratifications observed in the South China Sea.They found that a significant consequence of rotation was the formation of a secondary trailing wave packet associated with nonlinear steepening of the radiated inertia-gravity waves.They also conducted fully nonlinear simulations using the MITgcm (MIT General Circulation Model).Vlasenko and Stashchuk (2007) studied threedimensional shoaling in the Andaman Sea with a fully nonlinear nonhydrostatic numerical model using a continuous stratification.Results showed the transformation of ISWs of depression to waves of elevation as well as refraction effects.
The breaking of solitary waves travelling up slope has been addressed by several researchers.Helfrich and Melville (1986) and Helfrich (1992) considered the breaking criteria for a two-layered system.Vlasenko and Hutter (2002) performed simulations, in a continuously stratified fluid, of breaking solitary waves on slope-shelf topography and determined a parametrisation for the location of wave breaking for stratifications and bathymetry based on observations in the Andaman and Sulu seas.
The large internal solitary waves that shoal onto the Chinese continental shelf are highly energetic features that have implications for biological productivity and sediment transport (Wang et al., 2007;Reeder et al., 2011).St. Laurent (2008) concluded that these waves drive one of the most dissipative coastal regions of the world's oceans.Thus, it is important to understand their shoaling dynamics.Observations have shown that shoaling ISWs exhibit significant structural changes (Lynch et al., 2004;Duda et al., 2004;Orr and Mignerey, 2003;Ramp et al., 2004).Lien et al. (2012Lien et al. ( , 2014) ) observed waves in the vicinity of the ASIAEX experiments with recirculating cores.The generation mechanism of these waves is unclear.
In this paper we seek to better understand the evolution of shoaling waves in the South China Sea via numerical simulations.We investigate the sensitivity of the predicted shoaling behaviour to the numerical resolutions employed, finding that higher resolutions than previously used are necessary to accurately predict wave forms in shallow water.We also add to the limited understanding of the evolution of wave amplitude and currents during shoaling by determining the depths at which maximum amplitudes and currents are attained and via comparisons with adiabatic theory.The generation of mode-two waves is also clarified.
While ISWs often appear as packets, we consider the evolution of a single ISW for simplicity.In this first study we also ignore the important effects of background tidal currents which will modify the shoaling behaviour by introducing time varying currents and stratification.Instead, we focus on exploring the sensitivity of the evolution of a single shoaling ISW for a range of wave amplitudes to the underlying bathymetry, small changes in stratification and the effects of rotation.These simulations are based on observations from the ASIAEX experimental site (Orr and Mignerey, 2003;Duda et al., 2004;Ramp et al., 2004).
The model, bathymetries and stratifications used in this study are described in Sect. 2. Results are presented in Sect. 3 where sensitivity to the bathymetry, stratification, rotational effects and the potential implications of viscosity are considered.We also present results from resolution tests and consider the adiabaticity of the shoaling waves.Comparisons with observations from the ASIAEX and conclusions are presented in Sect. 4.

Numerical model
The two-dimensional fully nonhydrostatic internal gravity wave model (Lamb, 1994(Lamb, , 2007) ) is used in this study of shoaling ISWs.Flow in the along-shelf y direction is included but the flow is independent of y.The model uses the rigid-lid, traditional f plane, Boussinesq and incompressible flow approximations.The governing equations are where z is the vertical coordinate, (x, y) are the horizontal coordinates directed across and along the topography, u = (u, v, w) is the velocity field with components in the (x, y, z) directions, k is the unit vector in the z direction, ρ is the density, ρ 0 is the reference density, and p is the pressure.∇ is the 3-D vector gradient operator; however, all fields are independent of y.The subscript t denotes differentiation with respect to time.The gravitational acceleration is g = 9.81 m s −2 and for simulations including rotational effects the Coriolis parameter is f = 5.33 × 10 −5 s −1 corresponding to a latitude of 21.5 • N. ν and κ are a spatially varying viscosity and diffusivity; however, for most of our simulations these are set to zero.The model uses flux limiting which acts to control the overturns that occur in the simulations.In the ocean strong

Bathymetry
A total of 25 sample shelf-break bathymetries in the area of ASIAEX (see red square in Fig. 1), as extracted from the digital bathymetry 2 min resolution (DB2) database, are shown in Fig. 2. They contain small-scale variability indicative of the presence of small-scale, three-dimensional bathymetric features.The solid and dotted black curves in Fig. 2a show the base bathymetries used in our simulations.They are compared with two measured bathymetries in Fig. 2b and c and are referred to as h 0 (transect 0) and h 15 (transect 15) bathymetries, respectively.These are representative of transects with the smallest and largest slopes.Bathymetric slopes are indicated in the figure.Bathymetry h 0 has a bump of about 200 m amplitude at 2100 m depth.Bathymetry h 15 has a 10 km long, 200 m high bump at 700 m depth.Some simulations were done with this feature removed (see Fig. 2c).
In the simulations the grid is arranged so that the shelf starts at approximately x = 80 km.

Stratification and model initialization
Three stratifications have been used (Fig. 3).Four stratifications are shown in the first two panels: our base stratification ρ b (z), two fits to observed stratifications measured at the ASIAEX site (Orr and Mignerey, 2003), ρ 1 and ρ 2 , and ρ w , which is the density west of Luzon Strait used in Warn-Varnas et al. (2010).Profiles ρ 1 (z) and ρ 2 (z) are compared with the corresponding observed profiles in Fig. 3b.These profiles were chosen because they have the thermocline at the lower and upper range of observed thermocline depths.
We initialised our simulations with internal solitary wave solutions of the governing equations (with f = 0) obtained by solving the Dubreil-Jacotin-Long (DJL) equation (Lamb, 2002) which provides an initial ISW with a prescribed available potential energy (APE).Here and throughout the APE is for waves in an infinite domain computed using the background density as the reference density (Lamb and Nguyen, 2009).In simulations that include rotational effects the initial waves undergo a continual adjustment which results in a continual loss of energy from the leading ISW (Helfrich, 2007;Grimshaw et al., 2014).Because of this the shoaling behaviour will depend on the initial location of the ISW.We do not consider this.The initial waves are placed in the deep water.Because of the different lengths of the shelf slopes the initial location depends on the bathymetry: for bathymetries h 0 and h 15 the waves are started at x = −280 and −200 km, respectively.
We have used waves with several different amplitudes based on observed amplitudes (Ramp et al., 2004).The energies and wave amplitudes for several different initial waves are provided in Table 1.

Weakly nonlinear theory and conjugate flow amplitudes
Weakly nonlinear theory is often used to approximate observed internal solitary waves and it provides some insight into the expected behaviour of shoaling waves.Including cubic nonlinearity and ignoring damping and rotational effects, internal solitary waves in a fluid of constant depth can be modelled with the Gardner (or extended KdV) equation Table 1.Properties of Initial ISWs in 3000 m depth.ρ is the density stratification used, E a and E k are the available potential energy and kinetic energy of the wave, a is the wave amplitude (maximum isopycnal displacement), c the wave propagation speed and U max is the maximum horizontal current in the wave.
The KdV equation (α 1 = 0) predicts solitary waves of depression when α < 0 and of elevation when α > 0. Stratifications ρ b and ρ 1 have critical points, locations where α changes sign, at depths of about 91.8 and 120.0 m, respectively: in greater/shallower depths ISWs are waves of depression/elevation. Thus, for these two stratifications shoaling waves pass through a critical point on their way to the shelf which has significant implications for their evolution.Stratification ρ 2 , with its higher pycnocline, does not have a critical point for the chosen shelf depth of 80 m.When α 1 < 0, solitary wave solutions of the Gardner equation have a limiting amplitude of −α/α 1 .The solitary wave solutions become flat crested as the energy in the wave increases.For α 1 > 0, solitary waves of depression and elevation with unbounded amplitudes exist (Grimshaw et al., 2004).Using α 1 to predict the existence of solitary waves with both polarities is problematic because α 1 is not uniquely determined (Lamb and Yan, 1996).Choosing α 1 so that the second-order vertical structure function for the isopycnal displacement is zero at the depth where the leading-order eigenfunction has its maximum is a common choice (Grimshaw et al., 2002).This choice for selecting α 1 predicts that solitary waves of either polarity (i.e.elevation or depression) exist in depths greater than 235 m using our base stratification.In shallower water it predicts that only waves of one polarity exist.Fully nonlinear ISWs generally have limiting amplitudes.As this amplitude is approached they broaden and become horizontally uniform in their centre.The flow state in the centre of these flat-crested waves is called a conjugate flow.We computed conjugate flow solutions for depths between 50 and 400 m which provide the limiting asymptotic amplitude, propagation speeds and maximal currents for solitary waves as the wave energy goes to infinity (Lamb and Wan, 1998).The amplitudes of these solutions for the base stratification ρ b varied approximately linearly from −160 to 20 m as the depth decreased from 400 to 60 m.Only single solutions were found with the polarity predicted by the KdV equation.

Results
Table 2 provides information on the simulations that have been undertaken using a deep water depth of 3000 m.We begin by showing results for two simulations using an initial wave of amplitude 45.4 m (APE 100 MJ m −1 ), the base stratification ρ b and the two bathymetries (Cases 2 and 3).Following this we explore sensitivity to resolution before discussing more fully the complete set of simulations.
Figure 4 shows results from Case 2. Figure 4a shows the full continental slope and the initial wave, barely visible at x = −280 km.After 25 h the leading wave is at x = −64 km, where the water depth is about 650 m (Fig. 4b).A second solitary wave is visible approximately 15 km behind.At t = 42 h (Fig. 4c) the leading wave, now 1.5 km in length, has reached a depth of 250 m at x = 37 km.Three solitary waves are now visible.By t = 50 h (Fig. 4d) the leading wave has significantly deformed.It is now about 3 km long, spanning depths of 130-140 m.The slope of the wave front is less than it was at t = 42 h and in the centre of the wave the thermocline is almost parallel to the bottom (between 65.5 and 67.5 km).This is a well-know feature of shoaling waves when the thermocline is close to the bottom (Vlasenko and Hutter, 2002;Orr and Mignerey, 2003;Lamb and Nguyen, 2009).The rear of the wave is much steeper and behind the leading depression the thermocline has been raised above its equilibrium position (x = 66 km).This we refer to as the wave pedestal.Similarly raised pycnoclines can occur behind large internal lee waves generated by near-critical flow over a bump (Stastna and Peltier, 2005).There are now three distinct waves of depression trailing the leading wave.As the leading wave continues to shoal and passes beyond the critical point at 91 m depth (Fig. 4e), the front part of it becomes very long with a very shallow slope.On the shelf, solitary waves of depression in the ambient stratification do not exist so the leading depression gradually fades away losing energy to the trailing waves.The steep back of the wave never over- Table 2. Shoaling cases starting at 3000 m depth.Cases run using bathymetry h 15 without the bump are indicated with an optional nb and cases run with rotation turned on are indicated by an optional r.Three sets of simulations were run with and without viscosity.These are indicated by the (ν).Cases numbers with an optional ν indicate cases run with viscosity.For these cases three values of ν were used, the value of ν being indicated by a number, e.g.νn, where values of n = 3, 4, or 5 are for K = 10 −n m 2 s −1 .Thus, Case 4nb was run with the bump removed and Case 4_r was the same as Case 4 but with rotation turned on.Case 9ν4 was run with viscosity using turns and the second solitary wave, which has now propagated into the pedestal behind the leading wave, has the form of a square wave with the centre of the wave parallel to the bottom separating a very steep front and back (Fig. 4e).This wave is approximately at the critical point of the background stratification but it is propagating on the wave pedestal trailing the leading wave and this is the environment it propagates in.This environment consists of a background-sheared flow with negative vorticity (u z < 0) and a modified stratification.By t = 66 h the waves are on the shelf.Breaking has commenced behind the third depression.(c, d) Profiles in solitary waves of depression (solid) computed on background field given by profiles extracted from simulation in the pedestal at x = 92 km between the two leading waves of depression.Solution is compared with profiles taken from the second depression at x = 91 km (dots).tends from the back of the leading depression at x = 80 km to possibly x ≈ 74 km).In the two wave depressions (profiles D1 and 2 at x = 94 and 91 km) the horizontal velocities are similar.In the pedestal the horizontal currents are stronger behind the second depression (profile E2) than they are ahead of it (E1).Furthermore, behind the second depression the current at the bottom has a bulge with the maximum current of about 0.5 m s −1 occurring about 5 m above the bottom.This is equal to the estimated propagation speed of the rear of the second part of the pedestal (x ≈ 88.5 km).The density profiles (Fig. 6b) show that behind the second depression fluid in the lower 22 m is denser than the fluid ahead of the wave, a consequence of advection of dense water onto the shelf.This denser fluid can also be seen in Fig. 5 where a thin layer of this denser fluid can be seen to extend along the bottom into the depression ahead.A small density overturn can be seen near the bottom of profile E2 (Fig. 6b).The pedestal has an amplitude of over 20 m which is well in excess of the conjugate flow amplitude of 7.7 m for the ambient stratification on the shelf.We calculated ISWs of depression for a background field given by conditions in the pedestal at x = 92 km by solving a version of the DJL equation that includes background currents (Stastna and Lamb, 2002;Lamb, 2003).The computed waves had a maximum amplitude of 36.6 m.At this amplitude waves are flat crested, similar to the second depression centred at 90.5 km in Fig. 5.The flow in the centre of the computed solitary wave is the conjugate flow corresponding to the background conditions.Figure 6c and d compares the vertical profiles of the horizontal velocity and density field of the conjugate flow in the centre of the computed ISW with the vertical profiles from the second depression in the simulation (from x = 91 km).The two sets of profiles are virtually identical except near the bottom showing that the second square-shaped wave of depression is a flat-crested solitary wave riding on the background flow provided by the wave pedestal.The waves on the shelf have not however achieved a steady state.The amplitude of the leading depression slowly decays (solitary waves of depression do not exist on the shelf) and as it does the amplitude of the trailing pedestal also decays as do the waves behind it.
Figure 7 shows results from a simulation using bathymetry h 15 with the same initial wave (Case 3), now at x = −200 km.The results are similar to those using bathymetry h 0 but there are some differences.In both cases the shoaling wave fissions into two large solitary waves trailed by a number of smaller waves.For bathymetry h 0 the trailing waves consist largely of a couple of smaller solitary waves (at x ≈ 55 km in Fig. 4d) whereas for the steeper bathymetry h 15 the trailing waves are smaller and there are many more of them (e.g.Fig. 7c).The shelf slope for bathymetry h 15 is shorter than that for bathymetry h 0 which is likely responsible for the two leading waves being closer together when bathymetry h 15 is used (500 m apart vs. 1500 m).In addition to the differences in the bathymetric slopes h 15 has a bump at 700 m depth.This bump modifies the wave field resulting in the generation of higher mode waves (Fig. 7b) as can be seen by the internal wave beam sloping up from the bump.A zoom in of the higher mode waves at t = 29 h (same time shown in Fig. 7c) is shown in Fig. 8.The two dashed white boxes contain mode-two waves (the left box also includes higher mode waves).Downwelling at the bottom behind the bump after the wave has passed over it can be seen between x ≈ −17 and −15 km (Fig. 7b).In some cases this can result in larger vertical excursions than in the leading wave.Breaking behind the shoaling waves commences earlier for the steeper bathymetry (compare Figs. 4e and 7e).Significant wave reflection also occurs in this case due to the steeper bathymetry.
Theoretically, energy is conserved as waves shoal if viscosity and diffusion is ignored; however, in the simulations the total energy (kinetic plus available potential energy) changes due to numerical error.It can increase slightly in deep water due to numerical diffusion thickening the pycnocline.For Case 2 the wave field at t = 42 h (see Fig. 4c) contains 96 % of the initial energy, with 90 % of the initial energy contained in the wave field on the gently sloping part of the shelf slope (x > 90 km).The leading two solitary waves contain 50 and 16 % of the initial energy.At t = 50 h the leading depression contains 42 % of the initial energy with 7 % of the initial energy now residing in the pedestal between the two leading solitary waves.At t = 58 h only 84 % of the initial energy remains.The leading two waves of depression contain 19.9 and 10.9 % of the initial wave energy while between them the pedestal has 10.5 % of the initial energy.
For Case 3 at t = 29 h (Figs.7c, 8) the energy in the full wave field is 98.8 % of the original energy.The leading wave has about 36 % of the original energy, vs. 50 % for Case 2. The second solitary wave is not cleanly separated from the trailing waves so a precise estimation of its energy is not possible; however, about 17 % of the wave energy lies between 31 and 35 km.Approximately 15 % of the wave energy is in reflected waves to the left of the bump (x < −15 km).
In summary, the steep slope between depths of 2250 and 750 m for bathymetry h 15 results in significant reflection of incident wave energy.The leading solitary wave in shallow water contains less energy than for the case using bathymetry h 0 which has a gentler slope.For bathymetry, h 15 two large leading solitary waves are formed in shallow water which are trailed by a train of small-amplitude mode-one and higher mode waves that are generated by the interaction of the shoaling wave with the bump.For bathymetry h 0 there are 3-4 solitary waves and the leading solitary wave has about 50 % more energy that that in the other case.No trailing modeone wave train or higher mode waves are evident.When the waves arrive on the shelf a large pedestal forms which is conjugate to the flow in the leading depression.These second depression waves are much larger than the largest possible solitary wave for the ambient background conditions.showing mode-two waves (in white dashed boxes) and an internal wave beam (above the bump) generated when the solitary wave passes over the bump.

Resolution
It is important to have high enough resolution to adequately resolve solitary waves as they shoal.If the resolution is too low unphysical solitary-like waves can be simulated due to a balance between nonlinearity and numerical dispersion (Hodges et al., 2006;Vitousek and Fringer, 2011).Vitousek and Fringer (2011) considered waves propagating in a 2000 m deep domain using a stratification based on observations in the South China Sea and found that resolutions coarser than about 250 m resulted in waves that were too wide.
In the first set of tests, using an initial wave of 45 m amplitude and bathymetry h 15 , the horizontal resolution was varied ( x = 250, 100, 50, 33 and 16 m).The maximum time steps were 12.5, 5.0, 2.5, 2.5 and 2.5 s, respectively, and a uniform vertical resolution of J = 200 was used.For the highest-resolution case the time step decreased to approximately 1.4 s by the time the wave reached the shelf.Doubling the time step to 5.0 s in the x = 50 m simulation resulted in negligible differences.After travelling over 120 km to a depth of 2500 m the amplitudes in all five cases varied by about 2 %.At 800 m depth, just before the bump, the solitary wave in the x = 250 m simulation was about 5 % smaller than that in the highest-resolution case.After this stage the solutions diverged more rapidly.
Figure 9 compares the surface currents u(x, 0, t) when the leading wave is in water depths of 600, 395, 180 and 85 m.At 600 m depth (Fig. 9a) the leading wave has begun fissioning into several waves.The wave in the 250 m resolution case is significantly smaller than those in the other simulations while the wave in the 100 m resolution cases is almost indistinguishable from the high-resolution case.At 395 m depth (Fig. 9b) the amplitude of the leading wave in the 250 m resolution case is grossly underestimated and has the form of a small undular bore.In contrast, in the higher-resolution simulations two large leading solitary waves have emerged trailed by a train of smaller waves.The amplitude of the leading wave in the 100 m resolution simulation is now about 87 % of the highest-resolution case.At 180 m depth (Fig. 9c) this fraction has reduced to 77 %.The waves in the 50 m simulation are now slightly smaller than in the 33 m resolution case.In the 50 and 33 m resolution simulations the thermocline has been raised above its rest height behind the leading wave as indicated by the negative surface currents.This fea-ture is missing at lower resolutions.At later times the 50 and 33 m resolution simulations diverge as well, becoming quite different at 85 m depth (Fig. 9d).The back of the leading depression becomes very steep.It is steeper in the 33 m resolution case and is resolved by a 4-5 grid at both of these resolutions.In the 16 m resolution case the jump is better resolved (7-8 grid points).Wave breaking has commenced behind the two leading depressions just before this time.The 16 m resolution case shows much finer details and some weak overturns along the steep rear of the first square-wave depression.At this stage it is unlikely that the 16 m resolution case has converged.At depths greater than about 140 m the differences between the 50, 33 and 16 m resolution cases are not very significant.
Results from a simulation with x = 33 m and double the vertical resolution (J = 400) are indistinguishable at the scale of Fig. 9d.
For a larger wave (amplitude 115 m), with approximately half the width of the 45 m wave requiring higher horizontal resolutions, the dependence on resolution was qualitatively similar to first test set with the results diverging more rapidly.After travelling 8.5 km to a depth of 2800 the wave in the 250 m resolution simulation was already over 5 % smaller than that in the 33 m simulation, showing that at a depth of 3000 m the resolution needs to be significantly better than 250 m to accurately simulate waves of this amplitude.
While doubling the vertical resolution did not result in any significant changes in the leading waves, it did reduce noise that developed in the deep water behind the waves.When a resolution of 200 grid points was used weak (< 0.01 m s −1 ) near-surface currents (across the pycnocline at a depth of 50 m) developed over the steepest parts of the bathymetry, indicative of the pressure-gradient error that is intrinsic to sigma-coordinate models.With 400 grid points these currents were reduced by 50 %.Simulations with a resolution of 300 and 400 vertical grid points were similar and were very expensive.Because the leading waves were unaffected by increasing the vertical resolution from 200 to 400 we have for the most part used a vertically varying grid with 200 grid points in the vertical.The grid has almost uniformly spaced grid cells for depths greater than 700 m and less than 300 m in deep and shallow water with a transition region centred at 500 m depth.In deep water the resolution was between 21.1 and 21.4 m at depths greater than 700 m and decreases from 6 m in at 350 m depth to 5.35 m at the surface.The vertically varying grid better resolves the thermocline in deep water and improves energy conservation by reducing thickening of the pycnocline by numerical diffusion.In water with depths shallower than 350 m the resolution is almost uniform, with dz = 0.4 m on the shelf.This grid gave virtually identical results to the vertically uniform grid using 400 vertical grid points.

Adiabaticity of shoaling waves
Solitary waves are said to shoal adiabatically if they retain the form of a single solitary wave as they shoal, i.e. no fissioning or reflection occurs (Grimshaw et al., 2004;Vlasenko et al., 2005).This is an idealisation that can approximately occur if the water depth changes so slowly that the shoaling wave can adapt to its new depth without shedding significant energy.Adiabatically shoaling waves conserve energy to leading-order in amplitude (in the absence of rotational effects); however, they are in general trailed by a long, small-amplitude shelf which extends from the back of the solitary wave to the point reached by a wave, travelling at the long-wave propagation speed, generated at the point when the shoaling solitary wave first encounters a depth change (Grimshaw et al., 2004).Because energy is second-order in amplitude the energy in the trailing shelf is negligible, however the mass in the shelf, being first-order in amplitude, can be significant.
Here we investigate the adiabaticity of the shoaling waves by comparing the leading wave to ISW solutions of the DJL equation having the same total energy as the initial deep water wave.Both bathymetries are used.We also launch waves from different water depths using bathymetry h 0 levelled off at deep water depths of 3000, 1500 and 1000 m (see Fig. 2b).The base stratification was used except where noted.
ISW solutions of the DJL equation were computed over a range of water depths for 11 values of the APE: 25 MJ m −1 and 50-500 MJ m −1 in increments of 50 MJ m −1 .Figure 10a shows the total wave energy E (KE plus APE) as a function of depth for these 11 values.The ratio of kinetic to available potential energy (Fig. 10b) is always greater than one for solitary waves (Lamb and Nguyen, 2009).At 3000 m depth this ratio varies from slightly more than 1 to a maximum of 1.2 for the largest wave (amplitude 137 m) for which values are plotted.For larger waves (not shown) this ratio continues to increase, reaching 1.3 for an APE of 800 MJ m −1 and 1.4 for an APE of 2 GJ m −1 .These waves have amplitudes of 170 and 322 m and minimum Richardson numbers of just below 0.25 and 0.11, respectively.For the larger wave the ratio of the maximum current to the propagation speed is max(u)/c = 0.98.Waves with max(u)/c > 1 have cores of recirculating fluid.
As the water depth decreases from 3000 m with the APE fixed, the total wave energy increases due to the increasing KE / APE ratio, reaching a peak value at depths be- Coloured curves are maximum absolute amplitude and surface current and minimum bottom currents of shoaling waves as function of water depth.Blue -cases using bathymetry h 0 .Redcases using bathymetry h 15 .Black dash-dot -higher-resolution case using bathymetry h 0 .This curve overlies the lower-resolution result (blue curve, initial amplitude 115 m).Both are overlain by the results using bathymetry h 15 for depths greater than about 2600 m.
tween about 400 and 1200 m for small and large waves (KE / APE ≈ 1.16 and 1.3, respectively, Fig. 10b).As the depth decreases further the total energy rapidly decreases to twice the APE.This near equipartition of energy suggests waves of smaller amplitude.This is indeed borne out in Fig. 11 which shows the wave amplitude (maximum isopycnal displacement), maximum surface current and minimum bottom current (negative for rightward propagating waves of depression) as a function of depth for constant total energies E of 50 MJ m −1 and 100-1000 MJ m −1 in increments of 100 MJ m −1 .These values were obtained by interpolating results from the DJL solutions, for which the APE is prescribed, and will be referred to as the adiabatic curves: these are the curves a shoaling wave would follow if the wave preserved its total energy and maintained the form of a solitary wave for its depth.As the water depth decreases from 3000 m the wave amplitude increases significantly, reaching a maximum value in depths between 300 and 600 m.The largest wave (total energy of 1 GJ m −1 ) increases in amplitude by a factor of 1.4 from 129 to 183 m.For smaller waves the amplitudes increase by larger factors as the depth decreases (e.g.3.5, 2.9 and 2.3 for the three smallest waves).Ampli-tudes thereafter decrease rapidly to close to 60 m in water of 200 m depth.This is an indication of the conjugate flow limit being attained: the waves all have the same maximum amplitude and as the wave energy increases the waves get longer with negligible increases in amplitude.
The adiabatic curves are similar for the other two stratifications.Table 3 gives the initial and maximum amplitudes and the depth at which the maximum amplitude occurs for ISW solutions of the DJL equation for waves of fixed total energy for all three stratifications.The increase in amplitude of the waves is similar for stratifications ρ b and ρ 1 while waves using stratification ρ 2 show a larger amplitude increase, particularly for small waves.
A similar pattern is seen for the maximum surface currents except that the maximum currents occur in greater depths than do the maximum wave amplitudes and the range of depths at which near-maximum currents occur is much wider than those for near-maximum amplitudes.The bottom currents reach their minimum values in much shallower water and the increase in magnitude over their deep water values is much larger than the increases in wave amplitude and surface currents.
Figure 11 also shows the maximum wave amplitude, maximum surface current and minimum currents at the bottom from several simulations.For these simulations f = 0, as the adiabatic curves are only appropriate in the absence of rotation.Five simulations were done using bathymetry h 0 .Four of these used a horizontal resolution of about 33 m with 200 uniformly spaced grid points in the vertical (blue curves).These waves had APEs of 50, 100, 250 and 400 MJ m −1 .The fifth case (black dash-dot), also with an APE of 400 MJ m −1 , used double the vertical resolution.It largely overlies the corresponding lower-resolution case.Results from simulations using bathymetry h 15 are also shown (red curves).The initial APEs for waves launched from 3000 m depth were 100, 150, 250 and 400 MJ m −1 .The wave started at a depth of 1500 m had an initial APE of 375 MJ m −1 , chosen so the total wave energy is similar to that of the wave with an APE of 400 MJ m −1 started at 3000 m depth.Waves launched from 1000 m depth had APEs of 50, 60, 70 and 100 MJ m −1 .
Figure 11a and b show that the wave amplitudes and maximum surface currents tend to increase more slowly than the adiabatic curves in deep water and follow them more closely in shallow water.The depths at which they begin to follow the adiabatic curves more closely depend on the bathymetry and wave amplitude.
Consider the results obtained using bathymetry h 0 (blue curves) which has a bump at a depth of about 2000 m.When the wave launched at a depth of 3000 m with an APE of 400 MJ m −1 (initial amplitude 115 m, E = 847 MJ m −1 ) first reaches the crest of the bump (depth 1900 m) its amplitude is equal to that of a solitary wave with E equal to 750 MJ m −1 , while the maximum surface current corresponds to that for a wave with a total energy of 670 MJ m −1 .By the time the Table 3. Maximum increase in wave amplitude of ISW solutions of the DJL equation.First column is total energy (available potential and kinetic energy).a 3000 is the wave amplitude in a water depth of 3000 m. maxa is the maximum wave amplitude for the given energy, d max is the water depth where the maximum occurs.The three, separate values are values for densities ρ b , ρ 1 and ρ 2 .Given values are based on DJL solutions calculated at depth intervals of 50 m (see Fig. 11).wave has passed over the bump and shoaled to a depth of 1900 m for the second time the wave amplitude and, to a much greater degree, the maximum surface current have increased to the values corresponding to an ISW with energies closer to that of the original wave.The increase in the maximum surface current is particularly striking.This is an indication the wave is not shoaling adiabatically, i.e. it has not fully adjusted to the change in water depth.Results for this wave track an adiabatic curve for depths shallower than about 1600 m. Results for the other wave amplitudes are similar.In contrast to the adiabatic curves, the two smallest waves launched from a depth of 3000 m show little increase in the wave amplitude and maximum surface current until they have reached depths of about 750 m, where the bottom slope decreases from 0.03 to 0.004.At this point the amplitude and maximum surface current of the wave with an initial APE of 100 MJ m −1 (E = 204 MJ m −1 ) is well below those for the DJL waves with an initial APE of 50 MJ m −1 .After this the amplitudes and surface currents rise rapidly and fall tracking the adiabatic curves in depths shallower than about 700 m.This suggests significantly non-adiabatic shoaling over the steep slope and nearly adiabatic shoaling over the gentle slope between depths of 750 and 80 m.Similar behaviour is observed for the waves launched from a depth of 1000 m.

E
Results obtained using bathymetry h 15 show a similar pattern.Because this bathymetry is much steeper than h 0 for depths between 500 and 2300 m, the results with this bathymetry (red curves) do not follow the adiabatic curves until the waves reach much shallower depths of between 500 and 600 m.The effects of the bump at a depth of 600 m are clearly visible.The large loops in the wave amplitudes, particularly evident in the smaller waves, are a consequence of the largest vertical displacements occurring behind the bump after the solitary wave passes over it, rather than in the solitary wave itself, for a short period of time.
In water shallower than about 400 m the large shoaling waves have amplitudes that exceed the conjugate flow limit by as much at 50 %: at 350 and 200 m depths the largest shoaling waves have amplitudes of 142 and 95 m while the corresponding conjugate flow amplitudes are 137 and 61 m.Surface currents do not exceed the conjugate flow limits until depths of 250 m are reached.
Figure 11c shows the evolution of the strongest bottom currents.They greatly exceed those for the DJL solutions in water shallower than about 300 m for the large waves and about 200 m for the smaller waves.This is a consequence of the fluid below the thermocline being squeezed out from beneath the waves as the waves shoal which, in shallow water, results in enhanced bottom currents.This has important implications for the occurrence of instabilities, mixing and sediment resuspension beneath the shoaling wave.Also noteworthy is that for the large waves the bottom current increases relative to the adiabatic curves as the waves shoal to about 500 m depth.For the smaller waves they decrease as do the amplitudes and surface currents.For large waves in shallow water the bottom currents greatly exceed the surface currents.
In summary, the amplitudes and maximum surface and bottom currents of the shoaling waves qualitatively follow the behaviour expected of a wave that is shoaling adiabatically; however, there are some significant differences, a consequence of the timescale of depth changes experienced by the shoaling wave being shorter than the time required for the waves to adjust to their new depth.The amplitude and surface currents increase more slowly in deep water, particularly over the steep slope between depths of 2250 and 750 m for bathymetry h 15 .Part of the difference can be ascribed to the waves having not fully adjusted to the changing depth rather to fissioning as illustrated by the behaviour as the waves pass over the bump at a depth of 2000 m for bathymetry h 0 .For depths between about 1000 and 400 m, the properties of the shoaling waves track the adiabatic curves very well.In shallower water the waves can have amplitudes more than 50 % larger than any exact solitary wave at that depth while currents at the bottom can be twice as large.

Sensitivity to initial water depth
For reasons of numerical efficiency it is tempting to truncate the bathymetry and start waves from depths of, say, 1000 m.Here we briefly explore the implications of doing so.
Figure 12 compares the horizontal surface current and σ θ ≡ ρ − 1000 = 23.2kg m −3 isopycnal at different times for cases with initial APEs of 60 and 100 MJ m −1 launched from depths of 1000 and 3000 m, respectively, using bathymetry h 0 .The latter case is Case 2 illustrated in Fig. 4. The chosen isopycnal is at the depth of maximum buoyancy frequency in the undisturbed stratification.It is not the isopycnal which undergoes maximum vertical displacement in the wave which varies with wave amplitude and as the wave shoals.
These two waves are chosen because in shallow water the leading waves have nearly identical amplitude.At a depth of 750 m (Fig. 12a, b), where the shelf slope decreases from 0.03 to 0.004, the waves are about 15 and 150 km beyond the bottom of their shelf slopes.Both waves are asymmetric, with the wave launched from deeper water being much more asymmetric.By the time the wave has reached a depth of 560 m the wave launched from the deeper/shallower water has fissioned into three/two solitary waves (not shown).At a depth of 250 m (Fig. 12c, d) the multiple solitary waves produced via fissioning are well separated.The third wave in Case 2 has a similar amplitude and location to that of the second wave for the case started at 1000 m depth.The pedestal behind the leading wave starts to form.By the time the waves have reached depths of about 120 m (Fig. 12e, f) the third wave in the deep water case has now split into two waves.
As the waves pass the critical point and reach their final depth of 80 m (Fig. 12g, h) the leading wave continues to broaden and the currents decay in amplitude.The second wave in the deep water case advances into the wave shelf trailing the leading wave (the region of negative surface current), with the result that the wave that forms behind the leading depression is quite different for the two cases.
Truncating the deep water depth to 1000 m significantly modifies the fissioning processes resulting in fewer waves, reducing the relative amplitude of the second solitary wave and increasing the distance between the two leading waves.

Effects of small-scale bumps
The deep bump on bathymetry h 0 does not significantly affect the evolution of the shoaling waves.The much shallower bump on bathymetry h 15 has a significant impact.Figure 13 compares the wave fields at t = 29 and 34 h for simulations with and without the bump for an initial wave amplitude of 59 m (APE 150 MJ m −1 , Case 4).With the bump the leading three waves at t = 29 h are at x = 41.0,35.9 and 31.8 km.With the bump removed the waves are slightly ahead, at x = 42.8,37.3 and 33.7 km, which we attribute to the reduction in propagation speed as the wave passes over the bump.The bump also results in a reduction in wave amplitude, with the leading three waves containing 60 and 71 % of the initial wave energy for cases Cases 4 and 4nb.
The most striking difference between the two simulations is the wave field behind the leading three waves.For Case 4 a small-amplitude mode-one wave train trails the leading three waves.Following are higher mode waves, e.g. at time t = 29 (Fig. 13a), a concave mode-two wave is centred at x = 16 km and a convex mode-two wave (and higher mode waves) lies between x = 0 and 6 km with higher mode waves further behind and a wave beam sloping up to the right above the bump.At the later times mode-two waves in both simulations are clearly apparent, however with the bump there is a modetwo concave wave (x = 30 km) followed by a train of modetwo waves (between x = 16 and 20 km) whereas without the bump there is a single convex mode-two wave at x = 20 km.The bump also results in reflected waves containing about 6 % of the initial wave energy (less than the 15 % reflected for Case 3 discussed in Sect.3.1).In Case 4nb energy in the reflected waves is insignificant.

Sensitivity to initial wave amplitude
Figure 14 compares the shoaling behaviour for three initial waves (amplitudes 45, 83 and 115 m) using the base stratification ρ b and bathymetry h 15 (Cases 3, 6, and 7).Times for comparison are chosen so that the leading waves are at approximately the same location.For the smallest initial wave there are two large solitary waves at t = 35 h (Fig. 14b) whereas for the two larger waves there are three.The intermediate wave has a number of smaller solitary waves trailing the leading large waves that are absent in the two cases with smaller and larger initial waves.At the early times shown (panels a, c, e) the leading wave contains 36, 53 and 75 % of the initial wave energy showing a striking increase with the initial wave amplitude.The higher mode waves are similar for these three cases (not shown).

Effects of rotation
The effects of rotation on the evolution of the internal tide and on ISWs has been considered by many authors (e.g.Helfrich, 2007;Helfrich and Grimshaw, 2008;Grimshaw et al., 2014).Rotation makes long waves dispersive and as a consequence ISWs radiate long inertia-gravity (Poincaré) waves behind them and gradually decrease in amplitude.In some regions of parameter space the long radiated waves steepen and form new ISW packets (Helfrich, 2007).This process takes place on the inertial time period which in the South China Sea is about 32 h.This is comparable to the shoaling times in our simulations so it is clear that rotational effects will significantly affect waves launched from the deep water by the time they reach the shelf.
The effects of rotation on the shoaling of ISWs was investigated using bathymetry h 15 with the initial waves at x = −200 km.Starting the waves further from the shelf would increase the effects of rotation by reducing the amplitude of the wave at the bottom of the shelf slope and by introducing a long inertia-gravity wave behind the leading solitary wave.As a proxy for the first effect we can simply consider waves with different initial amplitudes.
Figure 15 compares the σ θ = 23.2kg m −3 isopycnals at different times for Cases 3(r) and 5(r).Figure 15a-c show results for an initial wave amplitude of 45.4 m (Cases 3 and 3r).At t = 23 h (Fig. 15a) the leading wave in the non-rotational case is just past the bump at 610 m depth.The isopycnal undergoes a maximum downward displacement of about 46 m.At this location the σ θ = 23.2kg m −3 isopycnal is not the isopycnal undergoing maximum displacement (the amplitude of the leading wave is about 53 m).In the rotational case the leading wave trails by almost 2 km and the isopycnal undergoes a maximum downward displacement of 26 m, almost half that in the non-rotating case.At this point the wave has been travelling for about three-quarters of an inertial period and has propagated 200 km.It is no surprise then that in the rotational case the leading wave has been significantly modified.The long small-amplitude wave of elevation between x = −25 and −12 km is slightly larger in the rotational case.The currents in the rotational case are much larger than those in the non-rotating case (not shown).After 35 h (Fig. 15b) the leading waves have reached a depth of about 200 m.At this depth the wave amplitude is strongly controlled by the water depth and the waves in the non-rotating and rotating cases have similar amplitudes, with the wave in the non-rotating case being much wider.In the non-rotating case three ISWs are present.There are only two in the rotating case.The wave of elevation between x = 30 and 50 km is larger in the rotating case (surface currents are twice as strong) and in the rotating case it subsequently steepens and forms a secondary packet of short waves (between x = 67 and 70 km at t = 45 h; Fig. 15c) which is not present in the non-rotating case.Such wave packets were reported by Grimshaw et al. (2014).
Figure 15d and e show results for larger initial waves (71.6 m, Cases 5 and 5r) and similar trends are observed.The secondary wave packet does not form in this case because its formation is disrupted by wave breaking behind the leading mode-one solitary waves.Starting the initial wave further away would give more time for the secondary wave packets to form, hence potentially making them more prominent.
Figure 16 shows the amplitude of the leading wave as a function of x as the waves shoal.Results for five pairs of simulations (with and without rotation) using different initial wave amplitudes are shown (Cases 3-7).When f = 0 the waves slowly increase in amplitude until the bump (x = −20 km) is reached.At this stage the wave fissions and the amplitudes of the three smallest waves momentarily decrease.All waves then rapidly increase in amplitude.The largest wave reaches its peak amplitude just before the crest of the bump is reached, then maintains its amplitude before starting to decrease in amplitude rapidly at x = 40 km (420 m depth).The other waves reach their maximum amplitude between x = 40 and 50 km (depths of 400-300 m) before rapidly decreasing until the shelf is reached at x = 80 km with smaller waves reaching their maximum amplitude in shallower water.
When rotational effects are included the wave amplitude immediately starts to decrease, approximately linearly with x.After passing over the bump at x = −20 km the effects of rotation on the leading wave no longer seem to be prominent, with the amplitude varying similarly to those in the non-rotating cases.

Sensitivity to stratification
We next consider the sensitivity of wave shoaling to variations in stratification similar to those observed over the course of the ASIAEX field program.The major difference in our stratifications is the depth of the thermocline which also acts as a proxy for the raising and lowering of the thermocline in response to locally generated internal tides and other subtidal currents which are not included in our simulations.Ramp et al. (2004) report on the influence of the internal tides on the wave forms in shallow water.
Figure 17 compares waves fields for the three stratifications (ρ b , ρ 1 , and ρ 2 ; see Fig. 3) using bathymetry h 15 , f = 0 and an initial wave APE of 100 MJ m −1 .These plots show the leading waves after they have reached the shelf, along with trailing mode-one and higher mode waves.For stratification ρ b (Fig. 17a) the shoaling solitary waves encounter a critical point and solitary waves are waves of elevation on the shelf.In this case the leading depression is very broad and decays with time.The pedestals will ultimately yield solitary waves of elevation.For stratification ρ 1 (Fig. 17b), the pycnocline is lower in the water column and the critical point is reached at greater depth.In this case breaking has resulted in almost complete destruction of the pedestals.For stratification ρ 2 (Fig. 17b) the pycnocline is above the middepth on the shelf and solitary waves of depression persist.In this case it can be seen that the leading depression has a steep front (at approximately x = 86 km), in contrast to the other cases.Ultimately, a broad square-shaped solitary wave of depression will form.Ramp et al. (2004) observed waves of depression at 120 m depth during times when the internal tide had raised the thermocline above the mid-depth (their Fig. 6, from 18:00 GMT 13 May to 06:00 GMT 14 May), some of which have the appearance of broad waves of depression (e.g. between 23:00 and 02:00).
In Case 3 (ρ b ; Fig. 17a) a packet of short small-amplitude mode-one waves is present between about 60 and 80 km that is absent from the other two cases, while in Case 12 (ρ 2 ; Fig. 17c) there is a short packet of mode-one solitary waves between 60 and 62 km that is not present for the other two stratifications.The mode-two wave field also shows some differences.For Case 9 (ρ 1 , Fig. 17b) a long convex modetwo wave is present between x = 44 and 50 km.For the other stratifications mode-two wave trains have been generated along with concave mode-two waves at x = 53 and 50 for Cases 3 and 12, respectively.
Rotation modifies the wave fields (Fig. 18).The leading mode-one waves on the shelf are similar but smaller in amplitude.The mode-two waves are also similar to those in the corresponding non-rotating cases.The striking differences between the rotating and non-rotating cases are the large  mode-one wave packets in the rotational cases (at x = 67-70, 69-73 and 60-64 km in Cases 3r, 9r, 12r, respectively) that have been generated by the nonlinear evolution of the inertia-gravity waves that form behind the shoaling solitary wave (Grimshaw et al., 2014).The mode-one wave packet was also present in Case 12 but is much larger in the rotating case.

Effects of viscosity and boundary layer separation
In the real world, shoaling waves are subject to a no-slip bottom boundary conditions and the effects of boundary layer instabilities as well as diffusion and dissipation related to other physical processes, in particular tidal currents.To illustrate the potential implications of these processes some simulations were done with vertical eddy viscosity/diffusivity, using a scale height of h s = 10.0 m.The dimensionless function f (z) is shown in the left panel in Fig. 19.The diffusivity/viscosity has a maximum value of K at the bottom and decreases by factors of 10 and 100 at approximately 30 and 53 m above the bottom.This form of the viscosity/diffusivity coupled with a no-slip bottom boundary condition results in flow separation and vortex shedding at the back of the shoaling ISWs (Lamb and Nguyen, 2009;Boegman and Ivey, 2009;Lamb, 2014).The functional form of the eddy viscosity/diffusivity is ad hoc.We take the point of view that physically what is important is a mechanism to create boundary layer separation and vortex shedding off the boundary and we do this in a simple way while confining the viscos- ity/diffusivity to a region close to the bottom boundary.Use of three different values of K illustrate the sensitivity of the results to its value.
Figure 19 compares results from an inviscid simulation (Fig. 19a, b) with results from simulations with K = 10 −5 , 10 −4 and 10 −3 m 2 s −1 .The latter is possibly unrealistically large.At the stage illustrated, the four solutions have just become noticeably different in form (when the lead waves are at x = 59 km all four results are very similar).For K = 10 −5 m 2 s −1 (Fig. 19c, d) the most noticeable difference is the enhancement of the positive near-bottom currents beneath the first part of the pedestal (x = 68.5 km).This is consistent with flow separation and the formation of a vortex behind the leading depression as has been observed in the laboratory and in lab-scale numerical simulations.This positive current behind the leading depression is strengthened as K increases.For K = 10 −4 m 2 s −1 (Fig. 19e, f) the waves are starting to deform.For the largest value of K the leading depression is significantly reduced in length and the trailing waves are smaller and show small-scale features.In shallower water the differences between the four cases grow.
These results and results from other simulations not shown suggest that viscous effects are not significant until the waves reach depths shallower than about 200 m; however, in deeper water the vertical resolution is reduced so these results are merely suggestive.In deep water the near bottom currents are much weaker, which would reduce the effects of the bottom boundary layer beneath the shoaling waves.Observations of sand dunes in the northern South China Sea reported by Reeder et al. (2011) at depths between 160 and 600 m suggest that bottom boundary layer processes could be important at depths greater than 200 m.Our simulations suggest that the presence of a turbulent bottom boundary layer beneath the waves could significantly affect shoaling solitary waves in water depths shallower than 200 m; however, these results are quite sensitive to the eddy viscosity.A more detailed investigation is needed.

Discussion and conclusions
Using a two-dimensional nonhydrostatic primitive equation model we have investigated aspects of the shoaling behaviour of internal solitary waves in the northeastern South China Sea using bathymetry and stratifications based on observations made in the vicinity of the ASIAEX site (Orr and Mignerey, 2003;Duda et al., 2004;Ramp et al., 2004).Sensitivity to the stratification, the bathymetry and the effects of rotation were considered.The majority of the simulations were done without explicit viscosity.A few runs were done to explore the potential implications of a no-slip bottom boundary condition and the presence of a turbulent bottom boundary layer.Consideration of the effects of concurrent barotropic tides on shoaling ISWs are left for future work.
While waves with amplitudes close to 200 m have been observed in the deep part of the South China Sea (Klymak et al., 2006;Lien et al., 2012Lien et al., , 2014)), we considered initial waves with amplitudes between 30 and 115 m in 3000 m deep water, which are more typical (Li and Farmer, 2011).The adiabaticity of the shoaling was considered by tracking the wave amplitude and maximum wave-induced currents as a function of water depth as the waves shoaled and by comparing them with values along the adiabatic curves (ISWs with constant energy).The amplitudes and surface currents during shoaling are significantly less than those for an adiabatically shoaling wave in deep water but can be much larger in very shallow water.Thus, the bathymetry in the vicinity of the ASIAEX site is steep enough to result in significant deviations of the shoaling wave from an ISW solution of the DJL equation with the same initial energy.Some energy is lost by fissioning of the initial wave into multiple waves (e.g.Case 2; Figs. 4, 5).Case 2 also illustrates the common occurrence of wave breaking in waves behind the leading depression (in this case at the back of the third wave of depression).Because of the very steep rear of the shoaling waves, high horizontal resolutions were necessary.For most of our simulations we used a horizontal resolution of 33 m.The trailing waves did depend on the initial water depths.Comparisons of the shoaling behaviour of solitary waves with different deepwater depths were done in which we tuned the amplitude of the wave starting in water of 1000 and 1500 m depths so that in water depths shallower than 1000 m the leading waves had the same amplitude.The wave starting in the deeper water had larger trailing waves, an indication of the non-adiabatic nature of the shoaling even in deep water.Thus, in order to make detailed comparisons with observations it will be important to start at an appropriate deep water depth.
The bumps on the shelf slope topography influenced the shoaling behaviour, particularly the shallower bump on bathymetry h 15 which is about 200 m high in 700 m deep water.As the shoaling wave passes over this bump higher mode waves are formed.As the bumps in the bathymetric transects are slices through three-dimensional features, this implies that three-dimensionality of the bathymetry is likely dergo the same rapid rise and fall of their amplitude as they shoal into shallower water.Ramp et al. (2004) presented time-depth temperature contour plots of internal waves at depths of 350, 200 and 120 m spanning an 11 day period.At 350 m depth they reported wave amplitudes (based on the 24 • isotherm, hence amplitudes potentially underestimated) ranging between 29 and 142 m which matches well with the range of amplitudes in our simulations.Many features of the observed waves are similar to those seen in our simulations although there is naturally a great deal of variability in the observed waves.At 350 m depth both the observed and simulated waves were fairly symmetric.At 200 m depth many of the observed waves and all of the simulated waves were quite asymmetric (see also Fig. 4 in Orr and Mignerey (2003) which shows acoustic images of an asymmetric wave close to this depth).As an example of a wave which is similar to one in our simulations consider the a wave observed at 08:00 GMT on 7 May (see Fig. 4 and Table III in Ramp et al., 2004).At 350 m depth it had an amplitude of 110 m and appears reasonably symmetric (note, however, the difference in stratifications ahead and behind the wave suggesting that the observed wave is propagating at the front of a long depression).It arrived 4 h later at the 200 m depth mooring, at which point it was quite asymmetric with the back of the wave being much steeper than the front of the wave.In our Case 6, which had an initial amplitude of 83.1 m, the leading wave had an amplitude of 117 m at 350 m depth, at which point it was about 1.5 km in width and had a propagation speed of 1.6 m s −1 .Ramp et al. (2004) report a propagation speed of 1.33 m s −1 and a width of 760 m (twice their reported half-width) for their observed wave.At 200 m depth the wave in our simulation had widened to 2 km and the propagation speed had decreased to 1.0 m s −1 .The wave was also asymmetric as is the observed wave.Using a propagation speed of 1 m s −1 the observed wave also has a width of about 2 km.By the time the observed waves reached a depth of 120 m they have already broken up into many smaller waves including waves of elevation superimposed on broad depressions.Our simulations with viscosity suggest that by this depth the waves may have been strongly modified by bottom boundary layer processes.
Many of the observations from the ASIAEX program show waves that are broad, with a gently sloping front and steep back, in water of 120 and 85 m depths (Ramp et al., 2004;Lynch et al., 2004;Duda et al., 2004), often trailed by a small number of narrow waves of elevation and square shaped depressions as in our simulations.Figure 3 in Duda et al. (2004) shows the presence of waves of elevation in which the thermocline is raised above its rest height in a water depth of 85 m, similar to the pedestals generated in our simulations.Ramp et al. (2004) show several examples (e.g.14:00-16:00 GMT on 7 May, 15:00-16:00 GMT on 9 May; their Fig. 6) in which near-bottom water is raised well into the water column in waves of elevation trailing a broad wave of depression in water of 120 m depth.
In our simulations we start with a single solitary wave in deep water.The shoaling wave fissions into several waves, with at least two well-separated solitary waves having been formed by depths of 600-700 m for both bathymetries, with the waves separating at greater depth for bathymetry h 0 which has the more gentle slope.For the steeper bathymetry the waves have reached shallower water before completely separating.In contrast, Lynch et al. (2004) and Ramp et al. (2004) report that large solitary waves at a depth of 350 m have split in two by the time they reach a depth of 200 m.
During the ASIAEX program only one mode-two wave was observed (Lynch et al., 2004).In a pilot program Yang et al. (2004) observed a single mode-two wave in 426 m deep water.Yang et al. (2009Yang et al. ( , 2010) ) analysed mooring data obtained during the VANS/WISE program which was conducted in 2005 and 2006 in the vicinity of the ASIAEX field program.Data from a mooring deployed along the 350 m isobath was analysed and showed the occurrence of many mode-two waves.They appeared in two forms: convex modetwo waves in which there is a bulge in the thermocline with isotherms raised/lowered in the upper/lower part of the water column, and concave mode-two waves with isotherms displaced in the opposite directions.Convex mode-two waves were by far the most common.The mode-two waves trailed mode-one waves suggesting that they were generated via the adjustment of shoaling mode-one waves (Yang et al., 2009).Stratification and seasonal variations played a role in generation of mode-two ISWs with mode-two waves more common in winter (Yang et al., 2009).
Our simulations also showed the formation of both convex and concave mode-two waves which evolved from the shoaling mode-one wave, as observed by Yang et al. (2009Yang et al. ( , 2010)).In the simulations they were present at the depths of the observed waves.Their presence and amplitude were sensitive to the stratification and to rotational effects.The shallow bump on bathymetry h 15 led to the generation of a modetwo wave train and made the mode-two ISWs more prominent (see Fig. 13).The generation mechanisms of mode-two ISWs, their frequency of occurrence, and locations on the continental shelf require further study.
An important process not considered in this study is the influence of barotropic tides and background currents on the shoaling of internal solitary waves.The advection of the shoaling wave by on-shelf and off-shelf tidal currents will decrease and increase the time the waves have to adjust -effectively increasing and decreasing the slope.In addition the on-and off-shelf motion will modify the stratification experienced by the waves and will generate other waves for the shoaling wave to interact with.Our simulations with near-bottom eddy viscosity suggest that the shoaling waves may be significantly modified by viscous effects in water shallower than 200 m, though our simulations are not well resolved in deep water so this remains an open question.Clearly, much remains for us to learn about the evolution of shoaling ISWs.

Figure 1 .
Figure 1.Luzon Strait and South China Sea region.Red square is the location of ASIAEX experimental site.Colours show depth in metres.

Figure 2 .
Figure 2. Sample-measured bathymetries at the ASIAEX location extracted from the digital bathymetry 2 min resolution (DB2) database inside the red square area in Fig. 1.(a) Comparison of 24 different measured bathymetries (grey) with two model bathymetries (solid and dotted black curves).(b) Comparison of first model bathymetry, h 0 , (black) with one of the measured bathymetries (grey).The bathymetries used for simulations starting at depths of 1000 and 1500 m are also shown.(c) Comparison of two other model bathymetries (solid and dashed black curves) with a second measured bathymetry (grey).The solid curve is bathymetry h 15 .The dashed curved is the same bathymetry with the bump at −100 km removed.In (b) and (c) slopes of parts of the bathymetry are indicated.

Figure 4 .
Figure 4. Density fields for Case 2 using bathymetry h 0 .Resolution: x = 33 m, J = 400.Initial wave amplitude 45 m.(a) 0 h (full depth shown).(b) 25 h.(c) 42 h.(d) 50 h.(e) 58 h.(f) 66 h.The initial wave at x = −280 km is barely detectable at this scale.(a) shows the full depth, (b) the upper 800 m and the remainder the upper 400 m of the water column.Density values are σ θ in kilograms per cubic metre.

Figure 5 .
Figure 5. Close up of shoaling waves from Case 2 at t = 66 h shown in Fig. 4f showing onset of breaking behind the leading two waves.Vertical dashed lines are locations of vertical profiles shown in Fig. 6.

Figure 6 .
Figure 6.(a, b) Vertical profiles of horizontal velocity and density for Case 2 at t = 66 h at locations shown in Fig. 5. (a) Horizontal velocity.(b) Density.Curves from second wave pedestal at x = 89 km (solid), the second depression at x = 91 km (dots), the first pedestal at x = 92 km (dashed), and in the leading depression at x = 94 km (dash-dot) (see Fig. 5 for profile locations).The light grey curve in (b) shows the background density profile.(c,d) Profiles in solitary waves of depression (solid) computed on background field given by profiles extracted from simulation in the pedestal at x = 92 km between the two leading waves of depression.Solution is compared with profiles taken from the second depression at x = 91 km (dots).

Figure 5 Figure 7 .
Figure5shows a close up of the leading waves and the wave breaking at t = 66 h.Vertical profiles of the horizontal velocity and density in the wave depressions and wave pedestal at the locations indicated by the vertical dashed lines are shown in Fig.6(the evolution of the shoaling waves suggests the formation of a single pedestal, which in Fig.4eex-

Figure 8 .
Figure 8. Density field for same case and time shown in Fig.7cshowing mode-two waves (in white dashed boxes) and an internal wave beam (above the bump) generated when the solitary wave passes over the bump.

Figure 9 .
Figure 9. Resolution test results.Shown are surface currents at four different depths for Case 3 which uses bathymetry h 15 and an initial wave amplitude of 45 m.Panels compare results from simulations with different horizontal resolutions using 200 grid points in the vertical.Leading wave at 610 m depth -see Fig. 7b.(b) t = 29 h.Leading wave at 370 m depth -see Fig. 7c.(c) t = 35 h.Leading wave at 180 m depth -see Fig. 7d.(d) Results at t = 42 h.First depression in water shallower than 82 m.Second depression at 90 m depth.(e)Zoom in on leading waves at t = 42 h.Vertical dashed lines in (a) indicate the location of the crest of bump and the depth maximum between the bump and the shelf.Depths are indicated by black marks along upper axis.

Figure 10
Figure 10.(a) Total ISW energy (APE plus KE) for ISWs with fixed APE as a function of water depth.Curves are for APE equal to 25 and 50-500 MJ m −1 in increments of 50 MJ m −1 .(b) Ratio of KE to APE for same waves.On right side the ratio increases with APE.

Figure 11 .
Figure 11.Wave properties as a function of water depth.(a) Wave amplitudes (absolute value), (b) maximum surface currents, (c) minimum bottom currents.Grey curves with diamonds are for ISW solutions of the DJL equation for waves with total energy E equal to 50 and 100-1000 MJ m −1 in increments of 100 MJ m −1 .Diamonds indicate computed values.Coloured curves are maximum absolute amplitude and surface current and minimum bottom currents of shoaling waves as function of water depth.Blue -cases using bathymetry h 0 .Redcases using bathymetry h 15 .Black dash-dot -higher-resolution case using bathymetry h 0 .This curve overlies the lower-resolution result (blue curve, initial amplitude 115 m).Both are overlain by the results using bathymetry h 15 for depths greater than about 2600 m.

Figure 12 .
Figure12.Surface current (left) and σ θ = 23.2kg m −3 isopycnal (right) at different times for waves with initial APE of 100 and 60 MJ m −1 started at depths of 3000 (black) and 1000 m (grey), respectively, using bathymetry h 0 .Times are different for the two waves, chosen so that the leading waves are at approximately the same location.Sample local depths are indicated along the upper axis.

Figure 14 .
Figure 14.Density fields for shoaling waves with different initial wave amplitudes using bathymetry h 15 and stratification ρ b .

Figure 15 .
Figure 15.Effects of rotation are illustrated by comparing the σ θ = 23.2kg m −3 isopycnals for two different initial waves.Base stratification with bathymetry h 15 .(a, b, c) Cases 3 and 3r: initial wave amplitude is 45.4 m (APE = 100 MJ m −1 ).(d, e) Cases 5 and 5r: initial wave amplitude is 71.6 m (APE = 200 MJ m −1 ).Solid black curves are results without rotational effects.Solid grey curves are results with f = 5.35 × 10 −5 s −1 .Vertical dashed lines indicate the location of the crest of the bump and the depth maximum between the bump and the shelf.Depths are indicated by black marks along upper axis.

Figure 16 .
Figure 16.Effects of rotation on the amplitude of the leading wave as a function of x.Bathymetry h 15 .Solid red curves: cases without rotation.Dashed blue curves: cases with rotation.Cases 3-7 and 3r-7r.The solid vertical grey lines indicate locations where water depth is 1000, 800, 600 (first occurrence), 350, 200 and 120 m.The vertical dashed lines indicate the locations of the top of the bump at x = −12.4km (480 m depth) and maximum depth (618 m) upshelf of it.