Two ocean models converge on coastal submesoscale fronts in idealized tests

Two ocean models built on fundamentally different numerical machinery produce strikingly similar pictures of the small, energetic fronts and eddies that churn coastal waters, according to a new idealized model comparison published in Geoscientific Model Development. The study, led by Kyle Hinson of Pacific Northwest National Laboratory, pitted MPAS-O, an unstructured-grid global ocean model developed for the Department of Energy’s Energy Exascale Earth System Model, against ROMS, a structured-grid regional model long used to study submesoscale dynamics. Across horizontal resolutions spanning 10 kilometers down to 100 meters, the two models represented submesoscale baroclinic instabilities so similarly that the choice of model mattered less than the choice of grid resolution.

The comparison matters because global, kilometer-scale simulations now realistically capture open-ocean submesoscale dynamics, yet whether global models can represent coastal submesoscales remains largely unexplored. Coastal flows are more heterogeneous than their open-ocean counterparts, shaped by complex coastlines and bathymetry, buoyancy inputs from rivers and estuaries, and interactions between coastal and offshore water masses. These factors produce a wider range of dynamical regimes that are harder to parameterize than the relatively homogeneous mixed-layer instabilities typical of the open ocean. The team chose the Mississippi-Atchafalaya River plume system in the northern Gulf of Mexico as a template for initial conditions, a region characterized by strong buoyancy forcing and a summertime, anisotropic, rapidly evolving eddy field with significant near-inertial motions.

The experimental design is deliberately simple. The researchers built a 100-kilometer by 300-kilometer re-entrant channel with a constant depth of 50 meters, free-slip walls along the north and south boundaries, and periodic conditions along the east and west boundaries. The configuration, which the authors call Bichan, a baroclinic instability-inducing channel, runs on an f-plane with an inertial period of roughly 17.4 hours. Despite its simplicity and flat bottom, the setup develops a complex ocean state with multiscale interacting submesoscale fronts and eddies, providing a dynamically rich environment for studying the full lifecycle of frontal eddies, from formation and frontogenesis to the onset of frontal instabilities, arrest, and eventual decay. The authors hope this easy-to-set-up test case, which generates flows highly sensitive to numerical configuration, could become a standard for comparing different models.

Initial conditions were tuned to observed shelf conditions. Temperature and salinity fields were defined by fixed vertical stratification and lateral density gradients, with an initial salinity of 15 grams per kilogram varying only in the cross-channel direction and temperature varying only vertically. The first deformation radius came out at approximately 5 kilometers, similar to conditions observed on the Texas-Louisiana shelf. An along-channel velocity of 0.5 meters per second, relaxed to zero toward the bottom, was seeded with a small 100-kilometer-period perturbation to give instabilities something to grow on. Near-inertial along-channel wind stress, slightly stronger than values observed in the northern Gulf, provided resonant forcing intended to sharpen the fronts and make the setup more similar to summer conditions there. Both models used 50 equally spaced vertical levels, giving 1-meter vertical resolution, so differences between them primarily reflected their numerical schemes and subgrid-scale parameterizations rather than vertical coordinate choice.

The team ran a series of 50-day experiments at resolutions of 10, 5, 2, and 1 kilometer and 500, 200, and 100 meters. Snapshots of surface relative vorticity on day 12 showed submesoscale fronts on the order of 1 to 10 kilometers and eddies of roughly 50 kilometers, about ten times the deformation radius, emerging across all resolutions. At 5 to 10 kilometers, which the authors consider eddy-permitting, the representation of submesoscale frontal processes was limited in both models. At 1 to 2 kilometers, considered eddy-resolving, relative vorticity magnitude increased and braided frontal instabilities emerged in both models. At sub-kilometer resolutions, cross-frontal processes became resolvable and frontal instabilities grew more complex, with the imprint of gravity waves propagating away from frontal convergence zones becoming apparent at 100 to 200 meters in both models.

Timing was one of the clearest differences. Low levels of kinetic energy persisted during the first 4 to 5 days in both models, but around day 5 the mean surface kinetic energy in ROMS rapidly increased, peaking around day 10 at 0.65 square meters per second squared for all sub-kilometer resolutions. The comparable peak in MPAS-O occurred 4 days later, at similar levels. ROMS attained higher surface kinetic energy roughly 4 days before MPAS-O across all resolutions, because baroclinic instabilities across the entire ROMS domain were generated nearly simultaneously, whereas MPAS-O tended to generate instabilities in a more spatially concentrated area near the center of the domain and more slowly at higher resolutions. Despite the timing differences, both models eventually converged to a similar solution over the remainder of the 50-day period, a feature the authors read as encouraging for potential long-term stability of model linkages.

