## Abstract

The formation of a recirculating subsurface core in an internal solitary wave (ISW) of depression, shoaling over realistic bathymetry, is explored through fully nonlinear and nonhydrostatic two-dimensional simulations. The computational approach is based on a high-resolution/accuracy deformed spectral multidomain penalty-method flow solver, which employs the recorded bathymetry, background current, and stratification profile in the South China Sea. The flow solver is initialized using a solution of the fully nonlinear Dubreil–Jacotin–Long equation. During shoaling, convective breaking precedes core formation as the rear steepens and the trough decelerates, allowing heavier fluid to plunge forward, forming a trapped core. This core-formation mechanism is attributed to a stretching of a near-surface background vorticity layer. Since the sign of the vorticity is opposite to that generated by the propagating wave, only subsurface recirculating cores can form. The onset of convective breaking is visualized, and the sensitivity of the core properties to changes in the initial wave, near-surface background shear, and bottom slope is quantified. The magnitude of the near-surface vorticity determines the size of the convective-breaking region, and the rapid increase of local bathymetric slope accelerates core formation. If the amplitude of the initial wave is increased, the subsequent convective-breaking region increases in size. The simulations are guided by field data and capture the development of the recirculating subsurface core. The analyzed parameter space constitutes a baseline for future three-dimensional simulations focused on characterizing the turbulent flow engulfed within the convectively unstable ISW.

## 1. Introduction

Internal solitary waves (ISWs) have long been associated with the transport of energy, mass, and momentum in stratified flows. These long nonlinear and nonhydrostatic waves adjust their waveform while propagating over shoaling topography. Enhanced by turbulence inside the wave, shoaling is a major mechanism by which mixing in the interior of the water column is intensified (Shroyer et al. 2011) and particulates are resuspended from the bed (Reeder et al. 2011). Upwelling and turbulent entrainment behind the wave may result in the commingling of planktons, squid, and fish, leading predators to follow the propagating ISWs (Moore and Lien 2007). Shoaling ISWs may profoundly change the properties of the water column, with broader implications for marine habitats and deep-sea exploration.

In our study, an ISW is regarded as a large-amplitude, mode-1 depression of the pycnocline, where nonlinearity is in balance with dispersion. As the ISW shoals, it may lose energy via dissipation or by dispersing into several solitary-like waves. In the absence of significant energy dissipation, the wave may be regarded as shoaling adiabatically—a process that is common over idealized and gently varying bathymetry (i.e., bottom slope *S* of less than 0.03) (Djordjevic and Redekopp 1978; Grimshaw et al. 2004; Vlasenko et al. 2005; Lamb and Warn-Varnas 2015).

When an ISW shoals over steeper slopes (i.e., *S* ≥ 0.03), the wave propagation speed, *c*, decreases below the maximum wave-induced horizontal velocity, *U*_{max}, inducing a steepening of the rear of the wave, overturning, and wave disintegration (Kao et al. 1985; Vlasenko and Hutter 2002; Aghsaee et al. 2010). However, over gentle slopes, the propagation speed may drop below *U*_{max}, followed by rear steepening and heavy fluid plunging forward, yet the wave will not disintegrate and the ISW is said to be convectively, or kinematically, unstable (Hodges 1967; Orlansky and Bryan 1969; Helfrich and Melville 1986). This heavy fluid becomes entrained above the location of the maximum isopycnal displacement, or trough, thereby being locked with the propagating wave. Such a trapped region can be described as a vortex core, or a region with closed streamlines (Derzho and Grimshaw 1997; Aigner et al. 1999).

A closed streamline core contains recirculating fluid. Core structure has been observed in ISWs in the field (Nakamura et al. 2010; Preusse et al. 2012; Lien et al. 2012, 2014; Zhang and Alford 2015; Zhang et al. 2015), experiments (Davis and Acrivos 1967; Grue et al. 2000; Carr et al. 2008; Luzzatto-Fegiz and Helfrich 2014), and simulations (Lamb 2002, 2003; Fructus and Grue 2004; Lamb and Wilkie 2004; Helfrich and White 2010; Soontiens et al. 2010; King et al. 2011; Lamb and Farmer 2011; Carr et al. 2012; Maderich et al. 2015, 2017; He et al. 2019). The core’s convectively unstable nature enhances turbulent mixing and energy dissipation in the water column. Because of its presence, the waves are no longer regarded as shoaling adiabatically. Concurrently, ISWs with a recirculating core may also transport mass across large distances [i.e., *O*(100 km)]. The process by which heavy fluid enters the core may be regarded as “breaking,” but it is not abrupt enough to cause a complete wave disintegration. As such, an ISW with a recirculating core has undergone convective breaking and remains convectively unstable due to overturning induced by the recirculating motion itself as it propagates over the gently varying bathymetry.

Simulations of shoaling ISWs over idealized bathymetry have highlighted the role of the preexisting water column properties, prior to wave passage, in the formation of a convective instability and possible recirculating core. For instance, Lamb (2002) considered the role of the background density, over idealized slope–shelf bathymetry, and argued that a recirculating core can form if there is stratification near the surface. Note that depending on the strength of the stratification, the waves may also be conjugate-flow limited, that is, they become horizontally uniform in their center, in which case a recirculating core may also exist (Lamb and Wilkie 2004). Ensuing work by Stastna and Lamb (2002) and Lamb (2003) examined the role of the baroclinic background current and its significance during the propagation of ISWs.

To date, theory, laboratory experiments, and simulations of recirculating cores in ISWs have focused on a class of such motion regarded as *surface* type, which may be different than that observed in the field. Surface cores reside at the top of the water column, above the trough. According to Lamb and Farmer (2011), a recirculating core forms because the background near-surface vorticity layer in the water column is stretched by the propagating wave, with *U*_{max} increasing past the wave propagation speed. The field observations of Lien et al. (2012) were the first to confirm that ISWs of depression can also support a *subsurface*-type core, located closer to the wave trough. From the solution of fully nonlinear-dispersive theory, based on the Dubreil–Jacotin–Long (DJL) equation (Long 1953; Turkington et al. 1991), He et al. (2019) argued that the primary criterion determining the presence of a subsurface recirculating core is the sign of the near-surface vorticity, associated with the preexisting baroclinic background current, not the density field. Their work further supported the core generation mechanism originally proposed by Choi (2006), using nonlinear asymptotic theory, who found that in a two-layer flow cores could form at the surface if the vorticity in the upper layer is positive (for a rightward-propagating wave of depression) or just above the interface if it is negative. Moreover, He et al. (2019) briefly explored the formation of a subsurface recirculating core in the context of a shoaling ISW, over an idealized bathymetry. Thus, a stratified near-surface layer may result in the formation of a recirculating core and, in the absence of near-surface stratification, a recirculating core may form so long as there is near-surface background shear.

