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 method to reduce the spin-up time of ocean models

The spin-up timescale in large-scale ocean models, i.e., the time it takes to reach an equilibrium state, is determined by the slow processes in the deep ocean and is usually in the order of a few thousand years. As these equilibrium states are taken as initial states for many calculations, much computer time is spent in the spin-up phase of ocean model computations. In this note, we propose a new approach which can lead to a very large reduction in spin-up time for quite a broad class of existing ocean models. Our approach is based on so-called Jacobian–Free Newton–Krylov methods which combine Newton’s method for solving non-linear systems with Krylov subspace methods for solving large systems of linear equations. As there is no need to construct the Jacobian matrices explicitly the method can in principle be applied to existing explicit time-stepping codes. To illustrate the method we apply it to a 3D planetary geostrophic ocean model with prognostic equations only for temperature and salinity. We compare the new method to the ‘ordinary’ spin-up run for several model resolutions and find a considerable reduction of spin-up time.     

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.     

A note on using the accelerated convergence method in climate models

Application of the accelerated convergence (or asynchronous integration) method of Bryan (1984) to climate problems with time‐dependent forcing is investigated using an ocean general circulation model (GCM), with an idealized box ocean and forced with an idealized seasonal restoring surface temperature. Numerical experiments consist of a control experiment, which is integrated synchronously for 9000 years, and 2 experiments with asynchronous integration, one with depth dependent acceleration and one with depth independent acceleration. The latter 2 cases were integrated synchronously for 9000 years after asynchronous equilibria are reached. It is found that a few thousand years of synchronous integration is needed to reach a new equilibrium after asynchronous equilibrium is obtained, consistent with a scaling argument. However, at the new equilibrium, temperature in the deep ocean only differs from that of the early stage of synchronous adjustment by about 0.01°C. So for practical purposes, 50 years of synchronous integration beyond asynchronous equilibrium is sufficient. A simple interpretation of the accelerated convergence method of Bryan is also presented.     

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.     

An ocean model for climate variability studies

Models of the time dependent ocean circulation can be simplified considerably by filtering out all short term, small scale motions which are unimportant for climatic processes. For time scales large compared with a day and space scales large compared with the internal Rossby radius of deformation (～50 km), the currents in most of the interior ocean can be determined diagnostically as quasi-equilibrium fields, so that only the salinity and temperature fields need be treated prognostically.Regions of closed f/h contours, however, represent exceptions. Here trapped vorticity gyres exist as free flow solutions without external forcing, and in the presence of forcing the barotropic velocity field must therefore be determined prognostically through a potential vorticity equation for the gyres.Lateral boundary layers and the equatorial regions also require separate treatment. These were not considered specifically, but it is suggested that integrated (parametrical) models analogous in structure to mixed-layer models or the integrated boundary layer models of aerodynamics may be the most appropriate technique for coupling these regions to the interior ocean in a comprehensive ocean model suitable for climate studies.A coupled multi-region model of the global ocean circulation based on these scale considerations could be sufficiently cost-effective to permit systematic investigation of the role of the oceanic heat storage and transport in climate variability studies over a wide spectrum of space and time scales.The analysis of the seasonal variations of the interior ocean circulation represents a simple example in which the filtered model yields considerably simpler and more readily interpretable results than a fully three-dimensional, unfiltered model.     

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.     

两种热通量边界条件对热带太平洋海温模拟的影响

利用中国科学院大气物理研究所大气科学和地球流体力学数值模拟国家重点实验室发展的气候海洋模式(LICOM),考察了两种热通量边界条件(牛顿冷却型边界条件和总体公式型边界条件)对热带太平洋海温平均态和年际变率模拟效果的影响.结果显示,在两种边界条件下,模式均能较好的再现海温的年平均空间分布特征和季节循环特征.对比分析发现,...     

A comparison of two formulations of barotropic–baroclinic splitting for layered models of ocean circulation

