| 11 Jun 2026
Internal wave–driven diurnal density stratification in the Eastern Arabian Sea
Abstract. The pycnocline and the surface mixed layer are the most prominent indicators of diurnal variability in the ocean. We configured a coupled framework to characterize these processes by combining the Price–Weller–Pinkel model resolving mixed-layer dynamics, with an internal wave model based on the Garrett–Munk spectra accounting for vertical displacements. This coupled system is implemented at three OMNI buoy sites (AD08 (12°N, 68.65°E), AD09 (8.22°N, 73.28°E), and AD10 (10.34°N, 72.58°E)) in the Eastern Arabian Sea. Available in-situ observations validate the simulated hourly density variability, showing good agreement with the observed temporal variability. The coupled model further captures diurnal variability and associated isopycnal displacements over spatial scales of up to 10 km, demonstrating its capability to resolve fine-scale stratified structures. Further, the study presents how sound speed is modulated by the isopycnals displacement. To extend its applicability, the model is initialized using output from the 3D MITgcm during representative months (January, April, July, and October) of different seasons. For shallow-water regions, the modified shallow-water Garrett–Munk spectral configuration is implemented, leading to reduced RMSE in hourly density simulations compared to the deep-water configuration. This MITgcm-PW-GM integrated model demonstrates that the 3D MITgcm model can provide initial vertical density conditions at any location within its domain, enabling localized simulations of diurnal density variability. These results highlight the capability of the coupled–integrated framework to efficiently resolve mixed-layer variability and internal-wave–induced isopycnal displacements associated with diurnal forcing. The ability to resolve diurnal variability of internal wave also makes the framework valuable for regional environmental applications further advancing the understanding of upper-ocean dynamics. In addition, it also supports applications such as simulation of acoustic field and its propagation losses which are crucial for naval operations.
General assessment
The manuscript presents a modelling workflow: coupling a 3D MITgcm background state, a 1D Price–Weller–Pinkel (PWP) mixed-layer model, and a Garrett–Munk (GM) internal-wave field, to reconstruct local sub-daily density and sound-speed variability where real-time profiles are unavailable. The multi-station, multi-season validation and the practical (acoustic/naval) motivation are appreciated. However, the central methodological claim, which MITgcm+PWP+GM can reconstruct a credible diurnal density field at locations without observations is never tested against the appropriate control, and several results appear affected by an amplitude artifact. Specifically: (i) the "diurnal" framing is inconsistent with the GM formulation and with the validation metric; (ii) the internal tide is likely double-counted in the integrated configuration; (iii) reported isopycnal-displacement (ID) amplitudes are physically implausible, exceed the plotted colorbars, and are consistent with an ill-conditioned definition (Eq. 9); (iv) the coupling is purely kinematic, so the "IW-driven" mechanism has no dynamical realization in the model; and (v) key parameters and the core 2D methodology are not documented. I recommend major revision.
Major comments
M1 The "diurnal" framing is inconsistent with the GM spectrum and with the validation metric.
The title and abstract characterize the reconstructed field as "diurnal variability of internal waves" (e.g., L27), but the GM spectrum is a broadband background continuum (f to N) with no diurnal periodicity. Genuine diurnal signals arise only from (i) PWP solar forcing and (ii) diurnal internal-tide lines (K1/O1), which enter through the MITgcm, but not GM. The validation metric, daily standard deviation (L318), mixes the diurnal band with all higher frequency GM variance and cannot isolate a diurnal component. Could authors provide diurnal band harmonic analysis or day composited diagnostics, and define explicitly which internal wave classes the framework represents. As currently formulated the framework represents a linear, broadband, random-phase background continuum, excluding nonlinear solitons and not explicitly carrying internal tide phase, so the "diurnal" terminology should be reconsidered.
M2 Probable double counting of internal tides in the integrated configuration.
The MITgcm background already resolves the internal tide. Superimposing a GM displacement field adds IW variance in the diurnal/semidiurnal bands on top of variance already present. Authors should justify why an external GM field is needed once the internal tide is explicitly resolved, and demonstrate that IW variance is not counted twice.
M3 The isopycnal-displacement definition (Eq. 9) is ill-conditioned below the pycnocline, and the reported amplitudes exceed the plotted colorbars.
The displacement is diagnosed as η = ρ′/(dρ̄/dz) (Eq. 9). This definition is ill-conditioned wherever the background stratification is weak. Below the pycnocline dρ̄/dz close to 0, so the small denominator amplifies any residual density anomaly, including numerical noise, into large apparent displacements. Moreover, in the present construction the density anomaly itself is generated as ρ′ ≈ ζ·(∂ρ_p/∂z) (Eqs. 6–7), which also vanishes as the gradient vanishes; η therefore approaches an indeterminate 0/0 form in the weakly stratified interior, where the result is controlled by noise rather than by a physical displacement.
This conditioning problem is consistent with the reported pattern. The maximum ID at AD09 in January (~34 m) occurs precisely at 200 m (L443), i.e. below the pycnocline in the weakly stratified layer where dρ̄/dz is smallest, exactly where Eq. (9) is most unstable, and where a genuine internal wave displacement would be expected to decay rather than peak. Such amplitudes are also large for background (GM-type) internal waves in the open ocean; while large displacements can occur at internal-tide generation sites, a ~34 m daily displacement in the deep interior warrants independent physical support (local N, GM energy level, or comparison with observed IT amplitudes) rather than being reported at face value.
Separately, and independent of the conditioning issue, the quoted amplitudes are inconsistent with the figures: values of 27–34 m (L443, L448, L454) and ~18 m (L554, L557, L600) exceed the Fig. 5 (±17.5 m) and Fig. 8 (±10 m) colorbars, so these regions are saturated and the reported values cannot be read from the panels.
M4 Purely kinematic coupling: the "IW-driven" mechanism has no dynamical realization.
The coupling is one way and offline: GM displacements are kinematically superimposed on PWP profiles (Eqs. 6–7) with no feedback of IW-induced shear or mixing onto the mixed-layer dynamics. Consequently the framework structurally cannot produce internal-wave-driven mixing, but the title claims an "IW-driven" process. In addition, Eqs. (6)–(7) linearize the displacement (first-order Taylor, T = Tp + ζ ∂Tp/∂z), while for ζ of 20–34 m across a sharp pycnocline this is invalid. Please (a) justify the pure kinematic superposition or discuss its limitations for stratification and acoustics, and (b) apply the finite displacement directly (T(z) = Tp(z − ζ)) or justify the linearization.
M5 Missing GM parameters and undocumented 2D methodology.
The GM reference energy E0, which sets the amplitude of all displacements and variances, and is therefore central to M3 and M7, is never specified, making the results non-reproducible. The shallow-water slope is left as "p = 3 or 4" (L219) without stating the value used. More critically, the stated novelty, a "revised approach for solving the IW boundary-value problem" enabling ID in time and space (L138–139) is not described: how horizontal wavenumbers enter, how random phases are assigned, and how the 2D field ζ(r,z,t) is synthesized are given only by citation (Yadidya et al. 2021; Evans 2000). If the 2D spatial extension is the main contribution, it must be presented explicitly, otherwise the study reads as an incremental application of Yadidya et al. (2021) to a new region. Please also report the mixing coefficients, sponge strength, and bottom-drag coefficient (only the 120 s time step is currently given).
M6 No absolute-density comparison; deterministic claims from a stochastic field.
Only the daily standard deviation is compared, no absolute density/temperature profile, and no single snapshot, is compared against observations. This is not incidental: in the integrated configuration the mean state of any reconstructed profile is essentially the MITgcm/ORAS5 field, and the GM term is a zero-range-mean perturbation, so an absolute comparison would test background-state accuracy, which the std metric removes. Moreover, ζ(r,z,t) is a single random phase realization, so point-by-point agreement with an observed snapshot is not expected by construction, the framework is meaningful only statistically. However, the text reports deterministic, location/time-specific extrema (e.g., "Mumbai … maximum ID of 14 m on 27 January," L551; the per-location, per-range readings of Fig. 8, L549–573). Please either provide absolute profile/snapshot comparisons to substantiate the deterministic claims, or restate all results as statistics (variance, spectra, distributions) and state that the framework does not provide deterministic predictions.
M7 Validation: weak metric, questionable magnitudes, and the wrong baseline.
(a) Weak metric. Standard deviation is the only skill measure and is insensitive to phase, frequency, and vertical structure; two very different time series can share the same std. No error bars and no quantitative per-depth metric (RMSE/correlation) are given for any profile in Figs. 2–4(a–j). If E0 is tuned to match observed std, the agreement is not independent.
(b) Discrepancies understated. Even on this lenient metric, clear mismatches exist, the coupled model departs at 150–200 m, and the observed spread at 50–100 m in Jun–Aug substantially exceeds the modeled envelope. "Good agreement" should be tempered.
(c) Axis scaling and magnitude. The k–m panels are capped at 0.5 kg/m³. A daily std of 0.5 kg/m³ implies day-to-day excursions of order 1.5–2 kg/m³, which is a large fraction of the surface-to-100 m density contrast, also the cap hides values that exceed it (e.g., September spikes in Fig.3k). Please convert the daily std to an equivalent vertical displacement (using local dρ̄/dz) and check consistency with the IDs of Figs. 5/8. The large 5 m std is attributed by the authors to hourly surface-flux forcing (L357–359), i.e., PWP solar forcing rather than IW dynamics, at odds with the "IW-driven" framing.
(d) Wrong baseline (Table 1). Table 1 compares the integrated model only against a bare MITgcm run, which contains no sub-daily variability by construction, a comparison designed to be won. The relevant control, integrated vs. the observation-initialized coupled model, is absent, thus the core claim (MITgcm substituting for observations) is never verified. Only July shows a clear reduction (0.62 to 0.44, still the largest).
M8 Shallow v.s. deep configuration (Fig.7) does not support the claim.
The density RMSE difference is only ~0.014 kg/m³ (~4%, shallow 0.336 v.s. deep 0.350), with interleaving profiles and a ~5-day observational record without error estimates, plausibly within noise. The buoyancy-frequency RMSE (shallow 4.702 / deep 6.246 cph) is comparable to the signal itself, so describing the stratification as "adequately represented" (L517) is not justified for either configuration. There is also a temporal/location mismatch: a May-2009 CTD at Station A is used to validate a framework configured for 2018 at different sites.
M9 Open boundary forcing, volume balance, and spin-up.
Four open boundaries are forced with daily T/S and velocity (L180). Please clarify whether net volume is balanced (inflow = outflow), particularly given BoB intrusion via the West India Coastal Current. With velocity forcing plus ten tidal constituents, is the analysis window sufficient for the internal-tide field and adjustment currents to reach statistical steadiness? Please provide a spin-up diagnostic (e.g., domain-integrated kinetic-energy time series).
M10 Are the four months continuous or independent short integrations?
The MITgcm is configured "separately for four representative months" (L296–299), each with 10 days of spin-up. If these are independent short runs rather than one continuous annual integration, then inter-month differences may conflate genuine seasonal signal with differences in initial conditions and spin-up state. Please clarify, and if independent, explain how seasonal signal is separated from initialization drift.
Minor comments
1.(L39–41) A reference is needed for "IWs are influenced by the structure of the mixed layer." Please clarify whether "structure" refers to mixed-layer depth, the T–S structure, or something else.
2.(L41–45) Please briefly explain the mechanism by which IWs lose energy to the mixed layer or bottom boundary layer under weak stratification.
3.(L51) Mu et al. (2026) attribute ~41% of near-inertial IW energy variability to MLD. Since in this framework the GM field enters almost entirely via MLD/stratification through the PWP profile, please comment on how the remaining ~59% of the driving (wind, topography, mesoscale) is represented.
4.(L80–81) A reference is needed for "sharp stratification contrasts, and energetic mesoscale and IW activity."
6.(L108–133) The literature review is dominated by self-citation and a narrow set of Indian-coast studies; broader international work on PWP+GM coupling and IW acoustics is underrepresented. Please broaden the review to position the paper's specific increment.
7.(L121–125) A reference is needed for the IT’s spatial and temporal variability in the AS. Authors only mentioned the eastern AS.
8.(L179–180) Both initial and boundary conditions derive from the relatively coarse ORAS5, which is later invoked to explain the integrated model's too-deep January MLD (Figs. 2–4). Given that this is a known systematic bias, please justify the choice or provide a sensitivity test rather than attributing the bias post hoc.
9.(L174–176) The eight-grid-point sponge layer is intended to suppress IW reflection. Please discuss whether it also damps the internal tides of interest, given that some are generated at sites close to the boundaries.
10.(L206) Please add references for the stability criteria (bulk/gradient Richardson numbers).
11.(L259–267) the density field is computed using the legacy UNESCO 1983 (EOS-80) formulation (Lines 259–261), whereas the sound speed follows the modern TEOS-10 standard (Line 267). While the numerical discrepancy is negligible for the final results, maintaining consistency within a single thermodynamic framework is recommended.
12.(L285–287) What ensures a smooth join between the PWP (0–500 m) and ORAS5 (>500 m) profiles? A discontinuity in T–S or N² at 500 m would contaminate the GM modal structure at deep sites, since this merged profile feeds the eigenvalue problem.
14.(Figs.2–4 L238) A 2.7 km hydrostatic MITgcm cannot self-resolve the short-wavelength IWs the GM spectrum represents. This motivates the external GM field but also means the MITgcm cannot serve to validate it; please discuss (see also M2).
15.(L292, L366) Please quantify the data exclusions (post-October records; the October 5 and 10 m near-surface data).
16.(L146) Please add the observation time period for the AD08, AD09, and AD10 buoys.
17.(L174-183) Key model parameters are not reported (horizontal and vertical mixing coefficients, sponge-layer strength, bottom-drag coefficient); only the 120 s time step is given.
18.(Eqs.4–5) The two spectra are mutually inconsistent: the deep-water Eq. (4) uses √(ω²+f²) with an ω³ denominator, whereas the shallow-water Eq. (5) uses √(ω²−f²) with ω². The standard band-limited GM spectrum uses ω²−f² and peaks near ω = f; the "+f²" in Eq. (4) removes the near-inertial peak and appears to be a typo. Please also check the 4/π prefactor and the normalization.
19.(Fig.6) Fig. 6 is not an independent result: T, ρ, and c are all generated from the same displacement ζ via Eqs.(6)–(7), so the ID–sound-speed correspondence (L469–474) is an algebraic consequence of the construction rather than a validated finding. Any ID amplitude bias (see M3) propagates directly into the sound-speed excursions.
20.(Figs.5, 8) Colorbar saturation: the plotted ranges (±17.5 m and ±10 m) do not span the amplitudes quoted in the text (see M3).
21.(Fig. 1) The "Maldives" label and pointer target the northern segment of the chain, please reposition toward the central archipelago for a more accurate geographic representation.
22.(L297) Please clarify whether the four monthly simulations were run from four different initial conditions or run from one simulations’ pickup file.
23.Inconsistent spellings throughout: Vellimalai/Vellimallai, Mangalore/Mangaluru, Khambhat/Khambat. Please standardize.
24.(L355) "couple model" should read "coupled model."
25.(L30, L134, L476, L608) The naval-operations motivation is repeated several times, in places nearly verbatim; please consolidate.