Subsurface recirculating cores have been observed to contain two counterrotating vortices (Lien et al. 2012) which contribute to the mixing of the fluid inside the wave and the dissipation of the turbulent kinetic energy of the wave. Lien et al. (2012) found that the dissipation may be approximately four orders of magnitude higher than that of the open ocean. In addition, the fluid inside the core was observed to be transported with the ISW. According to the aforementioned study, the associated instantaneous mass transport may depend on the size of the core and, in the South China Sea (SCS), it could be as high as 18 Sv (1 Sv ≡ 10^{6} m^{3} s^{−1}). Thus, a subsurface recirculating core is an important mechanism by which mass is transported, fluid is mixed, and energy is dissipated in the water column.

The field records of Lien et al. (2012, 2014) were based on single-point measurements near the Dongsha slope in the SCS, a region known for the active presence of ISWs and the sole location where subsurface recirculating cores have been observed. Questions still persist regarding the process by which shoaling ISWs reach the convective-breaking regime and on how the ensuing subsurface core formation occurs within the wave. Given the computational resources available nowadays, it is possible to complement the field observations on the Dongsha slope with high accuracy/resolution and fully nonlinear/nonhydrostatic simulations. Thus, an objective of our study is to simulate the formation of a convective instability via a high-order spectral multidomain penalty method (SMPM) (Diamessis et al. 2005; Joshi et al. 2016), by capturing the wave as it convectively breaks and to reproduce the formation of the subsurface recirculating core. The work can be guided by observations, including the recorded water column properties and the recorded water depth, thereby bridging the gap between localized observations and the full evolution of the shoaling process.

Given that shoaling ISWs with a subsurface recirculating core have only been observed near the Dongsha slope [albeit these may occur elsewhere (Zhang and Alford 2015; Zhang et al. 2015)], another objective of this study is to highlight the dominant mechanisms that lead to core formation. As such, in regions where large ISWs exist, core presence can be readily established based on the water column properties and, possibly, the local bathymetric profile. To this extent, the two questions guiding this study are 1) what are favorable conditions for the formation of a recirculating subsurface core in a ISW shoaling over gentle slopes and 2) how do variations in the properties of the water column and bathymetry impact subsurface core formation? Numerical simulations in two dimensions address the shoaling problem over a reduced section of the transect spanned by Lien et al. (2014). ISW properties are computed and presented, along with the dimensions of the convective-breaking region. The dissipation of kinetic energy and mass transport are not computed, as these will be the focus of a separate study. Our study aims to establish the foundation for future 3D simulations of a shoaling ISW with a subsurface recirculating core.

This paper is structured as follows: section 2 discusses the method, which includes the background field conditions, the governing equations of an ISW with a recirculating subsurface core, problem geometry, and simulation description; section 3 presents the results, detailing the wave properties for a given initial ISW. The effect of the maximum value of the Dongsha slope is also addressed, as this corresponds to the only region within the SCS where the subsurface cores have been observed. Section 4 explores the variations in the initial conditions where emphasis is placed on the initial ISW amplitude, along with the magnitude of the near-surface shear.

## 2. Method

### a. Field conditions

Figure 1 shows the bathymetry of the region of interest in the SCS, along with the bathymetry originally shown in Lien et al. (2012), in which ISWs were tracked. Within the SCS, Lien et al. (2012) tracked ISWs from 21.07° N, 118.49° E to 21.07° N, 116.50° E. These coordinates describe the track along which ISWs were found to propagate. The observed water depth along the black solid line shown in Fig. 1a is shown in Fig. 1b. The particular choice of bathymetry profile is crucial in dictating the physics of the shoaling problem. Thus, the characteristic bathymetry necessary for subsurface core formation is taken from the actual measured data of Lien et al. (2012), whereas General Bathymetric Chart of the Oceans (GEBCO) data are used to visualize the bathymetry over the greater region in Fig. 1.

The subsurface and surface moorings (Lien et al. 2014) were approximately 6 km apart, located at 21°N, 117.27°E and 21°N, 117.22°E, respectively (see Fig. 1), and cover the region of the steepest slope. Here, the profiles of temperature, density, and velocity were measured prior, during, and after the passage of the ISWs. The location of the moorings is denoted as the black cross marker in Fig. 1a and as the vertical black dashed lines in Figs. 1b and 1c.

The ISWs were tracked over a 2° west-oriented direction, or approximately 200 km. From a modeling perspective, such a long distance can be challenging given the broad range of scales that must be resolved to accurately capture the wave propagating with a subsurface recirculating core. Thus, in this study, a smaller region of interest is extracted, where the shoaling process can be thoroughly analyzed, thereby mitigating the computational overhead. Given that the location where the field observations were made is approximately 1.3°W from the start of the covered track, the shorter region of interest used in the simulations presented here, of approximately 80 km, is focused between 21°N, 117.8°E and 21°N, 117.0°E (see Fig. 1). The bottom slope corresponding to the gently varying bathymetry, for the shorter region of interest, is shown in Fig. 1c as the magenta solid line. The slope is computed along the direction of wave propagation. The maximum slope is 0.028 which occurs very close to the location of the subsurface mooring.

The governing equations discretized to simulate the shoaling process are formulated in a Cartesian coordinate system. To incorporate the observed bathymetry into the model, the SCS coordinates need to be converted from latitude–longitude to universal transverse Mercator. Since there is negligible latitudinal change for a given longitudinal displacement, as one moves along the observational path, the small deviations in the northing direction can be neglected. Thus, in a Cartesian framework, *x* is taken to represent the easting and *z* is the depthwise direction. The beginning of the SCS bathymetric transect can be used as a reference point to create a transect, initiating at zone 50 in the Northern Hemisphere 583 145 m northing, 232 235 5 m easting. The 80-km-long transect has a the water depth varying from 921 m at the deepest location to roughly 360 m at the shallowest location.

### b. Properties of the water column prior to the arrival of the nonlinear internal waves (NLIW)

The subsurface mooring, located at a water depth of approximately 525 m, had 1 upward-looking acoustic current Doppler profiler (ADCP), 10 temperature sensors, and 3 conductivity–temperature–depth (CTD) sensors. The surface mooring was deployed at a water depth of approximately 450 m and contained 2 ADCPs, 14 CTD sensors, and 3 temperature loggers. The spacing between CTDs varied between 10 and 30 m. For a more detailed description of the equipment, the reader is referred to section 2 of Lien et al. (2014). The moorings recorded data from 31 May 2011 to 3 June 2011. In our work, emphasis is placed on 2 June, because on this day the subsurface recirculating core was first observed.

Figure 2 shows the recorded time-averaged background profiles, prior to wave arrival, at the subsurface mooring. These fields are used as the background conditions to simulate wave propagation. Figures 2a–d show the background velocity, shear, density, and Brunt–Väisälä (BV) frequency, respectively. Field observations did not cover the upper 10 m of the water column (see Fig. 1 of Lien et al. 2014), and the data between the free surface and the water depth of 10 m are obtained with linear extrapolation. Our study simulates a wave propagating in the westward direction, toward the coast, and this direction is set to be positive. Note that this sign convention is opposite to that of Lien et al. (2014) as they defined the eastward and westward direction to be positive and negative, respectively. The BV frequency *N* is defined as $N=\u2061(\u2212g\rho o\u22121d\rho /dz)1/2$, where *g* is the gravitational acceleration and *ρ*_{o} is the reference density. The location of the pycnocline, defined as the depth at which the maximum BV frequency occurs, was observed to be at a depth of *z*_{o} = −22 m. At this depth, the density value was 1022.58 kg m^{−3}; the reference density is then set to *ρ*_{o} = 1022.58 kg m^{−3} for the present simulations.