In numerical models of ocean circulation, it is widespread practice to split the fast and slow motions into barotropic and baroclinic subsystems, respectively. In the case of the baroclinic equations, the dependent variables can either be (1) slowly-varying baroclinic quantities, obtained from splitting the original flow variables into barotropic and baroclinic components, or (2) the original unsplit variables, which can vary on both the fast and slow time scales. In the second case, the variables in each layer are adjusted after each (long) baroclinic time step to ensure compatibility with the results produced from the barotropic equations. The second approach can be applied to the layer thickness equation to ensure exact conservation of mass within each layer. In the case of the momentum equations, the second approach amounts to replacing unresolved fast portions of Coriolis and pressure forcing with time averages of well-resolved forcing from the barotropic system. In this study, both approaches for the momentum equations are evaluated, in several test problems, by comparing to analytical solutions or to solutions computed with an unsplit code that uses short time steps. The two methods give very similar results in some simple problems for which analytical solutions are known. However, in some eddying double-gyre simulations, the formulation with unsplit variables requires a significant reduction in the baroclinic time step in order to avoid numerical difficulties that include grid noise and inaccurate representation of the flow field. In contrast, the formulation with split variables does not display such difficulties, and in those same examples it can be used with zero explicit horizontal viscosity. All of these computations employ a two-level time-stepping method that was previously developed by the author.     

Efficient identification of ocean thermodynamics in a physical/biogeochemical ocean model with an iterative Importance Sampling method

Efficient identification of parameters in numerical models remains a computationally demanding problem. Here we present an iterative Importance Sampling approach and demonstrate its application to estimating parameters that control the heat uptake efficiency of a physical/biogeochemical ocean model coupled to a simple atmosphere. The algorithm has similarities to a previously-developed ensemble Kalman filtering (EnKF) method applied to similar problems, but is more flexible and powerful in the case of nonlinear models and non-Gaussian uncertainties. The method is somewhat more computationally demanding than the EnKF but may be preferred in cases where the approximations that the EnKF relies upon are unsound. Our results suggest that the three-dimensional structure of ocean tracer fields may act as a useful constraint on ocean mixing and consequently the heat uptake of the climate system under anthropogenic forcing.     

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.     

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.     

Numerical modeling of nonlinear resonance of semi-enclosed water bodies: Description and experimental validation

The purpose of this work is to a present a numerical model to solve a set of modified Boussinesq equations to analyse nonlinear resonance of semi-enclosed water bodies. The equations are solved on a finite element unstructured grid in order to achieve an optimal mesh resolution with the local geometry. The model is able to simulate long time lapses and realistic forcing in real bathymetries with a reasonable computational cost. To validate the numerical results, a set of experiments was carried out in a physical model of two adjacent elongated basins. Comparisons between numerical and experimental results for different geometries and nonlinear conditions show that the model is able to simulate with an excellent agreement the transient nonlinear resonant process.     

Linear and nonlinear wave diffraction using the nonlinear time dependent mild slope equations

Nonlinear wave diffraction is studied using the nonlinear time-dependent mild slope equation. The equations are solved using a combined Newton–Raphson and Crank–Nicolson finite difference scheme. The model results are verified for propagation of highly nonlinear waves over uniform depth and wave diffraction due to semi-finite breakwater and breakwater gap with different widths. Comparison between the numerical and experimental results indicates that the model is capable of simulating nonlinear wave diffraction. The model is applied to study the effect of the wave nonlinearity on the diffraction coefficient for a semi-infinite breakwater and a breakwater gap.     

Modeling the evolution of the world ocean ice cover in the 20th and 21st centuries

The current state of the simulation of sea ice cover as a component of new-generation global climate models is considered. Results from the model ensemble simulation of the observed world ocean ice cover, including its evolution in the 20th century, are analyzed, and projection of possible changes in the 21st century for three scenarios of anthropogenic forcing of the climate system are described. Unresolved problems and priorities for sea ice modeling are discussed.     

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.     

Weak and strong constraint data assimilation in the inverse Regional Ocean Modeling System (ROMS): Development and application for a baroclinic coastal upwelling system

