Primitive Equation Ocean Model Research Articles

Multistability and intermediate tipping of the Atlantic Ocean circulation.

Tipping points (TP) in climate subsystems are usually thought to occur at a well-defined, critical forcing parameter threshold, via destabilization of the system state by a single, dominant positive feedback. However, coupling to other subsystems, additional feedbacks, and spatial heterogeneity may promote further small-amplitude, abrupt reorganizations of geophysical flows at forcing levels lower than the critical threshold. Using a primitive-equation ocean model, we simulate a collapse of the Atlantic Meridional Overturning Circulation (AMOC) due to increasing glacial melt. Considerably before the collapse, various abrupt, qualitative changes in AMOC variability occur. These intermediate tipping points (ITP) are transitions between multiple stable circulation states. Using 2.75 million years of model simulations, we uncover a very rugged stability landscape featuring parameter regions of up to nine coexisting stable states. The path to an AMOC collapse via a sequence of ITPs depends on the rate of change of the meltwater input. This challenges our ability to predict and define safe limits for TPs.

Subgrid Parameterizations of Ocean Mesoscale Eddies Based on Germano Decomposition

AbstractOcean models at intermediate resolution (1/4°), which partially resolve mesoscale eddies, can be seen as Large eddy simulations of the primitive equations, in which the effect of unresolved eddies must be parameterized. In this work, we propose new subgrid models that are consistent with the physics of two‐dimensional flows. We analyze subgrid fluxes in barotropic decaying turbulence using Germano (1986, https://doi.org/10.1063/1.865568) decomposition. We show that Leonard and Cross stresses are responsible for the enstrophy dissipation, while the Reynolds stress is responsible for additional kinetic energy (KE) backscatter. We utilize these findings to propose a new model, consisting of three parts, that is compared to a baseline dynamic Smagorinsky model. The three‐component model accurately simulates the spectral transfer of energy and enstrophy and improves the representation of KE spectrum, resolved KE and enstrophy decay in a posteriori experiments. The backscattering component of the new model (Reynolds stress) is implemented both in quasi‐geostrophic and primitive equation ocean models and improves statistical characteristics, such as the vertical profile of eddy KE, meridional overturning circulation and cascades of kinetic and potential energy.

Anisotropic diffusion filters for flow-dependent variational data assimilation of sea surface temperature

It is well known that the background error covariance matrix in the data assimilation is anisotropic and varied with the flow field, that is, flow-dependent. However, it is generally assumed to be isotropic and static for ease of implementation. Therefore, to enhance the quality of the background error covariance, two easy-to-implement, flow-dependent three-dimensional variational analysis (3DVar) schemes, the advection–diffusion filter (ADF) and the Smagorinsky diffusion filter (SDF), are tentatively proposed. These serve as alternative schemes for the traditional methods that are based on the isotropic assumption. First, a correctness test of the gradient for the ADF is performed to verify the validity of the adjoint code. Then, the response analysis for the 1D advection–diffusion​ equation is investigated to validate the anisotropic characteristic of the ADF with various flow velocities. Second, the two flow-dependent 3DVar schemes are verified via single-observation assimilation experiments. Compared to the existing diffusion filter with a static diffusion coefficient, both new schemes show a flow-dependent capability in the assimilation experiment because they successfully spread the observational information anisotropically. Finally, a retroactive real-time forecast experiment framework based on a 3D baroclinic primitive equation ocean model is used to investigate the effects of the ADF/SDF on the analysis and forecast accuracy of sea surface temperature (SST). While the traditional variational assimilation method with an isotropic diffusion filter generally reduces initial field errors, both the ADF and SDF can further enhance the accuracy of the initial fields in strong current areas. In addition, better forecast skills can be obtained with the two flow-dependent schemes in the East China Sea.

On divergence- and gradient-preserving coarse-graining for finite volume primitive equation ocean models

We consider the problem of coarse-graining in the context of finite-volume fluid models. If a variable is defined on a high-resolution grid it may be coarse-grained so that it is defined on a grid of lower resolution. In general this will cause some information about the variable to be lost. In particular, horizontal divergences, gradients or other operators calculated on the coarse grid after projecting may differ from those calculated on the fine grid. In some cases we are able to choose averaging weights for coarse-graining such that the coarse-grid operators will give a result approximating that of the corresponding fine-grid operators applied on the fine grid. In this work we derive general conditions on the averaging weights that allow the divergence and gradient to be preserved. These conditions are applied to a regular triangular mesh with B-grid variable placement in which the fine-grid resolution is some integer multiple N of the coarse-grid resolution. For this case we find particular values for the averaging weights that preserve the divergence, and a set of different averaging weights that preserve the gradient. We observe that the vertical component of the curl is also preserved by the same coarse-graining that preserves divergence. These coarse-grainings are applied to data from FESOM2 simulations and we demonstrate that using the particular coarse-grainings derived herein gives an overall reduction in L1 error when compared with other methods.

Snowball Earth Bifurcations in a Fully-Implicit Earth System Model

There is now much geological evidence that the Earth was fully glaciated during several periods in the geological past (about 700[Formula: see text]Myr ago) and attained a so-called Snowball Earth (SBE) state. Additional support for this idea has come from climate models of varying complexity that show transitions to SBE states and undergo hysteresis under changes in solar radiation. In this paper, we apply large-scale bifurcation analyses to a novel, fully-implicit Earth System Model of Intermediate Complexity (I-EMIC) to study SBE transitions. The I-EMIC contains a primitive equation ocean model, a model for atmospheric heat and moisture transport, a sea ice component and formulations for the adjustment of albedo over snow and ice. With the I-EMIC, high-dimensional branches of the SBE bifurcation diagram are obtained through parameter continuation. We are able to identify stable and unstable equilibria and uncover an intricate bifurcation structure associated with the ice-albedo feedback. Moreover, large-scale linear stability analyses are performed near major bifurcations, revealing the spatial nature of destabilizing perturbations.

Effects of Adding Forced Near-Inertial Motion to a Wind-Driven Channel Flow

AbstractAn eddy-resolving primitive equation ocean model is used to examine energy transfers between frequency bands. Steady wind forcing is used to drive a geostrophic channel to which high-frequency winds are added. This excites near-inertial motion, which exerts a Reynolds stress on the slowly varying flow and acts to transfer kinetic energy between low and high frequencies. These transfers extract balanced energy primarily from the mesoscale. A frequency analysis of the transfers shows the bulk of the energy to be extracted from an intermediate range of frequencies that are large relative to the dominant kinetic energy–containing frequencies and small relative to the spectral gap separating high- and low-frequency bands. This phenomenon is robust and is found in systems spanning two orders of magnitude of kinetic energy. Direct calculation of potential energy transfers proved more difficult but nonetheless shows a similar low-to-high frequency transfer. For the parameter range considered, the ratio of potential-to-kinetic energy transfers is slightly larger than unity, and as such is consistent with balanced energy being extracted from horizontal scales that are somewhat larger than the relevant deformation radius.

Exponential time differencing for the tracer equations appearing in primitive equation ocean models

The tracer equations are part of the primitive equations used in ocean modeling and describe the transport of tracers, such as temperature, salinity or chemicals, in the ocean. Depending on the number of tracers considered, several equations may be added to and coupled to the dynamics system. In many relevant situations, the time-step requirements of explicit methods imposed by the transport and mixing in the vertical direction are more restrictive than those for the horizontal, and this may cause the need to use very small time steps if a fully explicit method is employed. To overcome this issue, we propose an exponential time differencing (ETD) solver where the vertical terms (transport and diffusion) are treated with a matrix exponential, whereas the horizontal terms are dealt with in an explicit way. We investigate numerically the computational speed-ups that can be obtained over other semi-implicit methods, and we analyze the advantages of the method in the case of multiple tracers.

Numerical Simulation of the Shallow Water Tides in Bohai Sea

The Princeton Ocean Model (POM), a three-dimensional primitive equation ocean model, is used to reproduce the tidal mechanism in the Bohai Sea. Firstly, the model results are validated by examining the harmonic constants of four major tides at 11 tide gauges, and then the shallow water tides are further focused on. According to the model results, it is found that both M4 and MS4 tides propagate in the similar way, which is characterized by five tidal amphidromic systems in Bohai Sea, and four of them rotate anticlockwise while another one rotates clockwise. It is also found that the M6 tide shows seven amphidromic points in Bohai Sea, among which four rotate anticlockwise and three rotate clockwise. The model results also show that the amplitudes of three shallow water tides are relatively large in shallow coastal waters, which is closely related to generation mechanism of the shallow water tides.

Ocean kinetic energy backscatter parametrizations on unstructured grids: Impact on mesoscale turbulence in a channel

We present a new energy backscatter parametrization for primitive equation ocean models at eddy-permitting resolution, specifically for unstructured grids. Traditional eddy parametrizations in terms of viscosity closures lead to excessive dissipation of kinetic energy when used with eddy-permitting meshes. Implemented into the FESOM2 ocean model, the backscatter parametrization leads to a more realistic total dissipation of kinetic energy. It maintains a reservoir of dissipated energy and reinjects this subgrid energy at larger scales at a controlled rate. The separation between dissipation and backscatter scales is achieved by using different-order differential operators and/or spatial smoothing. This ensures numerical model stability.We perform sensitivity studies with different choices of parameter settings and viscosity schemes in a configuration with a baroclinically unstable flow in a zonally reentrant channel with a horizontally uniform mesh. The best backscatter setup substantially improves eddy-permitting simulations at 1/4° and 1/6° resolution, bringing them close to a 1/12° eddy-resolving reference. Improvements are largest for levels of kinetic energy and variability in temperature and vertical velocity. A selected optimal default scheme is then tested in a mixed resolution setup – a channel with narrow transitions between an eddy-permitting and an eddy-resolving subdomain. The backscatter scheme is able to adapt dynamically to the different resolutions and moves the diagnostics closer to the high resolution reference throughout the domain.Our study is a first step toward using backscatter in global variable-mesh ocean models and suggests potential for substantial improvements of ocean mean state and variability at reduced computational cost.

Lateral Circulation in a Partially Stratified Tidal Inlet

Using a three-dimensional, hydrostatic, primitive-equation ocean model, this study investigates the dynamics of lateral circulation in a partially stratified tidal inlet, the Barataria Pass in the Gulf of Mexico, over a 25.6 h diurnal tidal cycle. Model performance is assessed against observational data. During flood tide, the lateral circulation exhibits the characteristics similar to those induced by differential advection, i.e., lateral flow consists of two counter-rotating cells and is convergent at the surface. The analysis of momentum balance indicates that, in addition to the pressure gradient and vertical stress divergence, nonlinear advection and horizontal stress divergence are also important contributors. During ebb phase, the lateral circulation is mostly toward the right shoal (when looking into the estuary) for the whole water column and persisting for almost the whole period. The surface divergence suggested by the differential advection mechanism lasts for a very short period, if it ever exists. The main momentum balance across most of the transect during ebb is between the along-channel advection of cross-channel momentum and pressure gradient. The sectional averaged lateral velocity magnitude during ebb is comparable to that during flood, which is different from the idealized numerical experiment result.

Components of Upper-Ocean Salt Transport by the Gyres and the Meridional Overturning Circulation

AbstractThe salt transport by the wind-driven gyres and the meridional overturning circulation (MOC) is studied in an idealized-geometry primitive equation ocean model. Two narrow continents, running along meridians, divide the model domain into two basins of different widths connected by a re-entrant channel south of 52.5°S. One of the continents, representing the Americas, is longer than the other, representing Europe/Africa. Two different configurations of the model are used: the “standard” one, in which the short continent is west of the wide basin, and the “exchanged” one, in which the short continent is west of the narrow basin. In both cases, deep water is formed in the basin to the west of the short continent. Most residual transport of the MOC’s upper branch enters this basin by flowing along open streamlines that pass westward south of the short continent before proceeding northward. The meridional salt transport in the upper ocean of the sinking basin is decomposed into two portions: transport along open streamlines and transport by closed streamlines (gyres). In the Northern Hemisphere of the basin in which deep water is formed, the total northward salt transport per unit width along open streamlines is larger in the standard configuration than in the exchanged configuration. This larger salt transport is caused by two factors: a larger northward advection of salt by the interbasin transport and a larger cross-streamline salt transport out of the subpolar gyre. It is concluded that increasing interbasin flow south of Africa would likely bring more salt into the Atlantic Ocean.

Pattern of vertical velocity in the Lofoten vortex (the Norwegian Sea)

Mean radial distributions of various dynamic characteristics of the permanently existing anticyclonic Lofoten vortex (LV) in the Norwegian Sea are obtained from an eddy-permitting regional hydrodynamic MIT general circulation model. It is shown that the model adequately reproduces the observed 3D thermohaline and dynamic structure of the vortex. The obtained radial distribution of the mean vertical velocity is found to form a complex structure: with the upward fluxes along the axis in and above the anticyclonically rotating LV core, compensated by the downward fluxes in the vortex skirt. These vertical motions maintain the vortex potential energy anomaly against dissipation. This secondary circulation is generated by the centrifugal force and, to a lesser extent, by the horizontal dispersion of the vortex energy, both intensified towards the sea surface. Below the vortex core, the maximum downward vertical velocity converges towards the vortex axis with depth. At these depth levels, the secondary circulation is forced by Ekman divergence in the bottom mixed layer. The theory of columnar vortices with helical structure, applied to the LV, relate the radial profiles of the vertical velocity with those of the horizontal circulation. The theoretically predicted the radial patterns of the mean vertical velocity in the LV were close to those, obtained from the primitive equation ocean model, when approximating the radial patterns of the azimuthal velocity with the Rayleigh profile.

Autonomous Sampling of Ocean Submesoscale Fronts with Ocean Gliders and Numerical Model Forecasting

ABSTRACTSubmesoscale fronts arising from mesoscale stirring are ubiquitous in the ocean and have a strong impact on upper-ocean dynamics. This work presents a method for optimizing the sampling of ocean fronts with autonomous vehicles at meso- and submesoscales, based on a combination of numerical forecast and autonomous planning. This method uses a 48-h forecast from a real-time high-resolution data-assimilative primitive equation ocean model, feature detection techniques, and a planner that controls the observing platform. The method is tested in Monterey Bay, off the coast of California, during a 9-day experiment focused on sampling subsurface thermohaline-compensated structures using a Seaglider as the ocean observing platform. Based on model estimations, the sampling “gain,” defined as the magnitude of isopycnal tracer variability sampled, is 50% larger in the feature-chasing case with respect to a non-feature-tracking scenario. The ability of the model to reproduce, in space and time, thermohaline submesoscale features is evaluated by quantitatively comparing the model and glider results. The model reproduces the vertical (~50–200 m thick) and lateral (~5–20 km) scales of subsurface subducting fronts and near-bottom features observed in the glider data. The differences between model and glider data are, in part, attributed to the selected glider optimal interpolation parameters and to uncertainties in the forecasting of the location of the structures. This method can be exported to any place in the ocean where high-resolution data-assimilative model output is available, and it allows for the incorporation of multiple observing platforms.

A potential method to accelerate spin up of turbulent ocean models

We demonstrate a simple method to quickly spin up (equilibrate) the temperature field of a turbulent idealised primitive equation ocean gyre model. We make use of the assumption that both velocity fields and non-linear eddy advection terms equilibrate on a shorter timescale than tracer fields.

Generation of Upwelling Circulation under Downwelling-Favorable Wind within Bottom-Attached, Buoyant Coastal Currents

AbstractA two-dimensional modeling study by Moffat and Lentz recently reported that downwelling-favorable wind can induce cross-shore upwelling circulation within a bottom-attached, buoyant coastal current. Here, we extend the problem to three dimensions. The driving mechanism and the sensitivity for the upwelling circulation are studied, using a primitive equation ocean model and an analytical model. After the initial downwelling adjustment that steepens the isopycnals and compresses the coastal current, the cross-shore flow can switch to steady upwelling circulation. This reverse circulation coincides with a vertically well-mixed water column and persists until interrupted by the arrival of river plume bulge from upstream. During the upwelling phase, the ageostrophic cross-shore flow follows the Ekman balance. The sense of cross-shore circulation is governed by a dimensionless parameter, the shear ratio, which measures the relative size of geostrophic shear and velocity shear supported by the wind in the shallow-water limit. Upwelling circulation occurs when the shear ratio is greater than one. This condition represents that, near the surface, the wind-intensified pressure gradient exceeds the maximum possible Coriolis force associated with the wind-forced, alongshore flow. The resulting upwelling circulation acts to slump the isopycnals to restore the geostrophic balance. Therefore, within a coastal current, decreasing wind stress in fact strengthens the upwelling circulation, as a weaker wind produces a weaker shear and thus increases the imbalance. This inverse relation holds until the wind is too weak to mix the water column. Based on the analytical model, a regime classification for the cross-shore circulation under downwelling-favorable wind is proposed. An observational example is given.

Compact symmetric Poisson equation discretization for non-hydrostatic sigma coordinates ocean model

In anticipation of relaxing the hydrostatic assumption in a sigma coordinates primitive equations ocean model, we show how a projection method can be designed with the use of a compact symmetric 15-point stencil for the Poisson equation. This is achieved by recognizing that, owing the non-orthogonality of the grid, the velocity has a contravariant and covariant set of components. The two sets play a different role in the primitive equations: the contravariant components enter the definition of the model fluxes, whereas the covariant components experience the forces and in particular the pressure gradient. By treating these two sets separately, the discretized gradient and divergence operators are simple finite differences. The two sets of components are related via a linear transformation, the metric tensor, which is entirely determined by the kinetic energy. We show how the spatial discretization of the kinetic energy fully controls the Poisson equation discretization, including its boundary conditions. The discretization of the Poisson equation is shown to converge at second order and to behave as well or better than alternative methods. This approach is a prerequisite to implement later an efficient Poisson solver, such as a multigrid algorithm.

Numerical study on the eddy–mean flow interaction between a cyclonic eddy and Kuroshio

A three-dimensional primitive equations ocean model (POM) is employed to study the momentum and energy balance of a moving cyclonic eddy (CE) during eddy–mean flow interaction. The CE generated by an idealized typhoon forms to the east of the Philippine islands. A momentum balance analysis shows that the dynamics of the CE are generally dominated by the geostrophic current throughout the life cycle of the CE. An energy analysis suggests that the eddy kinetic energy (EKE) and the eddy potential energy (EPE) decay rapidly after generation. The maximum EPE initially appears at the surface of the eddy center and gradually appears in the subsurface layer. The largest baroclinic instability (BCI) initially occurs at the surface. For a CE moving along a trajectory, the conversion from mean potential energy (MPE) to the EPE is positive (negative) in the front (rear) part of the trajectory, and then the eddy transfers its EPE forward along its trajectory by means of the front (rear) part of the eddy obtaining (losing) EPE from (to) the mean flow. During the interaction stage, the northward flowing Kuroshio interacts with the southward flow on the western side of the eddy and the inverse velocity shear between the Kuroshio and the eddy causes the EKE to gradually develop east–west asymmetry. The largest barotropic instability (BTI) is found in the interaction zone. Advection term, pressure work, and friction term play the dominating role in eddy decay in the eddy zone, while BTI only dominates in the interaction zone.

Improve the Simulations of Near-Inertial Internal Waves in the Ocean General Circulation Models

AbstractThe near-inertial wind work and near-inertial internal waves (NIWs) in the ocean have been extensively studied using ocean general circulation models (OGCMs) forced by 6-hourly winds or wind stress obtained from atmospheric reanalysis data. However, the OGCMs interpolate the reanalysis winds or wind stress linearly onto each time step, which partially filters out the wind stress variance in the near-inertial band. In this study, the influence of the linear interpolation on the near-inertial wind work and NIWs is quantified using an eddy-resolving (°) primitive equation ocean model. In addition, a new interpolation method is proposed—the sinc-function interpolation—that overcomes the shortages of the linear interpolation.It is found that the linear interpolation of 6-hourly winds significantly underestimates the near-inertial wind work and NIWs at the midlatitudes. The underestimation of the near-inertial wind work and near-inertial kinetic energy is proportional to the loss of near-inertial wind stress variance due to the linear interpolation. This further weakens the diapycnal mixing in the ocean due to the reduced near-inertial shear variance. Compared to the linear interpolation, the sinc-function interpolation retains all the wind stress variance in the near-inertial band and yields correct magnitudes for the near-inertial wind work and NIWs at the midlatitudes.

Energy budget-based backscatter in an eddy permitting primitive equation model

Increasing computational resources are starting to allow global ocean simulations at so-called “eddy-permitting” resolutions, at which the largest mesoscale eddies can be resolved explicitly. However, an adequate parameterization of the interactions with the unresolved part of the eddy energy spectrum remains crucial. Hyperviscous closures, which are commonly applied in eddy-permitting ocean models, cause spurious energy dissipation at these resolutions, leading to low levels of eddy kinetic energy (EKE) and weak eddy induced transports. It has recently been proposed to counteract the spurious energy dissipation of hyperviscous closures by an additional forcing term, which represents “backscatter” of energy from the un-resolved scales to the resolved scales. This study proposes a parameterization of energy backscatter based on an explicit sub-grid EKE budget. Energy dissipated by hyperviscosity acting on the resolved flow is added to the sub-grid EKE, while a backscatter term transfers energy back from the sub-grid EKE to the resolved flow. The backscatter term is formulated deterministically via a negative viscosity, which returns energy at somewhat larger scales than the hyperviscous dissipation, thus ensuring dissipation of enstrophy. The parameterization is tested in an idealized configuration of a primitive equation ocean model, and is shown to significantly improve the solutions of simulations at typical eddy-permitting resolutions.

Data assimilation experiments using diffusive back-and-forth nudging for the NEMO ocean model

Abstract. The diffusive back-and-forth nudging (DBFN) is an easy-to-implement iterative data assimilation method based on the well-known nudging method. It consists of a sequence of forward and backward model integrations, within a given time window, both of them using a feedback term to the observations. Therefore, in the DBFN, the nudging asymptotic behaviour is translated into an infinite number of iterations within a bounded time domain. In this method, the backward integration is carried out thanks to what is called backward model, which is basically the forward model with reversed time step sign. To maintain numeral stability, the diffusion terms also have their sign reversed, giving a diffusive character to the algorithm. In this article the DBFN performance to control a primitive equation ocean model is investigated. In this kind of model non-resolved scales are modelled by diffusion operators which dissipate energy that cascade from large to small scales. Thus, in this article, the DBFN approximations and their consequences for the data assimilation system set-up are analysed. Our main result is that the DBFN may provide results which are comparable to those produced by a 4Dvar implementation with a much simpler implementation and a shorter CPU time for convergence. The conducted sensitivity tests show that the 4Dvar profits of long assimilation windows to propagate surface information downwards, and that for the DBFN, it is worth using short assimilation windows to reduce the impact of diffusion-induced errors. Moreover, the DBFN is less sensitive to the first guess than the 4Dvar.