Figures 2a and 2b show the background current profile *U* and the vertical shear profile *U*_{z}. Note that near the surface *U* and *U*_{z} are negative because in this study the eastward direction is taken to be negative. As previously mentioned, the wave propagates in the westward direction such that the ISW-induced vorticity is positive. Thus, subsurface recirculating cores can be expected on this day because the background vorticity is opposite in sign to that associated with the wave (Choi 2006; He et al. 2019). In addition, to avoid hydraulic effects that may be associated with interactions of the background current with the gently changing bathymetry, the original values of *U* (black solid line) below 300 m are smoothed to zero when used in the simulation (blue solid line). Because the water column properties were reported at a single location, any horizontal variations are ignored. This approach assumes that the background current is driven mainly by the internal tides, propagating at similar speed as ISWs, therefore the background field can be treated as steady throughout the simulation.

### c. Governing equations

The governing equations for this modeling study are the two-dimensional incompressible Navier–Stokes equations under the Boussinesq approximation (Kundu et al. 2012). Prior to initializing the flow solver, the velocity field in the horizontal direction is decomposed into a perturbation *u*′(*x*, *z*, *t*) and a steady background field *U*(*z*). The *w*′(*x*, *z*, *t*) is used to describe the full vertical velocity. Per the Boussinesq approximation, the density field is decomposed into a reference value *ρ*_{o}, a background profile $\rho \xaf\u2061(z)$, and a perturbation field *ρ*′(*x*, *z*, *t*).

In vector form, for a fixed reference frame without rotation, the mass conservation and momentum equations are

in the along-wave, *x* direction, with

in the vertical *z* direction, where **u** is the two-dimensional velocity field [i.e., **u** = (*u*′ + *U*, *w*′)], *p*′(*x*, *z*, *t*) is the perturbation pressure with respect to the reference background state, *t* is time, *ν* is the kinematic viscosity, and *g* is the gravitational acceleration, aligned in the depthwise direction; rotation is neglected. During shoaling, the effects of changing water depth may dominate over rotational forces (Lamb and Warn-Varnas 2015). Nevertheless, rotation results in radiation of long inertia–gravity waves, which decreases ISW amplitudes over longer time scales than considered here (Helfrich and Melville 2006; Lamb and Warn-Varnas 2015).

The density equation is

### d. Numerical method

#### 1) Generating the initial conditions from fully NLIW theory

The isopycnal displacement *η*(*x*, *z*), driven by the fully nonlinear ISW and used to initialize the numerical model, is obtained by solving the DJL equation; it is a nonlinear eigenvalue problem derived from the steady incompressible Euler equations under the Boussinesq approximation, in a reference frame moving with the wave in which the flow is steady (Long 1953; Turkington et al. 1991). To solve the DJL equation, the pseudospectral numerical method developed by Dunphy et al. (2011) is employed. Obtaining a solution requires prescribing the background density and current field, along with a target value for the available potential energy (APE). The APE is defined as the energy released in bringing the density field to its reference state (Lamb 2008).

Once the solution of the DJL equation is obtained, the density field is computed via $\rho \u2061(x,z)=\rho \xaf[z\u2212\eta \u2061(x,z)]$. The wave velocity field is computed via spectral differentiation of the isopycnal displacement field. Therefore, the DJL equation provides the density, horizontal, and vertical velocity which are the initial conditions of the unsteady SCS shoaling simulation. More information on the DJL equation, and how to obtain the ISW velocity and density field from its solution, may be found in the appendix.

#### 2) Simulating the shoaling of the ISW

As the foundation of the numerical tool used in this study, the SMPM, originally developed by Diamessis et al. (2005), has been successfully applied to the study of small-scale stratified flow processes (Diamessis et al. 2011; Abdilghanie and Diamessis 2012; Zhou and Diamessis 2015, 2016), with minimal artificial dispersion and diffusion, including the propagation of ISWs in a uniform depth waveguide and two-dimensional studies of their interaction with a model no-slip sea floor (Diamessis and Redekopp 2006). Recently, Joshi et al. (2016) adapted the method to efficiently account for the nonhydrostatic effects with deformed boundaries while preserving high-order accuracy, thereby allowing the incorporation of gentle bathymetry over long domains.

Equations (1)–(4) are solved using the aforementioned two-dimensional deformed subdomain variant of the SMPM. A local Legendre polynomial expansion is used to approximate the solution at each node of a Gauss–Lobatto–Legendre grid in each element (Kopriva 2009). The points are distributed such that finer spacing is achieved near the element interfaces. Time integration is achieved via a stiffly stable third-order scheme (Karniadakis et al. 1991). Nonhydrostatic effects are handled efficiently through a mixed deflation Schur-complement-based pressure solver (Joshi et al. 2016).

The boundary conditions of the SCS shoaling problem are free slip at all four impermeable physical boundaries. In addition, at the left and right boundaries, an artificial Rayleigh-type damper, half ISW-width thick, is applied to eliminate any possible reflection from the incoming ISW (Abdilghanie 2011). For Eq. (4), no-flux boundary conditions are implemented in all four physical boundaries, along with the Rayleigh-type damper at the left and right boundary. Last, an exponential spectral filtering technique is applied to dissipate any numerical instabilities (Blackburn and Schmidt 2003).

### e. Simulation description

#### 1) Obtaining the fully nonlinear ISW field

The background fields discussed in section 2b are used to solve the DJL equation at the initial water depth of *H*_{i} = 921 m because this is the deepest point of the SCS bathymetric transect used in this study (see Fig. 1c). Since the background density was measured down to a water depth of 450 m, any value below this depth is assumed to be constant. Only the near-surface region of the background density is directly linked to the formation of recirculating cores (Lamb 2002). As such, the effects of the near-bottom stratification in the water column are not examined in our study.

The DJL-generated initial condition inserted into the SCS shoaling simulation corresponded to an ISW similar in amplitude to that observed by Lien et al. (2014) near the location of the moorings. Note that this solution may not be representative of the observed wave at a deeper location along the SCS transect as the lack of upstream field measurements of the observed wave impact the accurate representation of the initial wave. Thus, an initial condition resembling the observed wave at a location in shallower waters allows for the exploration parameter space of the shoaling problem which is the focus of this study. The DJL ISW has an initial amplitude of *A*_{i} = 143 m, a width of *L*_{w,i} = 1014 m, and a propagation speed of *c*_{i} = 1.925 m s^{−1}. Note that the subscript *i* is used to denote *initial*.

#### 2) Construction of the computational domain