We describe the development and preliminary application of the inverse Regional Ocean Modeling System (ROMS), a four dimensional variational (4DVAR) data assimilation system for high-resolution basin-wide and coastal oceanic flows. Inverse ROMS makes use of the recently developed perturbation tangent linear (TL), representer tangent linear (RP) and adjoint (AD) models to implement an indirect representer-based generalized inverse modeling system. This modeling framework is modular. The TL, RP and AD models are used as stand-alone sub-models within the Inverse Ocean Modeling (IOM) system described in [Chua, B.S., Bennett, A.F., 2001. An inverse ocean modeling system. Ocean Modell. 3, 137–165.]. The system allows the assimilation of a wide range of observation types and uses an iterative algorithm to solve nonlinear assimilation problems. The assimilation is performed either under the perfect model assumption (strong constraint) or by also allowing for errors in the model dynamics (weak constraints). For the weak constraint case the TL and RP models are modified to include additional forcing terms on the right hand side of the model equations. These terms are needed to account for errors in the model dynamics.Inverse ROMS is tested in a realistic 3D baroclinic upwelling system with complex bottom topography, characterized by strong mesoscale eddy variability. We assimilate synthetic data for upper ocean (0–450 m) temperatures and currents over a period of 10 days using both a high resolution and a spatially and temporally aliased sampling array. During the assimilation period the flow field undergoes substantial changes from the initial state. This allows the inverse solution to extract the dynamically active information from the synthetic observations and improve the trajectory of the model state beyond the assimilation window. Both the strong and weak constraint assimilation experiments show forecast skill greater than persistence and climatology during the 10–20 days after the last observation is assimilated.Further investigation in the functional form of the model error covariance and in the use of the representer tangent linear model may lead to improvement in the forecast skill.     

The effect of ocean tides on a climate model simulation

We implemented an explicit forcing of the complete lunisolar tides into an ocean model which is part of a coupled atmosphere–hydrology–ocean–sea ice model. An ensemble of experiments with this climate model shows that the model is significantly affected by the induced tidal mixing and nonlinear interactions of tides with low frequency motion. The largest changes occur in the North Atlantic where the ocean current system gets changed on large scales. In particular, the pathway of the North Atlantic Current is modified resulting in improved sea surface temperature fields compared to the non-tidal run. These modifications are accompanied by a more realistic simulation of the convection in the Labrador Sea. The modification of sea surface temperature in the North Atlantic region leads to heat flux changes of up to 50 W/m². The climate simulations indicate that an improvement of the North Atlantic Current has implications for the simulation of the Western European Climate, with amplified temperature trends between 1950 and 2000, which are closer to the observed trends.     

A diagnostic stabilized finite-element ocean circulation model

A stabilized finite-element (FE) algorithm for the solution of oceanic large scale circulation equations and optimization of the solutions is presented. Pseudo-residual-free bubble function (RFBF) stabilization technique is utilized to enforce robustness of the numerics and override limitations imposed by the Babuška–Brezzi condition on the choice of functional spaces. The numerical scheme is formulated on an unstructured tetrahedral 3d grid in velocity–pressure variables defined as piecewise linear continuous functions. The model is equipped with a standard variational data assimilation scheme, capable to perform optimization of the solutions with respect to open lateral boundary conditions and external forcing imposed at the ocean surface. We demonstrate the model performance in applications to idealized and realistic basin-scale flows. Using the adjoint method, the code is tested against a synthetic climatological data set for the South Atlantic ocean which includes hydrology, fluxes at the ocean surface and satellite altimetry. The optimized solution proves to be consistent with all these data sets, fitting them within the error bars.The presented diagnostic tool retains the advantages of existing FE ocean circulation models and in addition (1) improves resolution of the bottom boundary layer due to employment of the 3d tetrahedral elements; (2) enforces numerical robustness through utilization of the RFBF stabilization, and (3) provides an opportunity to optimize the solutions by means of 3d variational data assimilation. Numerical efficiency of the code makes this a desirable tool for dynamically constrained analyses of large datasets.     

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.