Fast dynamical spin-up of ocean general circulation models using Newton–Krylov methods

Numerical models of the ocean play an important role in efforts to understand past climate variability and predict future climate changes. In many studies, ocean models are driven by forcings that are either time-independent or vary periodically (seasonally) and it is often highly desirable or even essential to obtain equilibrium solutions of the model. Existing methods, based on the simple, expedient idea of integrating the model until the transients have died out, are too expensive to use routinely because the ocean takes several thousand years to equilibrate. Here, we present a novel approach for efficiently computing equilibrium solutions of ocean models. Our general approach is to formulate the problem as a large system of nonlinear algebraic equations to be solved with a class of methods known as matrix-free Newton–Krylov, a combination of Newton-type methods for superlinearly convergent solution of nonlinear equations, and Krylov subspace methods for solving the Newton correction equations. As an initial demonstration of the feasibility of this approach, we apply it to find the equilibrium solutions of a quasi-geostrophic ocean model for both steady forcing and seasonally-varying forcing. We show that the matrix-free Newton–Krylov method converges to the solutions obtained by direct time integration of the model, but at a computational cost that is between 10 and 100 times smaller than direct integration. A key advantage of our approach is that it can be applied to any existing time-stepping code, including ocean general circulation models and biogeochemical models. However, effective preconditioning of the linear equations to be solved during the Newton iteration remains a challenge.     

Bifurcation analysis of wind-driven flows with MOM4

In this paper, the methodology of bifurcation analysis is applied to the explicit time-stepping ocean model MOM4 using a Jacobian–Free Newton–Krylov (JFNK) approach. We in detail present the implementation of the JFNK method in MOM4 but restrict the preconditioning technique to the case for which the density distribution is prescribed. For a prescribed density field case, we present bifurcation diagrams, for the first time in MOM4, for the wind-driven ocean circulation. In addition, we show that the JFNK method can reduce the spin-up time to a steady equilibrium in MOM4 considerably if an accurate solution is required.     

Fast spin up of Ocean biogeochemical models using matrix-free Newton–Krylov

A novel computational approach is introduced for the efficient computation of equilibrium solutions of seasonally forced ocean biogeochemical models. The essential idea is to formulate the problem as a large system of nonlinear algebraic equations to be solved with a class of methods known as matrix-free Newton–Krylov (MFNK). MFNK is a combination of Newton-type methods for superlinearly convergent solution of nonlinear equations, and Krylov subspace methods for solving the Newton correction equations. The basic link between the two methods is the Jacobian-vector product, which may be probed approximately without forming and storing the elements of the true Jacobian. To render this approach practical for global models with O(10⁶) degrees of freedom, a flexible preconditioning strategy is developed. The result is an essentially “black-box” numerical scheme than can be applied to most existing biogeochemical models. The method is illustrated by applying it to find the equilibrium solutions of two realistic biogeochemical problems. Compared with the conventional approach of direct time integration, the preconditioned-MFNK scheme is shown to be roughly two orders of magnitude more efficient. Several potential refinements of the basic algorithm that may yield further performance gains are discussed. The numerical scheme described here addresses a fundamental challenge to using ocean biogeochemical models more effectively.     

A spectral barotropic model of the wind-driven world ocean

A global spectral barotropic ocean model is introduced to describe the depth-averaged flow. The equations are based on vorticity and divergence (instead of horizontal momentum); continents exert a nearly infinite drag on the fluid. The coding follows that of spectral atmospheric general circulation models using triangular truncation and implicit time integration to provide a first step for seamless coupling to spectral atmospheric global circulation models and an efficient method for filtering of ocean wave dynamics. Five experiments demonstrate the model performance: (i) Bounded by an idealized basin geometry and driven by a zonally uniform wind stress, the ocean circulation shows close similarity with Munk’s analytical solution. (ii) With a real land–sea mask the model is capable of reproducing the spin-up, location and magnitudes of depth-averaged barotropic ocean currents. (iii) The ocean wave-dynamics of equatorial waves, excited by a height perturbation at the equator, shows wave dispersion and reflection at eastern and western coastal boundaries. (iv) The model reproduces propagation times of observed surface gravity waves in the Pacific with real bathymetry. (v) Advection of tracers can be simulated reasonably by the spectral method or a semi-Langrangian transport scheme. This spectral barotropic model may serve as a first step towards an intermediate complexity spectral atmosphere–ocean model for studying atmosphere–ocean interactions in idealized setups and long term climate variability beyond millennia.     