Figure 3 shows the computational domain used in the shoaling simulation. The computational domain used includes a 20-km-long artificial plateau (i.e., constant water depth), demarcated by the solid red box in Fig. 3, with the SCS transect beginning at a range of 0 km. The plateau is included to allow the ISW to propagate without shoaling for approximately 10*L*_{w,i}. This approach eliminates any nonphysical changes to the waveform that would otherwise occur by placing the initial ISW over actual bathymetry. Aside from the artificial plateau, four distinct locations are also noted in Fig. 3: the initial position of the trough (location I; white dashed line), the surface and subsurface moorings (black dashed lines), and the shallowest portion of the transect (location II; yellow dashed line). These four locations will be used as reference in the subsequent analysis.

To numerically solve Eqs. (1)–(4) with the SMPM, the computational domain is partitioned into *m*_{x} subdomains in the streamwise (*x*) direction and *m*_{z} subdomains in the vertical (*z*) direction, with *n* points per element, or subdomain, in each direction. Together, the total number of degrees of freedom is *n*^{2}*m*_{x}*m*_{z}. The corresponding resolution used in the shoaling simulation analyzed here was determined via a grid-convergence study, where a test ISW was allowed to propagate until the subsurface mooring location and then visually examined for changes in the structure of the solution, as a function of the grid spacing. A test simulation was performed with a fixed *n* value and initial *m*_{x} = 400 subdomains in the horizontal direction and *m*_{z} = 20 subdomains in the vertical direction. Subsequent test simulations were performed up to *m*_{x} = 1600 and *m*_{z} = 30, and no changes were noted in the solution past *m*_{x} = 800 with *m*_{z} = 25. The values of *m*_{x} = 800, *m*_{z} = 25, and *n* = 15 are used in the runs reported in this study. As such, the computational domain has a total of 4.5 × 10^{6} degrees of freedom (DOF).

Table 1 shows the properties of the computational grid employed in this study. Given the total length of the computational domain, the minimum and maximum horizontal grid spacings are Δ*x*_{min} = 2.241 m and Δ*x*_{max} = 13.89 m, respectively, with an effective mean spacing of 9.212 m. In the vertical direction, the spacings varies not just per element but also per horizontal location given the progressive decrease in water depth. At the deepest location, the minimum and maximum vertical grid spacings are Δ*z*_{min,I} = 0.640 m and Δ*z*_{max,I} = 3.967 m, respectively. In contrast, at the shallowest location, the minimum and maximum vertical grid spacings are Δ*z*_{min,II} = 0.254 m and Δ*z*_{max,II} = 1.578 m, respectively. In both the horizontal and vertical direction, the largest grid spacing correspond to the center of the two-dimensional subdomain, or Gauss–Lobatto–Legendre element. No mesh refinement technique is applied throughout the simulations.

The time step size Δ*t* is chosen so as to respect the Courant–Friedrichs–Lewy limit for the initial velocity scale and the grid properties and the limit is set to 0.50 for both the *x* and *z* directions. An adaptive time-stepping method ensures that Δ*t* is adjusted during the shoaling simulation. For reference, Fig. 4 shows the isopycnals at location I, along with the SMPM grid superimposed in gray. There are approximately 60 points horizontally across the ISW and the maximum isopycnal displacement spans 60 points in the vertical. Across the transect, the vertical grid spacing decreases with decreasing water depth as noted in Table 1.

The computational domain is approximately 100*L*_{w,i} long; it is partitioned into overlapping windows which track the ISW as it shoals. Each window is approximately 16*L*_{w,i} long and the overlap region ranges from 6*L*_{w,i} to 7*L*_{w,i}, depending on the waveform since the ISW is adjusting to the gently varying bathymetry. Once the wave reaches the end of a window, a new window is generated that contains part of the original along with the next portion of the domain. The density and velocity fields are then copied inside the overlapping region, from the original to the new window. The total number of windows required for a SCS shoaling simulation is nine. This windowing technique decreases the DOFs to be solved by a factor of 6, thereby accelerating the simulation of the shoaling ISW along the transect.

#### 3) Choice of Reynolds and Schmidt number

