A Model Sensitivity Study of Ocean Surface Wave Modulation Induced by Internal Waves

Nonlinear internal waves are an upper ocean phenomenon that drives horizontal surface current gradients, which in turn modulate ocean surface waves. Under certain conditions, this wave‐current interaction creates ocean surface roughness heterogeneity, in the form of alternating rough/smooth bands. In this study, we investigate the sensitivity of the modulation effect to internal wave properties and develop sea states using simulations of individual internal wave solitons. We utilize a phased‐resolved two‐layer fluid model to capture the evolution of surface waves deterministically. We conduct extensive simulations with a wide range of parameters, including fluid layer density ratio, internal wave amplitude, and parametric wind speed. We use the ratio of the mean surface slope between the rough and smooth bands, which are identified in the simulated surface wave field, to systematically investigate their response to the internal wave forcing across all our simulation cases. Our results show that, among the internal wave parameters, the upper‐lower layer density ratio causes the strongest surface heterogeneity. Spectral analysis of the surface wave elevation and slope variance reveals that the wavenumbers above the peak are most impacted. We demonstrate that accounting for the internal wave‐induced modulation requires including wave steepness statistics, for example, when modeling air‐sea exchange using a surface roughness, z0, parameter. Currently, these statistics are not included in typically coupled modeling schemes, and these systems cannot account for the impact of internal waves, even if the solitary wave phenomena are resolved.

Nonlinear internal waves (NIWs) (λ > 100 m) are an upper ocean phenomenon that drives surface current gradients, which directly modulates surface waves contributing to z 0 . The primary mechanism for this internal-surface wave interaction has long been hypothesized (Dietz & Lafond, 1950;Ewing, 1950) and, subsequently, theoretically described (Phillips, 1980). Active remote sensing systems have greatly expanded our appreciation of NIWs prevalence and characteristics across the global oceans, both from satellite (Alpers et al., 1981;Apel & Gonzalez, 1983;Apel, 2002;Alpers, 1985;Brandt et al., 1999;Hsu & Liu, 2000) and ship-borne systems (Nyman et al., 2019). Recently, it has been confirmed that the roughness heterogeneity driven by NIW surface currents enhanced local wind shear and stress within the marine atmospheric surface layer (MASL) . Through a case study of NIWs and MASL variables observed from the R/P FLIP, the authors showed that wind shear within the lowest 10-m of the atmosphere adjusted (mechanically) to the NIW-induced surface roughness heterogeneity. The interaction stimulated air-sea momentum flux over the entire NIW packet (approximately seven individual solitons with surface signatures) comparable (in magnitude) to the momentum flux caused by the onset of diurnal sea-breeze. This result suggests that NIWs are an important air-sea exchange mechanism, akin to the recent findings of enhanced air-sea interaction near oceanic submesoscale fronts (Shao et al., 2019). With NIWs, there is the added component of transience, which may distribute the impacts across the lifetime of NIW events.
Hao and Shen (2020) (hereinafter denoted HS20) showed that the NIW-induced surface slope heterogeneity could be reproduced using a wave-phase-resolving two-layer hydrodynamic model. In the study, the free surface was initialized with a known wave spectrum (e.g., JONSWAP) and this broadband surface spectrum was allowed to interact downstream with the internal soliton wave form (Helfrich & Melville, 2006). The results of high-resolution simulation reproduced the surface wave heterogeneity (or "zebra stripes") associated with NIWs in the oceans. The direct simulation, based on the two-layer model, has thus been proven an effective tool for quantifying the surface slope variations at a reasonable computational cost.
The recent observational and numerical studies described above demonstrate that (a) the presence of NIWs induces surface roughness variability that directly influences the air-sea exchange process and (b) the NIW-surface wave interaction can be robustly simulated in a computational domain. The study in the present paper uses the model framework developed by HS20 to conduct a series of controlled experiments exploring the interaction between an NIW and surface wave spectrum under different regimes. The goal of the analysis is to understand under what conditions the NIW-induced surface roughness heterogeneity is produced and quantify its relative strength using the response of the surface wave field. Understanding the sensitivity of the ocean surface to NIW forcing will provide insights into the prevalence of the NIW-atmosphere interaction process first described by Ortiz-Suslow et al. (2019) and the role this mechanism plays in the overall atmosphere-ocean coupled system.