Bifurcation structure of a wind-driven shallow water model with layer-outcropping

The steady state bifurcation structure of the double-gyre wind-driven ocean circulation is examined in a shallow water model where the upper layer is allowed to outcrop at the sea surface. In addition to the classical jet-up and jet-down multiple equilibria, we find a new regime in which one of the equilibrium solutions has a large outcropping region in the subpolar gyre. Time dependent simulations show that the outcropping solution equilibrates to a stable periodic orbit with a period of 8 months. Co-existing with the periodic solution is a stable steady state solution without outcropping.A numerical scheme that has the unique advantage of being differentiable while still allowing layers to outcrop at the sea surface is used for the analysis. In contrast, standard schemes for solving layered models with outcropping are non-differentiable and have an ill-defined Jacobian making them unsuitable for solution using Newton’s method. As such, our new scheme expands the applicability of numerical bifurcation techniques to an important class of ocean models whose bifurcation structure had hitherto remained unexplored.     

Parametric vertical coordinate formulation for multiscale, Boussinesq, and non-Boussinesq ocean modeling

Two physical parameters are introduced into the basic ocean equations to generalize numerical ocean models for various vertical coordinate systems and their hybrid features. The two parameters are formulated by combining three techniques: the arbitrary vertical coordinate system of Kasahara [Kasahara, A., 1974. Various vertical coordinate systems used for numerical weather prediction. Mon. Weather Rev. 102, 509–522], the Jacobian pressure gradient formulation of Song [Song, Y.T., 1998. A general pressure gradient formation for ocean models. Part I: Scheme design and diagnostic analysis. Mon. Weather Rev. 126 (12), 3213–3230], and a newly introduced parametric function that permits both Boussinesq (volume-conserving) and non-Boussinesq (mass-conserving) conditions. Based on this new formulation, a generalized modeling approach is proposed. Several representative oceanographic problems with different scales and characteristics––coastal canyon, seamount topography, non-Boussinesq Pacific Ocean with nested eastern Tropics, and a global ocean model––have been used to demonstrate the model’s capabilities for multiscale applications. The inclusion of non-Boussinesq physics in the topography-following ocean model does not incur computational expense, but more faithfully represents satellite-observed ocean-bottom-pressure data. Such a generalized modeling approach is expected to benefit oceanographers in solving multiscale ocean-related problems by using various coordinate systems on the same numerical platform.     

Boussinesq-type modelling using an unstructured finite element technique

A model for solving the two-dimensional enhanced Boussinesq equations is presented. The model equations are discretised in space using an unstructured finite element technique. The standard Galerkin method with mixed interpolation is applied. The time discretisation is performed using an explicit three-step Taylor–Galerkin method. The model is extended to the surf and swash zone by inclusion of wave breaking and a moving boundary at the shoreline. Breaking is treated by an existing surface roller model, but a new procedure for the detection of the roller thickness is devised. The model is verified using four test cases and the results are compared with experimental data and results from an existing finite difference Boussinesq model.     

Pacific subtropical cell variability in coupled climate model simulations of the late 19th–20th century

Observational studies of the Pacific basin since the 1950s have demonstrated that a decrease (increase) in tropical Pacific sea surface temperatures (SSTs) is significantly correlated with a spin-up (slow-down) of the Pacific Subtropical Cells (STCs). STCs are shallow wind-driven overturning circulations that provide a pathway by which extratropical atmospheric variability can impact the equatorial Pacific thermocline and, through upwelling in the eastern equatorial Pacific, tropical Pacific SSTs. Recent studies have shown that this observed relationship between SSTs and STCs is absent in coupled climate model simulations of the late 19th–20th centuries. In this paper we investigate what causes this relationship to breakdown and to what extent this limits the models’ ability to simulate observed climate change in the equatorial Pacific since the late 19th century. To provide insight into these questions we first show that the NCAR Community Climate System Model’s simulation of observed climate change since the 1970s has a robust signal in the equatorial Pacific that bears a close resemblance to observations. Strikingly, absent is a robust signal in the equatorial thermocline. Our results suggest that the coupled model may be reproducing the observed local ocean response to changes in forcing but inadequately reproducing the remote STC-forcing of the tropical Pacific due to the underestimate of extratropical winds that force these ocean circulations. These conclusions are found to be valid in five different coupled climate model simulations of the late 19th–20th centuries (CCSM3, GISS EH, GFDL CM2.1, CSIRO-Mk3, and HadCM3).     