The shoaling process encompasses a broad range of scales extending from finescale motion due to convective breaking up to the ISW width and the propagation distance of *O*(100 km). As such, it is computationally prohibitive to simulate the shoaling problem with a field-value Reynolds number, $ReHi=ciHi/\nu $ and Schmidt number, Sc = *ν*/*κ* because of 1) the resolution required to straddle across the gently varying bathymetry over a long propagation distance, 2) the time-averaged profile of the background density and velocity field, 3) the ISW length scales as the wave shoals, and 4) finer scales limited to the formation of the subsurface recirculating core and any finer-scale structure within. The choice of $ReHi$ and Sc must consider these issues so that the two-dimensional parameter space exploration is economical in terms of memory and run-time costs. In this study, these parameters are set to $ReHi=2\xd7106$ and Sc = 1, which are both two orders of magnitude below the corresponding values of the open ocean.

The impact that the chosen $ReHi$ might have on the potential viscously driven deceleration of the ISW over long distances, has been explored by simulating an ISW propagating over a flat domain, for a distance of approximately 15*L*_{w,i}. Subsequently, the wave propagation speed was computed for both inviscid and viscous case, and compared with the theoretical DJL wave propagation speed. The relative difference was found to be *O*(10^{−3}), for both inviscid and finite $ReHi$, along the specified distance.

A three-dimensional study is required for examining the turbulent flow engulfed within the recirculating subsurface core, but not necessarily the propagation of the ISW. Since the field observations of Lien et al. (2014) indicated that, near the Dongsha slope, the waves propagate virtually along the same latitudinal coordinate, a two-dimensional approach to explore the shoaling process and the formation of the subsurface cores is justified since it can also provide insight into possible core dynamics. Although, ISW breaking is clearly an inherently three-dimensional process, the objective here is to explore the conditions that may lead to such breaking. Thus, simulating the shoaling problem with the specified $ReHi$ and a Sc of order unity may be reasonable for exploring the parameter space in this two-dimensional study.

## 3. Results

### a. Wave properties

#### 1) Obtaining the ISW location

The ISW may be tracked by locating the wave trough, which lies between the convergent and divergent zone, where *du*′/*dx* < 0 and *du*′/*dx* > 0, respectively (Chang et al. 2011). Because of the presence of a core, the center of the ISW is defined to be the location where *du*′/*dx* = 0 below the pycnocline. Figure 5a shows the location of the ISW trough as the solid black line with markers and it is captured at every 80 s throughout the simulation. At the initial position, this sampling time corresponds to changes in the trough position of approximately 135 m. As the water depth decreases and the ISW decelerates, the variations in position decrease to approximately 100 m. Error bounds at a given *x* location are obtained by locating *du*′/*dx* = 0 at different water depths. The relative error is found to be less than 1% suggesting that the approach to track the wave is reliable. The error bounds, characterizing the uncertainty in the displacement of the wave, are also included in Fig. 5a as error bars, but given the small difference these are minute and barely noticeable.

Two other regions within the ISW are identified and tracked: the front and lee of the wave. These are defined by first extracting the density profile for a given water depth, in the along-wave direction, then obtaining the streamwise location of the median density value for such profile. The front and lee are shown in Figs. 5a and 5b as the red dotted (lee) and cyan dotted (front) lines, respectively, along with the trough shown as the black solid line. The difference in the location of the lee and front, relative to the position of the ISW trough, is denoted as Δ. Tracking these two points may be a proxy for visualizing changes in the wave symmetry during shoaling. In Fig. 5c, the evolution of Δ along the SCS bathymetric slope is shown. The data indicates that the ISW is sensitive to the varying water depth, particularly at the subsurface mooring where it broadens while propagating over the maximum slope of the transect; this wave broadening is typically observed in the field (Helfrich and Melville 2006). Figure 5 highlights how full nonlinearity and hydrostaticity are effectively enforced in these simulations by relying on DJL theory for the initial wave and by solving the incompressible Navier–Stokes equations over the gently evolving bathymetry, further enhanced by a high-accuracy/resolution method.

#### 2) Determining the ISW propagation speed, amplitude, and width

The position of the ISW trough may be used to determine the propagation speed by performing a least squares fit using a linear model applied to subintervals of the data in Fig. 5a and computing the slope of the linear fit (Moum et al. 2007). For shoaling ISWs, changes in bathymetry may significantly impact the speed calculation because the wave slows down. The linear fit is applied over a subinterval which involves a sufficient range of wave positions, while considering the effects associated with the gentle change in water depth. To capture the propagation speed, the linear least squares fit is set to cover a subinterval of an ISW width, providing an accurate measure of the propagation speed especially near the moorings, where the slope change is more pronounced (see Fig. 1).

In this study, the ISW amplitude, *A*, is taken to be the maximum isopycnal displacement, obtained from *η*(*x*, *z*, *t*). The width, *L*_{w}, is computed by first, integrating *η* in the along-wave direction then dividing by the amplitude (Koop and Butler 1981) as

#### 3) Properties of the simulated ISW

Figure 6a shows the computed propagation speed *c*, along with the maximum ISW-induced horizontal velocity *U*_{max}. The amplitude *A* and width *L*_{w} of the ISW are shown in Figs. 6b and 6c. Figure 6d shows the water column depth to provide perspective of the SCS transect.

As the ISW shoals, the propagation speed and horizontal velocity decrease while the amplitude increases as an increase in amplitude leads to a decrease in width. The maximum amplitude is found to be approximately 153 m, occurring at a range of 55.68 km, at a water depth of approximately 478 m; here, the width is *L*_{w} = 775 m.

The simulated wave propagation speed at the location of the subsurface and surface mooring were 1.66 and 1.50 m s^{−1}, respectively. For comparison, the linear propagation speed for a two-layer water column with an upper layer subjected to constant shear may also be considered. The linear propagation speed may be computed from (Choi 2006)

where *h*_{1} corresponds to the thickness of the top layer, *g* is the gravitational acceleration, *h*_{2} = *H* − *h*_{1}, *ρ*_{1} and *ρ*_{2} are the density in the top and bottom layer, respectively, and *U*_{z,z=0} is the vertical shear at *z* = 0, taken to be constant throughout the top layer.

If one assumes that *h*_{1} = *z*_{0} and taking the average density value at the top and bottom layer of $\rho \xaf\u2061(z)$ as *ρ*_{1} and *ρ*_{2}, with *U*_{z,z=0} = 1.96 × 10^{−2} s^{−1} the linear propagation speed at the subsurface and surface mooring is approximately 1.204 and 1.197 m s^{−1}, respectively. These values are between 20% and 40% below the simulated values, suggesting that linear theory may not be appropriate to simulate the shoaling problem. Nevertheless, from Eq. (7), the effect of the shear in the background current is evident as it increase the propagation speed of the wave.

### b. Examining the presence of a convective instability

When the ISW reaches the location at the two moorings, the propagation speed has already decreased below *U*_{max} and the wave has entered the convective-breaking regime. Figure 7 shows colored isopycnals of the shoaling ISW 1) prior to becoming convectively unstable, 2) at the subsurface mooring location, 3) at the surface mooring location, and 4) at location II where the SCS transect is the shallowest; Fig. 7e shows the transect. In Fig. 7a, the maximum horizontal wave-induced velocity remains below the propagation speed; no convective breaking is noted. In Fig. 7b the isopycnals indicate the presence of convective breaking and subsequent overturning in the water column and that a heavy-over-light fluid configuration has been established. Once the ISW reaches the shallowest part of the transect, the ISW is propagating with an enclosed isopycnal region. Heavy fluid appears to be trapped inside the wave, suggesting the presence of a recirculating core. The snapshots in Fig. 7 demonstrate that overturning is due to the shoaling of the wave.

Figures 6 and 7 indicate that the condition *U*_{max}*c*^{−1} > 1 precedes the formation of the convective instability and, possibly, the recirculating core. That is, once the wave propagation speed decreases below the maximum horizontal wave-induced velocity following a short transitional window, a convective overturn ensues and the formation of a region with enclosed isopycnal subsequently occurs. These findings are consistent with the simulations of Lamb (2002), as well as the field observations of Lien et al. (2012, 2014), where *U*_{max}*c*^{−1} > 1 always preceded the generation of a recirculating core. Note that, visualizing the density field may not be a clear indicator of recirculating fluid. A more robust approach would be to examine the streamline pattern in a reference frame moving with the ISW, which is addressed in section 3c.

#### Comparing simulated ISW properties with observations

The wave properties reported by Lien et al. (2014) include the wave amplitude *A*, the wave width, the propagation speed *c*, and maximum wave-induced velocity *U*_{max}, as found in their Table 1. The observed *U*_{max} was obtained from 1-min averaging of the ADCP data at the subsurface and surface mooring location. The observed propagation speed at the subsurface mooring location was computed by considering the wave arrival time between the mooring and the deployed bottom pressure sensor, separated by approximately 1 km. The propagation speed at the surface mooring was computed as the distance between subsurface and surface moorings divided by the time for the center of ISW to pass the two moorings. The observed amplitude was taken to be the maximum isopycnal displacement, while the observed width of the wave was defined as one-half amplitude following the maximum vertical displacement.

At the subsurface and surface mooring, the observed propagation speed was 2.20 and 1.71 m s^{−1}, respectively, and the corresponding simulated values are 1.66 and 1.50 m s^{−1}. The simulated wave did not decrease in speed as dramatically as the observed wave as it propagated along the location of the moorings (~11% vs ~28%). The simulated value of *U*_{max} in the upper layer is 1.70 m s^{−1}, while the observed value was 2.23 m s^{−1}. The maximum simulated wave amplitude is approximately 153 m and occurred at a depth of 478 m. The observed amplitude was 137 m.

Differences between the observed and simulated wave are expected since no observational input on the upstream conditions is available to select a more representative initial condition. In addition, we have extrapolated observed values of the background current and stratification into the upper 10 m of the water column. The observed wave by Lien et al. (2014) presumably had a different amplitude at the initial water depth of the present study. The range of all possible and stable DJL solutions, with the fields shown in Fig. 2 used as initial conditions, do not yield an ISW with wave properties similar to those observed by Lien et al. (2014) at the mooring sites.

Furthermore, the assumption of a steady horizontally homogeneous background current and stratification, may not be realistic near the Dongsha slope. The ratio *U*_{max}*c*^{−1} could change considerably if a different background current profile is used in deeper waters. Even though the background current is steady, it may have a strong dependence in the normal-to-isobath direction, which can impact the velocity field of the wave as it shoals. Last, the direction of propagation of the observed ISWs near the Dongsha slope may not be straight over the bathymetry shown in Figs. 1 and 3. Wave propagation can vary up to *O*(10°) in the westward direction (see Fig. 4 of Lien et al. 2014), thus impacting the propagation speed and amplitude. Given that, in the present study, any variation in the propagation direction is not captured, the modeled and observed wave properties may differ.

Both the simulated evolution of *U*_{max} and *c* and the properties observed on 3 June shown in Fig. 8a of Lien et al. (2014), exhibit a decreasing trend in value up to the surface mooring. This is the location where the largest difference between observed *U*_{max} and *c* exists and this feature is captured in the simulation. After the surface mooring, the simulated evolution of velocity and propagation speed differs from the field data, and the simulations do not exhibit similar values of *U*_{max} and *c*. Such differences could be attributed to the present study being two dimensional, where there is no physical mechanism by which energy can be dissipated, so that once *U*_{max} increases past *c*, the recirculating core forms and persists as the ISW continues shoaling over SCS bathymetric transect. Alternatively, since the field data are averaged over a specific period, typically *O*(1 min), a close match between simulation and field data may be difficult to achieve. Thus, this simulation captures the essential qualitative aspects of the formation of the convective instability and recirculating subsurface core, but does not quantitatively match the observed wave.

### c. Visualizing the subsurface recirculating core

Figure 8 shows the simulated ISW at the subsurface (Figs. 8a–c) and surface (Figs. 8d–f) moorings, using three different definitions to visualize the fluid inside the wave. In Figs. 8a and 8d, the visualized isopycnal range is saturated to resolve the convectively unstable fluid. Figures 8b and 8e show the contour where *U*_{max}*c*^{−1} = 1 and Figs. 8c and 8f show the streamlines for an observer in a reference frame fixed with the wave. The streamlines are obtained from the streamfunction

where *u*(*x*, *z*) = *u*′(*x*, *z*) + *U*(*z*); arrows are included to denote the flow movement across the ISW. The wave propagates with speed *c*, in the rightward direction, as denoted by the black arrows below the trough.

Figures 8b and 8e may be used to examine the size of the convective-breaking region. At the subsurface mooring, the length and height of the region are found to be *l*_{c} = 180 m and *h*_{c} = 28 m, while, at the surface mooring, the length and height of the core are found to be *l*_{c} = 370 m and *h*_{c} = 45 m, respectively. The streamline characteristics shown in Fig. 8f are similar to those from the observations made by Lien et al. (2012). The simulated ISW has two counterrotating regions which was a distinct feature of the observed ISW with a subsurface recirculating core near the Dongsha slope. Thus, the present simulation captures the formation of a subsurface recirculating core in a shoaling ISW over a gentle bathymetry. Note that from Fig. 8, not all the fluid that is convectively unstable or contained within the region circled by *U*_{max}*c*^{−1} = 1 seems to be effectively recirculating. Recirculating cores are known to leak the trapped fluid into the ambient during ISW propagation. A more robust definition of the core and its boundary can be based on Lagrangian coherent structures (Luzzatto-Fegiz and Helfrich 2014). This method is not currently explored and is left for future simulations of subsurface recirculating cores.

## 4. Discussion

### a. Evaluating the effect of slope near the moorings

The maximum slope along the SCS transect, shown in Fig. 1c, is approximately 0.028 and occurs near the moorings. Realizing that the local slope may play a pivotal role in core formation, a separation simulation is performed with a bathymetric transect where the maximum slope is attenuated. The result is a modified SCS transect with a maximum slope of approximately 0.015, as shown in Fig. 9a, which increases the water depth from 461 to 485 m. The original and modified slope are shown in Fig. 9b.

The maximum horizontal velocity and propagation speed are shown in Fig. 9c. The black and blue line correspond to the original and modified transect, respectively. With the modified transect, the wave still attains the convective-breaking regime since *U*_{max}*c*^{−1} > 1 occurs. However, the location where *U*_{max}*c*^{−1} is the largest changed from approximately 60 km, or *H* = −425 m, for the original slope to 65 km, or *H* = −390 m, for the modified slope. Both simulations suggest that the convective breaking persists with *U*_{max}*c*^{−1} > 1 throughout the rest of the transect. Therefore, a wave propagating over the measured transect is expected to exhibit convective breaking earlier possibly due to the sudden change in water depth, as opposed to a more gradual change in the slope. Figures 9d and 9e show the amplitude and width, respectively, for both original and modified SCS transect. No changes in the ISW length scales are noted.

Figures 10b and 10d show the streamline pattern associated with the ISW, for the modified transect at the location of the two moorings, respectively. Closed streamlines are not present for the modified case, although Fig. 10d suggests that the subsurface recirculating core is in the process of forming; the core eventually forms in shallower water (not shown). In addition, the size of the convective-breaking region may be obtained by considering the contour over which *U*_{max}*c*^{−1} = 1; it is shown as the solid magenta line in Figs. 10b and 10d. At the location of the subsurface mooring, its dimensions are *l*_{c} = 120 m and *h*_{c} = 19 m, while at the surface mooring *l*_{c} = 330 m and *h*_{c} = 40 m. These values are approximately 10%–50% smaller than those from the original transect, suggesting that the amount of heavy fluid plunging forward, into the ISW, is influenced by the presence of the maximum slope and possibly it is unique to the region. Thus, ISW experiencing convective breaking that results in recirculating subsurface cores may be occurring elsewhere, but are not as noticeable as those near the Dongsha slope where the local bathymetric slope enhances core formation.

### b. Sensitivity to initial ISW amplitude

The ISW used to simulate the shoaling problem in section 3, with an initial amplitude of 143 m, may not have corresponded to the observed wave in deeper waters. As such, shoaling simulations are conducted with larger initial waves. Larger initial amplitudes are favored over smaller because Lien et al. (2012) observed a convective-breaking ISW, shoaling over the Dongsha slope, with an amplitude of 180 m at a water depth of 600 m.

One objective of this study is to examine how the properties of the shoaling ISW, and subsurface recirculating core, vary with initial ISW amplitude. Considering the observed water column properties *U*(*z*) and $\rho \xaf\u2061(z)$, along with a water depth of 450 m, from the solution of the DJL equation, Eq. (A2), the minimum amplitude of an ISW with a subsurface recirculating core near the surface mooring is 140 m. On the other hand, the DJL solution corresponding to an ISW with an amplitude of 167 m, at the deepest region in the transect, is the smallest convectively unstable wave and this value is taken as the upper bound of desired initial amplitudes in this study. The DJL solutions encompassing an initial amplitude higher than 143 m and less than 167 m are then selected to simulate the shoaling problem.

Figure 11 shows the ISW properties, as a function of the initial wave amplitude, for the new simulations. The observed ISW properties are also included; the values are obtained from Table 1 of Lien et al. (2014). In addition, given that the field data for 3 June shows a relatively constant *U*_{max} at the location of both moorings (see Fig. 8 of Lien et al. 2014), in the present discussion, the observed *U*_{max} for 2 June it is assumed to be constant between the location of the moorings.

Regardless of initial amplitude, convective breaking occurs for all cases. Furthermore, considering the location where *U*_{max}*c*^{−1} = 1, larger ISWs achieve convective breaking earlier during shoaling. Figure 11 also shows the ISW amplitude and width of the simulations along with the bathymetry.

No simulated internal solitary waves are found to reach the observed values of *U*_{max}*c*^{−1} at the location of the moorings, although all waves attain their maximum value close to the surface mooring. In addition, the simulations also do not reach the observed amplitude. Thus, the initial ISW may not necessarily need to be larger than the originally selected 143-m amplitude wave. It may be possible that the underlying assumptions of this study cannot fully describe the field conditions near the Dongsha slope.

The wave properties at the subsurface and surface mooring locations are shown in Table 2. The data also include the height and length of the convective-breaking region, which are characterized by a height *h*_{c} and a length *l*_{c}. Simulation results indicate that a larger initial ISW, results in a greater isopycnal displacement near the moorings, and a larger breaking region. Waves with an initial amplitude larger than 143 m contain a breaking region larger than that observed. As such, the present two-dimensional simulations, with the recorded background conditions, describe the process by which the ISW experiences convective breaking and a subsequent subsurface recirculating core forms but they do not reproduce the observed wave properties.

### c. Variations in near-surface background shear

The sensitivity in the formation of the subsurface recirculating cores to the background current is explored in this study by modifying the near-surface background current in the upper 20 m. Figure 12a shows the original baroclinic background current profile, *U*(*z*) (solid blue line), along with two new profiles: *U*^{r}(*z*) (dashed blue line) and *U*^{l}(*z*) (dotted blue line). The background shear for all profiles is shown in Fig. 12b. The shear at the surface is negative for all profiles. Note that, subsurface recirculating cores may be expected so long as the background current vorticity is opposite to that of the wave.

The initial waves, with the modified background currents, had the same initial amplitude of 143 m as the base case. Figures 12c–e show the evolution of the wave properties of the ISW, along the shoaling track, with the original and modified background velocity profiles. In Fig. 12c, the profile of the ratio *U*_{max}*c*^{−1} for both modified profiles follows that of the original; all three cases achieve *U*_{max}*c*^{−1} > 1. No significant difference in the amplitude, width, propagation speed, *U*_{max} is noted for the modified profile cases.

Table 3 shows the properties of the convective-breaking region. The length scales slightly vary at the surface mooring location, across cases. For instance, the height, *h*_{c}, is larger for the simulation with *U*^{l}(*z*) than that of *U*(*z*) and *U*^{r}(*z*), with the difference between approximately 6% and 30%, respectively. Thus, the magnitude of the shear at the free surface influences the size of the convective-breaking region and, possibly, the size of the subsurface recirculating core: the larger the magnitude of the background shear, the larger the breaking region.

### d. The potential of a convective instability as a function of the ISW amplitude

According to Helfrich and Melville (1986) and Helfrich (1992), the potential convective breaking of an ISW may be determined by comparing the incident wave amplitude, *A*, with the thickness of the bottom layer of the water column, for a given bed slope. The bottom layer thickness is obtained by subtracting the pycnocline depth *z*_{o} from the total water depth *H*. In this context, three different regimes that describe breaking during shoaling have been proposed: convective breaking if *A*/(*H* − *z*_{o}) > 0.4, shear breaking if 0.3 < *A*/(*H* − *z*_{o}) < 0.4, and no breaking (stable) for *A*/(*H* − *z*_{o}) < 0.3. Vlasenko and Hutter (2002) more recently proposed a convective-breaking criterion, as a function of bottom bed slope, based on shoaling ISW simulations such that the wave amplitude required for overturn may be readily obtained. However, the criterion did not consider gentle slopes, where ISWs have also been observed to experience convective breaking. The aforementioned studies recognized that the ISW-induced horizontal velocity exceeds the wave propagation speed immediately before the wave breaks and that when breaking occurs, waves may transport mass upslope. Examining the wave amplitude for a given bed slope may indicate convective breaking and perhaps hint at the presence of a convectively unstable ISW.