Methods
A major challenge to previous numerical studies of the surface response to internal waves was simulating the free surface gravity wave field. Generally, there are two different approaches to simulation-based studies of the short surface waves: (a) phase-averaged models, where the spatial propagation of the surface waves is tracked similarly to geometric optics and (b) phase-resolved models where the free surface undulations are directly captured. In the phase-averaged models, the surface motions induced by the NIWs are simplified by using a spatially-varying surface current (Alpers, 1985;Donato et al., 1999;Peregrine, 1976 is described by an advection equation, also known as the radiative transfer equation (Young, 1999). This framework is limited by the key assumption that surface waves change slowly. While this assumption holds for waves longer than ∼10 m, it fails for the short waves, which make up the surface roughness and are critical to the interaction with the MASL. These short waves can be blocked and reflected by the internal wave induced surface current field. Consequently, the phase-resolved perspective has seen more applications in recent studies (e.g., Jiang et al., 2019;Taklo & Choi, 2020).
Practically, integrating the Navier-Stokes equations directly is too computationally demanding for this particular problem because of the vast divergence of scales between the NIWs and the free surface waves, if trying to resolve the short scale wave response. Alternatively, the two-layer ocean model is based on the potential flow assumption and is thus less computationally expensive. The nonlinear interaction between surface waves and internal waves has been studied by several simplified models derived from the two-layer model, including the nonlinear Schrödinger equation coupled to the Korteweg-de Vries (KdV) equation (Craig et al., 2012) and the second-order approximation to the original governing equations (Jiang et al., 2019;Taklo & Choi, 2020). These simplifications mean that the short wave dynamics are not fully resolved. Recently, the numerical evidence of the surface slope variations induced by NIWs was presented in HS20, where they directly solved the governing equations of the two-layer ocean model using a high-order spectral method (Alam et al., 2009). This development enabled the two-layer model to be a robust and effective tool for studying this process, without incurring an unreasonable computational cost. Here, we describe the two-layer model and internal wave form used for the sensitivity analysis in the present study. The methodology used for this study builds on recent work and further details regarding the model framework can be found in Alam et al. (2009) and HS20.

Two-Layer Model
In the two-layer model, the fluid is assumed to be inviscid, irrotational, and incompressible. Thus, the fluid motions in each layer are governed by the Laplace equation of the velocity potential. Also, it is assumed that the fluids are immiscible at the interface. The governing equations of the system are: where ϕ u and ϕ l denote the velocity potential in the upper and lower layers, respectively, x and y are the horizontal axis, and z is the vertical axis. We define h u and h l as the upper and lower layer depths, respectively, and η u and η l as the displacements of the surface wave and the internal wave, respectively. The kinematics and dynamics boundary conditions at the surface (z = η u ), the interface (z = −h u + η l ), and the bottom (z = −h u − h l ) are.
[ where ρ u and ρ l are the fluid density in the upper and lower layers, respectively, the horizontal gradient operator is ∇ = (∂/∂x, ∂/∂y), the subscript z denotes the spatial derivative in z − axis, and R = ρ u /ρ l is the density ratio between the layers.
Following Alam et al. (2009), we define two interim variables at the surface and the interface: ) . Then the boundary conditions shown above are transformed to four evolution equations.
The eigenmodes of the two-layer fluid system are used for generating the initial condition in our simulations. By neglecting the nonlinear terms in the evolution equations, we can obtain the eigenmodes of the system.
where a = a(ω, k x , k y ) is the surface wave amplitude, ω is the wave angular frequency, k x and k y are the wavenumbers in x-and y-axis, respectively, φ = k x x + k y y − ωt is the phase function, and = √ 2 + 2 is the wavenumber magnitude. The dispersion relation of the two-layer system is ( + coth coth ) 4 − (coth + coth ) 2 + (1 − ) 2 2 = 0.
For each k, the dispersion relation yields two solutions for the frequency, ω. For the higher frequency solution, the surface wave and the internal wave components are in-phase. For the lower frequency solution, the surface and internal wave components are out of phase. These two different types of eigenmodes are known as the (barotropic) surface wave mode and the (baroclinic) internal wave mode, respectively (Apel, 1988).
The nonlinear wave dynamics in the two-layer system are uniquely determined by four surface/interface quantities: η u , η l , , and ψ i . We construct the initial condition by the superposition of the eigenmodes shown above. At each step, the vertical derivatives, ϕ u,z and ϕ l,z , are calculated by iteratively solving the boundary-value problems at each perturbation order from the surface and interface quantities. Then the linear terms are calculated in Fourier space and the nonlinear terms are calculated in the physical space. The evolution equations are advanced in time using a fourth order Runge-Kutta scheme. More details of the numerical methods are given in HS20.

Surface Wave Field
The initial surface wave field is constructed from the Joint North Sea Wave Project (JONSWAP) spectrum (Hasselmann et al., 1973) where α is the Phillips parameter, ω is the wave frequency, ω p is the peak wave frequency, λ p is the peak wavelength, T p is the peak wave period, γ = 3.3 is the peak-enhancement parameter, σ = 0.07 for ω ≤ ω p , σ = 0.09 for ω > ω p , and D(θ) = 2 cos 2 (θ)/π with ∈ [ − ∕2, ∕2 ] being the angular spreading function. The parameter α and wave frequency ω p are calculated as = 0.076 where g is the gravity acceleration as 9.8 m s −2 , U 10 is the wind speed at 10 m above the mean sea surface, and F is the fetch.
The corresponding directional wavenumber spectrum is where k = (k x , k y ) is the wavenumber vector. The initial surface wave field is constructed as the superposition of wave components with amplitude where Δk x and Δk y are the discretizations of wavenumber k x and k y and Ψ is a random variable uniformly distributed over [−π, π].

Internal Solitary Wave
In the two-layer ocean model, the dynamics of the internal solitary wave are determined by several key parameters, including the amplitude, the density ratio, and the layer depths. In our set-up, these parameters are used to construct the initial soliton, based on the KdV theory, which assumes that the amplitude is smaller than a vertical length scale, for example, h u ( Figure 1). When the amplitude is too large, the perturbation expansion used in our numerical scheme may fail to converge. To avoid numerical instability, we restrict our discussions to the weakly/ moderately nonlinear internal waves comparable to those in HS20, such that both the surface and interface quantities can be expanded in terms of the typical wave slopes. Note that this is a necessary simplification in the phase-resolved models to ensure that the surface wave dynamics can be captured together with the internal wave (Craig et al., 2012;Hao & Shen, 2020;Jiang et al., 2019;Taklo & Choi, 2020).
The initial internal wave profile is described according to the solitary wave solution of the KdV equation. The elevation of the interface is defined by the amplitude a iw , the internal wave width λ iw , and the upper and lower layer depths h u and h l (Kodaira et al., 2016): 10.1029/2022EA002394 where c iw is the velocity of the internal wave, and c 0 is the linear long wave celerity.

Quantifying the Response at the Ocean Surface
HS20 demonstrated that a two-layer model system can be used to capture the surface wave response to internal wave forcing: a pair of bands of smooth/rough surface propagating together with the internal soliton. In field observations, these surface features have been well-documented and directly linked with NIWs (Dietz & Lafond, 1950;Ewing, 1950;Ortiz-Suslow et al., 2019). A similar phenomenon was also recreated in laboratory (Kodaira et al., 2016). Theoretically, it is recognized that these bands are caused by the internal wave-induced surface currents featured by the convergence and divergence zones in phase with the internal wave form (Phillips, 1980).
To quantify the surface response, surface slope, S(x, y), is defined as: We use the spanwise-averaged local slope to identify the smooth and rough regions as where L y is the domain size in the spanwise direction. Using a low pass filter of ( ) (with cut-off wavenumber k = 0.14 rad/m to filter out the high-wavenumber oscillation associated with surface waves of k p = 0.63 rad/m), we can better distinguish the smooth and rough regions associated with the internal wave ( Figure 2). The low-pass filtered spanwise averaged slope is denoted as ̃ . HS20 used visual inspection to identify the smooth and rough bands, but for this study, we develop a systematic approach. First, we normalize the filtered, span-wise averaged slope ̃ by the initial value S 0 , that is, ̃∕ 0 . Then we search over x for the local minimal and maximal of ̃∕ 0 , of which the locations are denoted by x min and x max , respectively ( Figure 2c). We assign these locations as the centers of the smooth and rough bands, respectively. As shown in Figure 2, the width of the smooth and rough bands denoted as 2ϵ are approximately 5λ p , where λ p is the peak surface wavelength. From this, we can compare the mean surface slope within the rough and smooth bands for a given simulation time step by the ratio: This quantity provides a simple proxy for the surface slope heterogeneity induced by the internal wave-posited by Ortiz-Suslow et al. (2019) as the main driver of NIW interaction with the air-estimated directly from the simulated free surface without having to employ a particular model of z 0 .

Results
In this section, we first illustrate the qualitative features of the surface response to the internal wave and then present the surface sensitivity from two aspects: the free surface responses to the variations in the internal wave properties and the surface wave properties. For the former, we use the same initial surface waves as HS20 and systematically vary the characteristics of the underlying soliton. For the latter, we hold the internal wave properties the same, and vary the parameters of the JONSWAP spectrum used for constructing the initial surface wave field to simulate the case of a ramping up wind speed from 0.5 to 5 m s −1 .

An Example of the Free Surface Response to an Internal Wave
In this study, we analyzed several simulations with varying internal and/or surface wave characteristics. To highlight the key surface response to internal waves, we will first present an example of the time evolution of the surface elevations and slope for the case: ρ u /ρ l = 0.997, a iw /h u = −0.80, and initial surface waves with parameters U 10 = 2.82 m s −1 and F = 23.7 km (Figure 3). The former two parameters are the density ratio of the upper and lower layers in the two-layer model, and the ratio of the internal wave amplitude to upper layer depth, respectively. The details of the other simulations will be discussed in the next section, but the overall time evolution is qualitatively similar across simulations.
At the initial time instant t/T p = 0 (T p is the peak surface wave period), the surface wave field is spatially homogeneous (Figures 3a and 3b). As time increases to t/T p = 10, a distinct heterogeneity in S(x, y) emerges-this is the development of the smooth/rough band pair over the soliton trough at x = 0 in the figure. Visually, this is more difficult to discern in the elevation field. By t/T p = 40, the smooth and rough bands have formed and are evident in both the elevation and slope fields. In Figures 3e and 3f, the smooth (rough) band appears constrained to the region −5 < x/λ p < 0 (0 < x/λ p < 5).
The ratio given in Equation 29, R s , provides a summary of the temporal evolution of the surface slope ( Figure 4). The evolution of R s has two distinct phases: (a) a fast ramp-up from the initial value R s = 1 (i.e., S is homogeneous) and (b) an equilibrium state where R s is quasi-stationary in time. In this case, we found that phase one occurs within ∼20T p . This evolution is consistent with the simulation analyzed in HS20 using the two-layer model and was typical of the simulations analyzed in the present study. The equilibrium value reached in phase two appeared to correlate with the initial internal wave and surface wave characteristics. Thus, to provide a single metric for analyzing the various simulations, we defined 〈R s 〉 as the time average of R s during phase two. The quasi-stationary subrange of time duration used to calculate 〈R s 〉 is determined in three steps: (a) use a seventh order polynomial filter to smooth the noisy R s , (b) find the starting time of quasi-stationary subrange as the first time instance where the filtered ∂R s /∂t < 0, and (c) find the ending time of quasi-stationary subrange as the next time instance where R s exceeds ±5% of the R s value associated with step (2). The final 〈R s 〉 was calculated as the time-averaged R s during the identified quasi-stationary time duration, as shown in the orange line in Figure 4. The phase after the quasi-stationary subrange in Figure 4 is caused by normal numerical fluctuations.
To illustrate this, in Figure 4b, we have shown the evolution profile of another case of similar parameters with ρ u /ρ l = 0.997, a iw /h u = −0.71 which shows an opposite falling trend in the third phase compared to Figure 4a. We anticipate that the fluctuations can be mitigated by an ensemble run of random surface wave fields in future study.

Effect of Varying Internal Wave Profile
In this section, we analyze the sub-set of cases where the internal wave form is varied beneath the same initial surface wave field. In this sub-set, the initial JONSWAP spectrum is set up with U 10 = 2.82 m s −1 and F = 23.7 km, resulting in a surface wave field with peak wave properties ω p = 2.5 rad/s, λ p = 10 m, and T p = 2.5 s. The computational domain size is set to L x × L y = 500 × 125 m, with 6144 and 1536 grid points in the x and y directions, respectively, and the computational time step is 0.07 s. The shortest resolved surface wave has a wave length of 0.25 m and a wave period of 0.4 s. Therefore, the computational time step 0.07 s is sufficiently small to simulate the surface wave evolution. These initialization parameters are consistent with HS20. For all simulations, the upper layer depth (h u ) is held constant and for this analysis, we have limited our scope to focus on the simulation time when R s is quasi-stationary (e.g., the orange line in Figure 4).

Case Setup
We performed 14 simulations with different density ratio ρ u /ρ l and normalized internal wave amplitude a iw / h u , as shown in Figure 5. The case from HS20 (ρ u /ρ l = 0.997, a iw /h u = −0.71) is marked with a red cross. The density ratio ρ u /ρ l varies from 0.975 to 0.997, which corresponds to an upper layer density varying from the density of pure water to seawater (i.e., 1,003.0-1,025.6 kg/m 3 ). The normalized internal wave amplitude a iw / h u spans −0.8 to −0.2.
The internal wave driven surface current is approximated to the leading order using mass conservation (Kodaira et al., 2016): Equation 30 is the Eulerian current obtained by considering the mass conservation of the flow at infinite and the flow at the cross-section of the internal wave trough with the free surface displacement η u neglected, and therefore it does not include the residual current such as the Stokes drift. Based on Equation 30, Figure 5 shows the different simulations over the internal parameter space and contoured by the maximum stream-wise gradient, that is, |dU s / dx| max . It is evident that the maximum gradient occurs for large internal wave amplitudes at the interface of fresh and seawater (i.e., a relatively large density difference).
In Figure 6, we plot the internal wave profiles η l , the surface current U s , and the current gradient dU s /dx. For clarity, we only plot the nine cases that share  . Density difference has little impact on the internal wave form itself (Figure 6a) but does strongly impact the surface currents and their maximum gradients during the subrange of time when R s reaches equilibrium (Figures 6b and 6c), as compared to internal wave amplitude. As a result, we expect the surface elevation and slope fields to respond more strongly to varying ρ u /ρ l . For example, for a constant a iw /h u = −0.71, a change in ρ u /ρ l from 0.997 to 0.99 (a <1% change) results in a doubling of the U s over the internal wave trough calculated from Equation 30.

Surface Wave Modulation
To qualitatively view the effect of internal wave properties on the surface response, we plot surface wave elevation η u and surface wave slope S during the quasi-steady time range for three representative cases in Table 1, as Note. The three cases, that is, IW-A, IW-B, and IW-C, are the smallest, median, and largest 〈R s 〉 among all the cases shown in Figure 5. shown in Figure 7. The three cases, namely, IW-A, IW-B, and IW-C, represent the smallest, median, and largest 〈R s 〉 among all the cases shown in Figure 5. In the case IW-A (a iw /h u = −0.2), the surface field is insensitive to the presence of the internal wave, with homogeneous elevation and slope fields (Figures 7a and 7b). For increased wave amplitude in the case IW-B (a iw /h u = −0.71), the internal wave is evident in S but remains difficult to visualize simply in the elevations. Upon close inspection, Figure 7c does exhibit some stretching and bunching of the surface waves near x/λ p = 0, but it is not immediately obvious. By changing ρ u /ρ l to 0.975 (a net 2.2% decrease) in the case of IW-C the typical smooth/rough banding above the internal wave trough is visually clear in both the elevation and slope fields. In particular, the high wavenumber contributions are almost entirely wiped out in S(x, y) within ∼5 surface peak wavelengths from the internal wave trough.
To investigate the response of the surface waves at different length scales, we analyze the surface elevation variance density and slope spectra within the rough and smooth bands separately. As a reference, we compare these local spectra to the spectrum calculated over the entire domain, that is, ignoring the presence of the surface wave heterogeneity caused by the internal wave. This is akin to field data collection, where internal wave-induced surface signatures may be unknowingly buried within data, often time series of surface elevations.
For each time instant, the omnidirectional spectrum E(k) can be written as the integral of the directional spectrum E(k, θ) with the spreading angle θ as The mean slope spectrum is calculated as We choose a domain with a width of 5λ p in the rough band and the smooth band, respectively, to calculate the local spectra, similar to that shown in Figure 2.
In Figure 8, we plot the time-averaged omnidirectional wavenumber spectra in the quasi-stationary time duration for the three representative cases IW-A, IW-B, and IW-C in Table 1. Note that the existence of the peak in the spectrum of the entire domain at the near-zero wavenumber is induced by the internal wave, which is consistent with recent experimental observations (Lenain & Pizzo, 2021). Specifically, the abrupt change around k = 0.5k p is caused by the surface displacement induced by the internal wave at low wavenumbers. It appears only in the result of the entire domain because the lengths of the smooth and rough bands used for calculating the spectra are 5λ p , which cannot resolve the near-zero wavenumber mode. Therefore, there is no abrupt change in the spectra of the smooth and rough bands. In general, the spectral response to the internal wave is constrained to wavenumbers above the peak, k/k p > 1. An exception to this is the case of large amplitude and large density difference (Figures 8c and 8f), where the soliton-induced modulations are strong, k p within the smooth band is shifted to lower wavenumbers and there is an overall loss in energy at all scales. For IW-A (Figures 8a and 8d), we found a slight increase in elevation and slope for k > k p within the rough band, as compared to the smooth band and entire domain (the latter two being effectively equivalent). This phenomenon changes for the median case (Figures 8b and 8e), where the loss of spectral energy within the smooth band is larger (in absolute value) and distributed across all k > k p . For the subrange k p < k < 2k p , it appears that energy is lost within the rough bound, relative to the reference spectrum across the entire domain. This is counter-intuitive, but the difference is quite small and it is difficult to determine if this represents a physical signal or the result of numerical artifacts. Table 2 further quantitatively shows the statistical significance of the difference between the spectra of the smooth/rough band and the entire domain quantified by P-values for cases IW-A, IW-B, and IW-C. The spectral differences in the smooth/rough and the entire domain are more statistically significant in the higher wave number region indicated by the   smaller P-values. Specifically, the spectral differences are significant at a 99% confidence interval for k > 2k p for all three cases.
The mean slope ratio 〈R s 〉 can be used to systematically analyze the impact of internal wave-induced surface currents on surface waves in the rough and smooth bands. This quantity should be correlated with the mean-squared slope, which is known to be important for modeling the aerodynamic surface roughness parameter z 0 (Anctil & Donelan, 1996). Thus, without making any assumptions, 〈R s 〉 is derived directly from the simulated wave fields and should be indicative of the internal wave impacts on the air-water coupled system.
In Figure 9a, we plot the variation of 〈R s 〉 of the cases with varying density ratios for different internal wave amplitude a iw /h u = −0.71, −0.60 and −0.40, and in Figure 9b, we present the trend of 〈R s 〉 for varying internal amplitude with different density ratio ρ u /ρ l = 0.997, 0.990, and 0.975. A linear trend between 〈R s 〉 and ρ u /ρ l for a iw / h u = −0.71 is observed, which follows from the importance of density ratio to determine the maximum surface current magnitude and gradient above the internal wave. For small density difference, that is, ρ u /ρ l = 0.997, 〈R s 〉 is insensitive to varying a iw /h u until a threshold value of between −0.4 and −0.2. Below these values, 〈R s 〉 is effectively a constant ∼2 regardless of the internal wave amplitude. However, with increasing density difference to ρ u /ρ l = 0.990 and 0.975, 〈R s 〉 increases with increasing internal wave amplitude.
It is well known that the horizontal current gradient is central to the surface response (Alpers, 1985;Peregrine, 1976). Here, we analyze the relationship between internal wave properties, that is, ρ u /ρ l and a iw , and the surface current gradient, and its impact on ⟨ ⟩ . As shown in Figures 6b and 6a larger internal wave amplitude and a larger density difference lead to a larger surface current and also a larger surface current gradient. In Figure 10a, we interpolate ⟨ ⟩ for the 14 cases to a contour plot. Note that a strong correlation between ⟨ ⟩ and |dU s /dx| max exists with a correlation coefficient of 0.96 and a P-value of 6.1 × 10 −8 calculated from the 14 cases, as shown in Figure 5. Figure 10b shows the plot 〈R s 〉 as a function of the maximal surface current gradient |dU s / dx| max . In all cases, the values of 〈R s 〉 fall into the range between 1 and 5. The theoretical lower limit of 〈R s 〉 is 1 when the surface wave field is spatially homogeneous in the absence of the internal wave. As expected, a stronger surface current gradient induced by the internal wave leads to a larger surface wave slope variation. The relation between 〈R s 〉 and |dU s /dx| max is approximately linear. Non-representative results from individual simulations (e.g., |dU s /dx| max ∼ 0.006, ⟨ ⟩ ∼ 3.75 ) may be due to random numerical fluctuations, and it is likely we would find better convergence by using ensemble average of simulations using different initial wave fields.

Effect of Varying Surface Wave Field
This series of experiments are designed to test the limit of the internal wave induced surface response with increasing surface wave development. Using the wind speed-dependent JONSWAP spectrum enables a straightforward means of evaluating this limit. A corollary can be drawn to an oceanic domain where the internal wave 14 of 19 forcing is quasi-stationary and the local wind is ramping up (e.g., for the afternoon sea breeze): at what point do we expect these smooth/rough bands to disappear? This scenario is commonly observed and directly reflects qualitative observations made from the FLIP during the 2017 CASPER experiment and the analysis presented in Ortiz-Suslow et al. (2019).

Case Setup
We conduct the numerical experiments with different initial surface wave fields constructed from increasing wind speed U 10 ranging from 0.5 to 5 m s −1 with a fixed fetch of 50 km, as shown in Table 3. The resulting peak wavelengths of surface wave fields range from 5.2 to 24.1 m. The underneath internal wave profile is the same for all the cases with a iw /h u = −0.71 and ρ u /ρ l = 0.997. The same grid numbers N x = 6,144 and N y = 1,536 and the same simulation time step Δt = 0.07 s are applied to all the simulations. The domain length for each case is chosen to satisfy the periodic boundary condition for surface waves and maintain the grid resolution Δ e comparable to each other. For the analysis below, we consider only the time within the quasi-stationary subrange used to calculate 〈R s 〉.

Surface Wave Modulation
To qualitatively view the effect of wind speed on the surface response, as shown in Figure 11, we plot the surface elevation η u (x, y) and surface wave slope S(x, y) for the smallest, medium, and largest wind speed among the cases in Table 3, that is, U 10 = 0.5 m s −1 , 2 m s −1 , and 5 m s −1 . For small wind speeds as U 10 = 0.5 m s −1 , the rough and smooth bands are distinct in the surface elevation and slope fields (Figures 11a and 11b). By 2 m s −1 (Figures 11c  and 11d), the rough and smooth regions are evident in both the elevation and slope fields, with the former being slightly harder to visualize as compared to the very light wind conditions. As expected, for 5 m s −1 (Figures 11e and 11f), the elevation field is essentially homogeneous, however, the internal wave-induced heterogeneity remains evident in the slope field out to a few λ p from the soliton trough.
In Figure 12, we further plot the omnidirectional wavenumber spectra for the cases in Figure 11. Consistent with the observation in Figure 11, the difference in spectral energy for both elevation and slope variance, in the rough and smooth regions, becomes less obvious with increasing U 10 . In Table 4, we further show the statistical significance of the difference between the spectra of the smooth/rough band and the entire domain quantified by P-values for the three cases. Similar to Table 2, the spectral differences between smooth/ rough bands and the entire domain are statistically significant at higher wavenumbers, and are at a 99% confidence interval for k > 2k p for all three cases. Specifically, for the case with U 10 = 5 m s −1 (Figures 12c and 12f), the difference between the spectra in the smooth and rough bands is statistically  Note.
Here Δ e,x = 1.5L x /N x and Δ e,y = 1.5L y /N y are the effective grid sizes after considering dealiasing. In all cases, the grid numbers are N x = 6144 and N y = 1536, and the time step is Δt = 0.07 s. The same fetch F = 50 km is used for all the surface wave fields. 10.1029/2022EA002394 15 of 19 significant at a 99% confidence level but physically very small, which is consistent with the indistinguishable smooth and rough bands in the instantaneous surface wave field as shown in Figure 11e. For the very high wavenumbers, minor differences between local spectra across the smooth and rough bands may be evident. For the median case (2 m s −1 , Figures 12b and 12e), the local spectra diverge for all scales above k p . It is also found that the differences between the rough band and entire domain spectra are larger in magnitude-especially at high kas compared to similar contrasts with the smooth band. For the very light wind scenario (Figures 12a and 12d), we find an even larger divergence in spectral energy for elevation and slope variance across all k > k p . We do not find any shift in k p for the smooth (or rough) bands, as is noted in the analysis of the internal wave-varying-subset of simulations in Section 3.2.2.
In Figure 13, we plot 〈R s 〉 and the peak surface wavelength λ p as a function of the wind speed U 10 . The peak wavelength is calculated using Equation 20, which increases with a larger U 10 . As expected, there is a strong negative relationship between 〈R s 〉, that is, the surface slope ratio between rough and smooth bands, and wind speed. We also found that uncertainty in 〈R s 〉 increases with decreasing wind speed. Kodaira et al. (2016) provide a simplified explanation for the existence of the rough and smooth bands based on the variation of wavenumber caused by the surface current. Specifically, the wavenumber variation Δk is more significant when the original k is larger (i.e., the wavelength is shorter). Our results in Figure 13 are consistent with Kodaira et al. (2016) by showing that the significant surface response is associated with smaller U 10 resulting in a smaller wavelength of surface waves.

Impacts on the Parameterized Aerodynamic Roughness z 0
There are numerous approaches to parameterizing 0 (Equation 1) in coupled model systems. Typically, they take the form of z 0 = 0.11ν/ * + β 0 , where β 0 is the aerodynamic roughness and 0 = 0.011 2 * ∕ is the canonical form following Charnock (1955). In the widely used coupled ocean-atmosphere response experiment algorithm (Edson et al., 2013), the leading coefficient (0.011) is prescribed as wind or surface wave dependent using bulk statistics, for example, wind speed, significant wave height, etc.
Using multiple linear regression on measurements of the directional wave spectra of shoaling waves, Anctil and Donelan (1996) proposed a wave height and steepness dependence: where σ η is the root-mean-square of surface elevation η u (x, y) and σ s is the root-mean-square of surface slope S(x, y), respectively. We utilize this formula for our analysis because it includes higher order wave statistics more sensitive to the internal wave modulations (Figure 12). We found that parameterizing β 0 using lower order statistics (e.g., peak wave age or significant steepness) is insensitive to the internal-surface wave interactions captured by the simulations. Note that there are a variety of choices to parameterize β 0 . Therefore, future studies are needed to systematically examine the accuracy of β 0 related to the different parameterizations, such as those obtained by Janssen (1989Janssen ( , 1991.
In line with our analysis of 〈R s 〉, we compare β 0 in the rough and smooth bands: We observed that the spectral variations induced by the internal wave are concentrated at higher wavenumbers (k/k p > 1), see Figure 12. Therefore, we expect that the root-mean-square surface slope term above should dominate the comparison in Equation 34:   , the root-meansquare of the surface steepness ratio has the highest correlation coefficient of 0.999 with a P-value of 2.55 × 10 −5 . Figure 14 confirms this, showing the dependence of the surface roughness ratio and the root-mean-square surface slope ratio on U 10 . Both ratios have a maximum at the lowest wind (U 10 = 0.5 m s −1 ) and decrease monotonically as the wind increases to U 10 = 5 m s −1 . By 2-3 m s −1 , β 0 in the rough band is still an order of magnitude larger than the smooth band. Interestingly, at the maximum wind, β 0,rough /β 0,smooth ≠ 1. Thus, even when the surface modulation by the internal wave is not obvious in the direct observation of sea surface elevation (Figure 11), there is still significant surface roughness heterogeneity caused by the internal waves.

Conclusions
In this study, we have performed a controlled numerical experiment to evaluate the modulation of surface waves by a propagating nonlinear internal soliton using a deterministic phase-resolved two-layer fluid model with a free surface. This modulation takes the form of the "banding", or alternating regions of smooth and rough ocean surface readily visible on the ocean surface under certain conditions. For the internal wave, we tested the surface response to varying density ratios between the two layers and the soliton amplitude. For the surface waves, we evaluated how a developing sea state responds to a constant internal wave forcing using the JONSWAP spectrum, which has a prescribed wind speed dependence. Our simulations feature the realistic setting of the interaction between a broadband surface wave field and a typical oceanic internal wave.
The surface response induced by the internal wave is featured by the visual observation of a rough band following a smooth band, which can also be reflected in the difference in the local omnidirectional wave energy and slope spectra. We have quantified the surface response by using a systematic approach to evaluate the mean surface slope ratio between the rough and smooth bands, that is, ⟨ ⟩ . Our results show that ⟨ ⟩ is more sensitive to the density difference of the upper and lower layer fluids, while the effect of internal wave amplitude is less significant. Through spectral analysis of the surface wave energy and slope densities, we found that the internal wave primarily impacts wavenumbers above the peak, but can affect waves near the peak in the case of very large density differences between the two layers. In general, we found a strong linear response between ⟨ ⟩ and the maximum surface current gradient, regardless of the specific internal wave parameters.
We investigated the impact of an internal wave on a developing sea by varying the parametric wind speed in the initial surface wave field (JONSWAP-type) from U 10 = 0.5-5 m s −1 . At the upper wind limit, we found that the internal wave modulation of the surface elevation diminishes, but remains weakly evident in the slope field. Consistent with the varying internal wave scenario, the elevation and slope spectra are most impacted by the internal wave at k > k p , and there is a clear negative (nonlinear) dependence in ⟨ ⟩ with developing sea state. We applied our findings to a model of aerodynamic roughness, β 0 , and found that the surface roughness parameter is sensitive to the internal wave modulations only when accounting for wave steepness statistics, that is, mean-square-slope. This approach is not conventionally used in bulk parameterizations and suggests that conventional coupled modeling cannot account for the impact of internal waves, even if the phenomenon is resolved in Figure 13. Mean slope ratio ⟨ ⟩ and surface wavelength λ p as functions of U 10 for cases in Table 3. The error bars represent the range of standard deviations ± 3 within the quasi-stationary time duration.   the model scheme. Furthermore, we found that even at a 5 m s −1 sea statewhen the internal wave-induced banding is invisible in the surface fieldsthe internal wave still induces a ∼2-fold increase in β 0 . At 2-3 m s −1 , the roughness varies by an order of magnitude. These numerical results provide independent corroboration to the hypothesis of Ortiz-Suslow et al. (2019): that the mechanical adjustment of the airflow above internal wave-induced surface roughness heterogeneity is significant enough to locally increase near-surface wind shear and enhance air-sea momentum flux over the internal wave.
Finally, we remark on some limitations of this study. For example, we have neglected the effects of the stratification associated with the continuous density profile, turbulence, thermal effects, wind input, and the stronger nonlinearity in the internal-interface. As a result, the quantitative relationships found here provide suggestions to the natural phenomena, but are not directly translatable. Future studies would benefit from employing a full Navier-Stokes solver and extending our approach to include internal wave packets or collections of individual solitary waves, which may have a cumulative impact on the surface.