On the solution of the carbonate chemistry system in ocean biogeochemistry models

We present a simplified method for solving the local equilibrium carbonate chemistry in numerical ocean biogeochemistry models. Compared to the methods typically used, the scheme is fast, efficient and compact. The accuracy of the solution is dictated by the number of species retained in the expression for alkalinity and there is almost no computational penalty for retaining minor contributions. We demonstrate that this scheme accurately reproduces the results of the commonly used method in the context of a three-dimensional global ocean carbon cycle model. Using this model we also show that neglecting the regional variations in surface dissolved inorganic phosphorus and silicic acid concentrations can lead to significant systematic bias in regional estimates of air–sea carbon fluxes using such models.     

Euler–Lagrange equations for the spectral element shallow water system

We present the derivation of the discrete Euler–Lagrange equations for an inverse spectral element ocean model based on the shallow water equations. We show that the discrete Euler–Lagrange equations can be obtained from the continuous Euler–Lagrange equations by using a correct combination of the weak and the strong forms of derivatives in the Galerkin integrals, and by changing the order with which elemental assembly and mass averaging are applied in the forward and in the adjoint systems. Our derivation can be extended to obtain an adjoint for any Galerkin finite element and spectral element system.We begin the derivations using a linear wave equation in one dimension. We then apply our technique to a two-dimensional shallow water ocean model and test it on a classic double-gyre problem. The spectral element forward and adjoint ocean models can be used in a variety of inverse applications, ranging from traditional data assimilation and parameter estimation, to the less traditional model sensitivity and stability analyses, and ensemble prediction. Here the Euler–Lagrange equations are solved by an indirect representer algorithm.     

Distinguishing coupled ocean–atmosphere interactions from background noise in the North Pacific

When considering physical mechanisms for decadal-timescale climate variability in the North Pacific, it is useful to describe in detail the expected response of the ocean to the chaotic atmospheric forcing. The expected response to this white-noise forcing includes strongly enhanced power in the decadal frequency band relative to higher frequencies, pronounced changes in basin-wide climate that resemble regime shifts, preferred patterns of spatial variability, and a depth-dependent profile that includes variability with a standard deviation of 0.2–0.4°C over the top 50–100 m. Weak spectral peaks are also possible, given ocean dynamics. Detecting coupled ocean–atmosphere modes of variability in the real climate system is difficult against the spectral and spatial structure of this ‘null-hypothesis’ of how the ocean and atmosphere interact, especially given the impossibility of experimentally decoupling the ocean from the atmosphere. Turning to coupled ocean–atmosphere models to address this question, a method for identifying coupled modes by using models of increasing physical complexity is illustrated. It is found that a coupled ocean–atmosphere mode accounts for enhanced variability with a time scale of 20 years/cycle in the Kuroshio extension region of the model's North Pacific. The observed Pacific Decadal Oscillation (PDO) has many similarities to the expected noise-forced response and few similarities to the model's coupled ocean–atmosphere variability. However, model deficiencies and some analyses of observations by other workers indicate that the possibility that part of the PDO arises from a coupled ocean–atmosphere mode cannot be ruled out.     

Sensitivity of AUV response to variations in hydrodynamic parameters

The sensitivity of the response of a typical AUV to changes in hydrodynamic parameters is examined. The analysis is primarily performed using a computer model of an axi-symmetric vehicle typical of many AUVs in service today. The vehicle used is the Canadian Self-Contained Off-the-shelf Underwater Testbed (C-SCOUT), designed and built by graduate and work term students. The fully nonlinear computer model is based on Newton–Euler equations of motion, and uses the component build-up method to describe the excitation forces. The hydrodynamic parameters are varied in a series of simulations with the computer model; the response being analyzed for specific performance indicators.     