Statistical comparisons relied on flow invariants, quantities independent of coordinate choice, namely total strain, relative vorticity, and divergence, computed as probability distributions over the surface layer during a 10-day window of maximum surface kinetic energy. Inter-model differences were consistently smaller than differences across resolutions within a single model. Coarser resolutions produced fewer fronts and lower levels of divergence and convergence, while both models showed a marked change in divergence distributions as resolution dropped below 1 kilometer. At 100 meters, the percentage of strongly divergent and convergent flows along fronts reached 55 percent in MPAS-O and 50 percent in ROMS. Joint distributions of vorticity against strain and divergence showed similar behavior in both models, with increasing resolution strengthening the relationship between positive vorticity and strain rates along a line of pure shear, and a clear divide at negative vorticity equal to minus the Coriolis parameter remaining visible in both, though more distinctly in ROMS.

Not everything was clean at the finest scales. Grid-scale noise emerged in MPAS-O at sub-kilometer resolutions, organizing into wave-like structures reminiscent of Gibbs phenomenon, much weaker at 500 meters and growing with increasing resolution. The authors attribute this to the biharmonic viscosity, scaled to a grid Reynolds number of 20 that was tuned for coarser-resolution simulations, requiring retuning at sub-kilometer resolutions for the model’s new second-order Adams-Bashforth time-stepping scheme. A simple rerun of the 500-meter case with double the biharmonic viscosity, corresponding to a grid Reynolds number of 10, sufficiently reduced the noise. A systematic retuning was beyond the scope of the study given the computational cost; the 100-meter simulation required nearly six days of wall-clock time on 4096 cores and produced 11 terabytes of output. The authors suggest flow-adaptive horizontal mixing schemes, such as those of Leith and Smagorinsky, as a preferable path forward.

The starkest quantitative divergence lay in mixing. Using an online discrete variance decay diagnostic implemented in MPAS-O for this study and verified against a tracer variance budget in a linear advection test, the team separated numerical mixing, the spurious tracer variance destruction caused by discretizing advection, from physical mixing by vertical diffusion. ROMS generally showed more numerical mixing and less physical mixing than MPAS-O across all resolutions. Bulk numerical mixing decreased by 90 percent in ROMS from 10 kilometers to 100 meters, and by 89 percent in MPAS-O from 10 kilometers to 200 meters. Physical mixing increased with resolution in both models, by 52 percent in MPAS-O from 10 kilometers to 100 meters and by nearly 200 percent in ROMS. At 5 to 10 kilometers, MPAS-O exhibited over 50 percent more physical mixing than ROMS, though differences diminished at finer resolutions. The ratio of bulk numerical to physical mixing fell from 42 percent at 10 kilometers to 2.8 percent at 200 meters in MPAS-O, rising again to 40 percent at 100 meters, while in ROMS it fell from 131 percent to 4 percent.

The authors trace these mixing contrasts to specific schemes. MPAS-O’s flux-corrected transport tracer advection exhibits a stronger variance-producing response than the MPDATA scheme used in the ROMS simulations, and volume-averaged numerical mixing in MPAS-O oscillated between positive and weakly negative values after the frontogenic phase, implying an effective diffusivity that becomes negative, with the strongest variance production in the 100-meter run between days 15 to 20 and 25 to 30. ROMS, by contrast, stayed positive definite at all resolutions. Differences in physical mixing likely arise from the models’ differing configurations of the K-Profile Parameterization, with MPAS-O using the CVMix library and ROMS its internal scheme. The elevated vertical physical salinity mixing in MPAS-O, the authors hypothesize, accelerated frontal arrest and decay, explaining why MPAS-O showed noticeably less frontal structure than ROMS beyond the halfway point of the simulations. The similarities found here, the team concludes, hold promise for future coupling between the two models, and a companion paper extends the comparison to realistic simulations of the Mississippi-Atchafalaya plume.

Subject of Research: An idealized model comparison of coastal submesoscale baroclinic instabilities between the MPAS-O global ocean model and the ROMS regional model across resolutions from 10 kilometers to 100 meters.

Article Title: Towards high-fidelity simulations of coastal submesoscale baroclinic instabilities with MPAS-O (vE3SM3.0.0) – Part 1: Idealized experiments

Article References: Hinson, K. E., Schlichting, D. R., Hetland, R. D., Engwirda, D., Smith, K. M., Petersen, M. R., & Uyeda, K. (2026). Towards high-fidelity simulations of coastal submesoscale baroclinic instabilities with MPAS-O (vE3SM3.0.0) – Part 1: Idealized experiments. Geoscientific Model Development, 19(21), 10321-10344. https://doi.org/10.5194/gmd-19-10321-2026

Image Credits: AI Generated

DOI: 10.5194/gmd-19-10321-2026

Keywords: MPAS-O, ROMS, submesoscale, baroclinic instability, ocean modeling, numerical mixing, frontogenesis, Mississippi-Atchafalaya plume, grid resolution, tracer advection, K-Profile Parameterization, idealized experiments