In this study, the above amplitude-based breaking criterion is applied to all simulations, using the background conditions shown in Fig. 2. Figure 13 shows the results for both the original and modified bed slope; the original slope data are shown as the colored cross markers while the modified slope data are shown as the black circle markers. Note that these have a different slope values at the corresponding mooring locations. The threshold values proposed by Helfrich (1992) and Vlasenko and Hutter (2002) are specified as the red and black solid lines, respectively. All simulated ISWs with the original slope reach the approximate convective-breaking value of 0.4 at the surface mooring location, though not at the subsurface mooring location. The field observations of Lien et al. (2014) reported an approximate value of 0.4 near the moorings for their 2 June wave.

However, the modified slope results do not reach the proposed convective instability limit. At the subsurface mooring location the wave is considered stable while at the surface mooring location it is within the shear instability region. Nevertheless, as shown in Figs. 9 and 10, the wave does experience convective breaking since *U*_{max}*c*^{−1} > 1. Hence, convective breaking depends on the preexisting background current and density field (Lamb 2002; Stastna and Lamb 2002; Lamb 2003; Soontiens et al. 2010), and not only on the wave amplitude.

## 5. Conclusions

The shoaling of an internal solitary wave of depression, over realistic bathymetry and with actual field background conditions, has been simulated by solving the incompressible Navier–Stokes equations under the Boussinesq approximation in two dimensions, with a high-order spectral multidomain penalty method. The bathymetry, density, and background current fields in the water column were measured by Lien et al. (2014), near the Dongsha slope in the South China Sea. In this study, particular emphasis has been placed on the formation and evolution of recirculating subsurface cores during the shoaling process, within the constraints of two-dimensional dynamics where there is no turbulence-driven viscous dissipation and fluid mixing; this being inherently three dimensional. The subsurface cores form because of the sign of the relative vorticity associated with the background current, which is opposite to that of the propagating ISW. This study examined variations in the wave properties that support the existence of these cores, as observed in the SCS.