A numerical method for the two layer shallow water equations with dry states

A numerical method is proposed for solving the two layer shallow water equations with variable bathymetry in one dimension based on high-resolution f-wave-propagation finite volume methods. The method splits the jump in the fluxes and source terms into waves propagating away from each grid cell interface. It addresses the required determination of the system’s eigenstructure and a scheme for evaluating the flux and source terms. It also handles dry states in the system where the bottom layer depth becomes zero, utilizing existing methods for the single layer solution and handling single layer dry states that can exist independently. Sample results are shown illustrating the method and its handling of dry states including an idealized ocean setting.     

An analytical approach for the calculation of random wave forces on submerged tunnels

This paper deals with the random forces produced by high ocean waves on submerged horizontal circular cylinders. Arena [Arena F, Interaction between long-crested random waves and a submerged horizontal cylinder. Phys Fluids 2006;18(7):1–9 (paper 076602)] obtained the analytical solution of the random wave field for two dimensional waves by extending the classical Ogilvie solution [Ogilvie TF, First- and second-order forces on a cylinder submerged under a free surface. J Fluid Mech 1963;16:451–472; Arena F, Note on a paper by Ogilvie: The interaction between waves and a submerged horizontal cylinder. J Fluid Mech 1999;394:355–356] to the case of random waves. In this paper, the wave force acting on the cylinder is investigated and the Froude Krylov force [Sarpkaya T, Isaacson M, Mechanics of wave forces on offshore structures, Van Nostrand Reinhold Co.; 1981], on the ideal water cylinder, is calculated from the random incident wave field. Both forces represent a Gaussian random process of time. The diffraction coefficient of the wave force is obtained as quotient between the standard deviations of the force on the solid cylinder and of the Froude Krylov force. It is found that the diffraction coefficient of the horizontal force C_do is equal to the C_dv of the vertical force. Finally, it is shown that, since a very large wave force occurs on the cylinder, it may be calculated, in time domain, starting from the Froude Krylov force. It is then shown that this result is due to the fact that the frequency spectrum of the force acting on the cylinder is nearly identical to that of the Froude–Krylov force.     

A spectral barotropic model of the wind-driven world ocean

A global spectral barotropic ocean model is introduced to describe the depth-averaged flow. The equations are based on vorticity and divergence (instead of horizontal momentum); continents exert a nearly infinite drag on the fluid. The coding follows that of spectral atmospheric general circulation models using triangular truncation and implicit time integration to provide a first step for seamless coupling to spectral atmospheric global circulation models and an efficient method for filtering of ocean wave dynamics. Five experiments demonstrate the model performance: (i) Bounded by an idealized basin geometry and driven by a zonally uniform wind stress, the ocean circulation shows close similarity with Munk’s analytical solution. (ii) With a real land–sea mask the model is capable of reproducing the spin-up, location and magnitudes of depth-averaged barotropic ocean currents. (iii) The ocean wave-dynamics of equatorial waves, excited by a height perturbation at the equator, shows wave dispersion and reflection at eastern and western coastal boundaries. (iv) The model reproduces propagation times of observed surface gravity waves in the Pacific with real bathymetry. (v) Advection of tracers can be simulated reasonably by the spectral method or a semi-Langrangian transport scheme. This spectral barotropic model may serve as a first step towards an intermediate complexity spectral atmosphere–ocean model for studying atmosphere–ocean interactions in idealized setups and long term climate variability beyond millennia.     

Free vibrations of three-dimensional extensible marine cables with specified top tension via a variational method

This paper presents a model formulation that can be used for analyzing the three-dimensional vibration behaviours of an inclined extensible marine cable. The virtual work-energy functional, which involves strain energy due to axial stretching of the cable and virtual work done by external hydrostatic forces is formulated. The coupled equations of motion in the Cartesian coordinates of global systems are obtained by taking into account the difference between Euler’s equations and equilibrium equations. The method of Galerkin finite element is used to obtain the mass and stiffness matrices which are transformed into the local coordinate systems. Then the eigenvalue problem is solved to determine its natural frequencies and corresponding mode shapes. The model formulation developed herein is conveniently applied for the cases of specified top tension. The numerical investigations are carried out to demonstrate the validity of the model and to explore in details the influence of various parameters on the behaviours of marine cables. Results for the frequency avoidance phenomenon, maximum dynamic tension and coupled transverse mode shapes are presented and discussed.     