The initial conditions representing the ISW were obtained from the solution of the DJL equation using observed density and background velocity fields. The shoaling simulation indicated the presence of subsurface recirculating cores, as observed in the field. The background fields are steady in time and only vary with depth.

The location of the subsurface and surface mooring, deployed by Lien et al. (2014) along the SCS transect, was used to guide the subsequent analysis. The simulations explored the sensitivity of core formation to the initial wave amplitude, near-surface background shear, and maximum slope within the SCS bathymetric transect. Rotational effects were not included, as it has been shown by Lamb and Warn-Varnas (2015) that the effects of the changing water depth is the dominant factor during the convective-breaking process.

For an initial ISW with an amplitude of 143 m, results indicated the presence of a heavy-over-light fluid configuration and a convective-breaking region with a height slightly below the observed value. Streamline visualization showed the presence of two regions of closed streamlines, characteristic of recirculating fluid driven by two counterrotating vortices. This streamline configuration was observed in the field by Lien et al. (2012), for an ISW with a subsurface recirculating core.

In addition, the impact of the maximum slope at the SCS bathymetric transect was explored by performing a simulation with an attenuated slope value, near the location of the subsurface and surface moorings. Results showed that an earlier arrival in shallower water expedited the formation of the subsurface recirculating core. The core’s unique streamline pattern was not obtained for the simulation with the modified slope at the locations of interest, although it was in the process of forming; no changes in the length scale of the wave were noted. Thus, the sign of the near-surface vorticity, in combination with the local bathymetric profile, results in the formation of subsurface recirculating cores that are easily noticeable near the Dongsha slope.