A tree-based solver for adaptive ocean modelling

The development of an adaptive (in space and time) ocean model from an existing adaptive finite-volume Navier–Stokes model is described. A flexible and efficient quadtree spatial discretisation is used which requires collocation of all variables (i.e. an A-grid discretisation). We demonstrate that the use of an approximate projection method allows for implicit damping of instabilities generally associated with the A-grid, at the expense of a relatively small amount of numerical energy dissipation, while accurately preserving dispersive properties and geostrophic balance. The finite-volume formulation also maintains second-order spatial accuracy at all solid boundaries. Test cases demonstrate the efficacy of the adaptive ocean model, and the advantages it has in terms of efficient representation of multi-scale behaviour within a single model. The model is freely available as open-source code.     

A 3-D model based on Navier–Stokes equations for regular and irregular water wave propagation

A spatial fixed σ-coordinate is used to transform the Navier–Stokes equations from the sea bed to the still water level. In the fixed σ-coordinate system only a very small number of vertical grid points are required for the numerical model. The time step for using the spatial fixed σ-coordinate is efficiently larger than that of using a time dependent σ-coordinate, as there is substantial truncation error involved in the time dependent σ-coordinate transformation. There is no need to carry out the σ-coordinate transformation at each time step, which can reduce computational times. It is important that wave breaking can be potentially modeled in the fixed σ-coordinate system, but in a time-dependent σ-coordinate system the wave breaking cannot be modeled. A projection method is used to separate advection and diffusion terms from the pressure terms in Navier–Stokes equations. The pressure variable is further separated into hydrostatic and hydrodynamic pressures so that the computer rounding errors can be largely avoided. In order to reduce computational time of solving the hydrodynamic pressure equation, at every time step the initial pressure is extrapolated in time domain using computed pressures from previous time steps, and then corrected in spatial domain using a multigrid method. For each time step, only a few of iterations (typically six iterations) are required for solving the pressure equation. The model is tested against available experimental data for regular and irregular waves and good agreement between calculation results and the measured data has been achieved.     

Local and far-field surface elevation around a vertical cylinder in unidirectional steep wave groups

The problem of diffraction of a unidirectional incident wave group by a bottom-seated cylinder is considered. We assume the amplitude of the incoming wave to be small in comparison with other linear scales of the problem, and develop the corresponding second-order perturbation theory. We use the Fourier transform to treat time variation and separate spatial variables when solving the non-homogeneous second-order problem. The resulting set of non-homogeneous Bessel equations is solved numerically.Solutions for various types of incoming wave spectrum are obtained including the Gaussian spectrum and the Pierson–Moskowitz spectrum. To validate the method, problems with gradually decreasing bandwidth of Gaussian spectrum are solved and it is shown that the corresponding solution approaches that for the monochromatic case. The Pierson–Moskowitz spectrum with a set of realistic physical parameters is used as an example of extreme wave interaction with an offshore structure. The corresponding first- and second-order solutions are obtained and the effect of non-linearity on the solution is discussed with the emphasis on the growth of maximum free-surface elevation on the cylinder’s surface and generation of high frequency free radiated waves.     

Parameter identification of a compliant nonlinear SDOF system in random ocean waves by reverse MISO method

The determination of the drag and inertia coefficients, which enter into the wave force model given by Morison's equation, is particularly uncertain and difficult when a linear spectral model is used for ocean waves, and the structure is compliant and has nonlinear dynamic response. In this paper, a nonlinear System Identification method, called Reverse Multiple Inputs–Single Output (R–MISO) is applied to identify the hydrodynamic coefficients as well as the nonlinear stiffness parameter for a compliant single-degree-of-freedom system. Four different types of problems have been identified for use in various situations and the R–MISO has been applied to all of them. One of the problems requires iterative solution strategy to identify the parameters. The method has been found to be efficient in predicting the parameters with reasonable accuracy and has the potential for use in the laboratory experiments on compliant nonlinear offshore systems.