Variations in the initial conditions were also studied by changing the initial amplitude and inserting larger waves. Larger waves resulted in larger convective-breaking regions, although no field observation exists that could be used to corroborate these results. Furthermore, none of the simulated ISWs matched the amplitude, maximum horizontal velocity, and wave propagation speed observed by Lien et al. (2014) at the surface and subsurface mooring or from shipboard measurements; the amplitude was overpredicted while the velocity scales where underpredicted. This inconsistency may be attributed to the choice of the initial wave for the present study, which may not correspond to the upstream conditions that resulted in the observed wave by Lien et al. (2014).

The effect of the shear, associated with the profile of background current, was studied by modifying its near-surface magnitude while preserving its sign. These simulations indicated that while the ISW’s amplitude and velocity scales did not vary from the simulation with the original background current profile, the size of the convective-breaking region near the surface mooring did. This resulted in a larger magnitude of the near-surface background shear and a larger convective-breaking region and potentially recirculating subsurface core.

This study has concluded that subsurface recirculating cores may exist wherever large-amplitude ISWs propagate in the presence of a background current that has near-surface shear opposite in sign to the ISW-induced vortical motion in the water column. Given that the simulations suggested a persistent convective-breaking region, future studies will incorporate a spanwise direction, thereby allowing for three-dimensional turbulent dynamics, to examine the dissipation of kinetic energy and mixing inside the recirculating subsurface core. Field observations indicate that the wave’s interior eventually becomes stable, since *U*_{max} < *c* is recovered along the SCS transect. In addition, the ISWs have been observed to decrease in amplitude as they propagate over the Dongsha slope, with the subsurface core, potentially due to energy dissipation. Last, the definition of the recirculating subsurface core boundary will be addressed via incorporation of a particle tracking technique. While visualizing the density field or the streamfunction indicates the presence of trapped fluid, it does not provide an accurate method of determining core boundary. As such, the mass transport capacity of the subsurface recirculating cores in shoaling ISWs over gentle slopes can be more systematically quantified.

## Acknowledgments

The authors thank Marek Stastna (University of Waterloo), Yangxin He (University of Waterloo), Kraig Winters (Scripps Instititution of Oceanography, UCSD), Frank Henyey (University of Washington Applied Physics Laboratory), and Kristopher Rowe (Cornell University) for their insightful comments and suggestions regarding recirculating cores and this work. The authors also acknowledge the thorough and candid feedback of the reviewers. Financial support is gratefully acknowledged from the Cornell University Graduate School through the Provost Diversity Fellowship and National Science Foundation Division of Ocean Sciences (OCE) Grant 1634257. This work is dedicated to the memory of Sumedh M. Joshi.

### APPENDIX

#### Dubreil–Jacotin–Long Equation

For a reference frame defined in *ξ*–*ζ* coordinates, moving with the ISW under consideration, steady, incompressible, and inviscid flow, the DJL equation (Long 1953) can be derived from the incompressible Euler equations under the Boussinesq approximation (Turkington et al. 1991). If the effect of a background velocity profile *U*(*ζ*) is included, the DJL equation can be expressed as

where *η*(*ξ*, *ζ*) is the isopycnal displacement, *c* is the wave propagation speed, and *N*^{2}(*ζ*) is the squared BV frequency. With the isopycnal displacement, the density *ρ*(*ξ*, *η*), horizontal velocity *u*(*ξ*, *ζ*), and vertical velocity *w*(*ξ*, *ζ*) fields of the ISW can be obtained from

respectively. The equation can be solved iteratively by prescribing a steady background current and background density profile, along with a range of APE values. For every APE, there is a corresponding isopycnal displacement field, from which Eqs. (A2)–(A4) can be applied to obtain the ISW-induced fields.

## REFERENCES

*Theor. Comput. Fluid Dyn.*

*J. Fluid Mech.*

*Fluid Dyn. Res.*

*J. Comput. Phys.*

*Phys. Fluids*

*Phys. Fluids*

*J. Atmos. Oceanic Technol.*

*Phys. Fluids*

*J. Fluid Mech.*

*Phys. Fluids*

*J. Phys. Oceanogr.*

*J. Comput. Phys.*

*J. Fluid Mech.*

*J. Phys. Oceanogr.*

*Nonlinear Processes Geophys.*

*J. Fluid Mech.*

*J. Phys. Oceanogr.*

*J. Fluid Mech.*

*J. Fluid Mech.*

*J. Fluid Mech.*

*J. Fluid Mech.*

*Annu. Rev. Fluid. Mech.*

*Nonlinear Processes Geophys.*

*J. Geophys. Res.*

*J. Comput. Phys.*

*J. Fluid Mech.*

*J. Comput. Phys.*

*J. Fluid Mech.*

*J. Fluid Mech.*

*Implementing Spectral Methods for Partial Differential Equations*

*Fluid Mechanics*. 5th ed. Academic Press, 920 pp

*J. Fluid Mech.*

*J. Fluid Mech.*

*J. Fluid Mech.*

*Phys. Fluids*

*J. Phys. Oceanogr.*

*Nonlinear Processes Geophys.*

*J. Phys. Oceanogr.*

*J. Phys. Oceanogr.*

*Tellus*

*J. Fluid Mech.*

*Fluid Dyn. Res.*

*Nonlinear Processes Geophys.*

*Mar. Mamm. Sci.*

*J. Phys. Oceanogr.*

*Cont. Shelf Res.*

*J. Geophys. Res.*

*PLOS ONE*

*Mar. Geol.*

*J. Geophys. Res.*

*Phys. Fluids*

*Phys. Fluids*

*Stud. Appl. Math.*

*J. Phys. Oceanogr.*

*J. Geophys. Res.*

*J. Geophys. Res.*

*J. Geophys. Res.*

*Phys. Fluids*

*J. Fluid Mech.*