The GFDL Global Ocean and Sea Ice Model OM4.0: Model Description and Simulation Features

We document the configuration and emergent simulation features from the Geophysical Fluid Dynamics Laboratory (GFDL) OM4.0 ocean/sea ice model. OM4 serves as the ocean/sea ice component for the GFDL climate and Earth system models. It is also used for climate science research and is contributing to the Coupled Model Intercomparison Project version 6 Ocean Model Intercomparison Project. The ocean component of OM4 uses version 6 of the Modular Ocean Model and the sea ice component uses version 2 of the Sea Ice Simulator, which have identical horizontal grid layouts (Arakawa C‐grid). We follow the Coordinated Ocean‐sea ice Reference Experiments protocol to assess simulation quality across a broad suite of climate‐relevant features. We present results from two versions differing by horizontal grid spacing and physical parameterizations: OM4p5 has nominal 0.5° spacing and includes mesoscale eddy parameterizations and OM4p25 has nominal 0.25° spacing with no mesoscale eddy parameterization. Modular Ocean Model version 6 makes use of a vertical Lagrangian‐remap algorithm that enables general vertical coordinates. We show that use of a hybrid depth‐isopycnal coordinate reduces the middepth ocean warming drift commonly found in pure z* vertical coordinate ocean models. To test the need for the mesoscale eddy parameterization used in OM4p5, we examine the results from a simulation that removes the eddy parameterization. The water mass structure and model drift are physically degraded relative to OM4p5, thus supporting the key role for a mesoscale closure at this resolution.


Framing the OM4 Development
Numerical models are an essential part of the arsenal of science tools used to mechanistically understand how the ocean works within the earth climate system. Our goal for this paper is to scientifically document the foundations for the new generation of global ocean/sea ice models developed at the Geophysical Fluid Dynamics Laboratory (GFDL) of the National Oceanic and Atmospheric Administration (NOAA). These new models are collectively labeled the fourth generation of ocean-sea ice models (OM4), with OM4 providing the ocean and sea ice components for GFDL's global coupled and Earth system models, CM4 and ESM4, respectively. These models are used in Version 6 of the Coupled Model Intercomparison Project (CMIP6; Eyring et al., 2016), including the Ocean Model Intercomparison Project (OMIP) (Griffies et al., 2016). This paper serves as a foundation for future studies using OM4 in basic science research as well as for applications such as seasonal forecasting and climate projections.
GFDL has a long history of developing global ocean-sea ice models and corresponding coupled climate models to study the Earth system on seasonal-to-centennial time scales (e.g., Delworth et al., 2002Delworth et al., , 2006Delworth et al., , 2012Dunne et al., 2012;Gnanadesikan et al., 2006;Griffies et al., 2005Griffies et al., , 2011Manabe et al., 1991;Rosati et al., 1997). Model development is a process rather than an event, with ongoing development motivated by the need to improve simulation fidelity and capability. Three strategies for improving simulations from ocean-sea ice climate models include (i) enhanced formulations of the dynamical core, such as treatment of the vertical coordinate, (ii) improved parameterization of unresolved physical processes, such as mesoscale eddies and boundary layers, and (iii) refined numerical resolution to better capture both resolved and unresolved physical processes. We describe in this paper the result of employing all three strategies, making use here of a new ocean model dynamical core, improvements to physical parameterizations, and refinements in the horizontal and vertical grid spacing. In the remainder of this introductory section we present our strategies used to design and construct OM4.

Ocean Model Formulation and Spurious Diapycnal Mixing
The formulation of ocean circulation models, in particular the method used to discretize the vertical direction, has been identified as a source for spurious diapycnal mixing that can manifest as spurious heat uptake and modified water masses Legg et al., 2006;Winton et al., 1998). We hypothesize that this issue remains a primary source for model drift found in the traditional Eulerian geopotential coordinate models, particularly those that admit large tracer variance at the grid scale. The accumulation of tracer variance at the grid results from either maximizing the grid Reynolds number or by using a resolution that incompletely resolves those transient mesoscale eddies that are admitted Ilicak et al., 2012;Lee et al., 2002). Recent studies have shown that Lagrangian vertical coordinates can alleviate the problem (Leclair & Madec, 2011;Petersen et al., 2015) by better constraining the transport across vertical coordinate surfaces. Indeed, Lagrangian methods as employed by isopycnal models can eliminate spurious mixing for idealized test cases (Ilicak et al., 2012), thus providing a suitable numerical tool for studying adiabatic dynamics. However, isopycnal coordinates have poor vertical resolution in weakly stratified regions commonly present in high latitudes, thus making isopycnal coordinates problematic for realistic global climate modeling.
As anticipated by , a hybrid coordinate approach has emerged as a promising method to reduce spurious mixing in realistic ocean climate models. In the approach favored for climate (Bleck, 2002), and as pursued here, the chosen hybrid coordinate combines isopycnal coordinates in the ocean interior with a mixed layer resolving z * coordinate in unstratified regions. (The vertical coordinate z * is very close to geopotential but allows for the incorporation of sizable free surface undulations. This vertical coordinate was used by the MOM5 code in the ESM2M Earth system model of Dunne et al. (2012). See Stacey et al. (1995) and Adcroft and Campin (2004) for details of z * .) This approach has been demonstrated to improve the vertical structure of tracer fields and to in turn reduce spurious heat uptake (Chassignet et al., 2003;Megann, 2018;Megann et al., 2010). We make use of such a hybrid coordinate as implemented using a variant of the Arbitrary Lagrangian-Eulerian (ALE) method, with details given in section 2.1 and Appendix A. We demonstrate a significant reduction in ocean heat uptake and reduced model drift when choosing a hybrid z * -isopycnal relative to the use of z * everywhere.

Horizontal and Vertical Grid Spacing
The computational cost of an ocean model is primarily a function of the grid spacing, physical parameterizations, and the extent of ocean biogeochemical and/or ecosystem components. We here discuss considerations for both horizontal and vertical grid spacing chosen for OM4p5 and OM4p25. 1.2.1. Horizontal Grid Finite computational resources have meant that GFDL has used nominally 1°zonal resolution within its trunk climate models since CM2.1 (Delworth et al., 2006), ESM2M and ESM2G (Dunne et al., 2012), and CM3 . This grid spacing does not admit transient ocean mesoscale activity. Hence, the simulations have reduced ocean variability at the corresponding weekly to monthly time scales of the eddies. Their simulation integrity is furthermore quite sensitive to details of the mesoscale parameterizations. Delworth et al. (2012) provided evidence for improvements in certain climate features upon better resolving horizontal heat transport by transient mesoscale eddies in the 1/4°grid GFDL/CM2.5. Further benefits for ocean heat uptake and model drift were realized by even further resolving mesoscale eddies with 1/10°grid spacing in the GFDL/CM2.6 climate model (Delworth et al., 2012;Griffies et al., 2015). Marzocchi et al. (2015) also showed continued improvement with refined resolution in a North Atlantic regional model at comparable resolutions. Resolving mesoscale eddies yields a simulation that is more representative of the real ocean in many facets and further enhances the ability to represent fine scale air-sea interactions to determine a robust vertical grid algorithm suitable for the World Ocean, including the distinct temperature/salinity stratifications found in the North Atlantic and Southern Ocean.

Subgrid Scale Parameterizations
The parameterizations of unresolved processes are critical for the mean state and transient behavior of ocean climate models. We made many innovations in developing OM4, including improvements to parameterizations of the surface boundary layer, the bottom boundary layer, submesoscale mixed layer instabilities, and transient mesoscale eddies. We set out to use the same parameterizations across all resolutions, relying on scale awareness within the parameterizations to account for changes in grid resolution within a single model and between models. This ambition was not fully realized. We thus optimized a selection of parameterizations to obtain viable solutions across the two distinct resolutions of OM4p25 and OM4p5. Nonetheless, we contend that the two OM4 configurations are more similar than different. Indeed, we consider the similarity in simulations to be encouraging for the goal of using scale aware parameterizations that can be employed independently of resolution.

MOM6 Code Development and OM4 Parameterizations
The development of OM4 occured in parallel to the development of the Modular Ocean Model version 6 (MOM6), with the needs of OM4 setting priorities for MOM6 development. Some of the MOM6 algorithm was ported from a prototype code known as GOLD (Dunne et al., 2012). Even so, nearly all physical parameterizations and analysis code were newly developed to work in MOM6's new general coordinate framework. This approach prompted us to first target the OM4 configuration relying on the fewest essential parameterizations, namely OM4p25. Although more costly than OM4p5, the initial focus on OM4p25 allowed us to conduct relevant model development before developing code for the mesoscale eddy parameterization. The early lack of code for other parameterizations also necessitated prioritizing the order of parameterization implementation and development. We chose to order the implementation based on complexity and time scales of the impacted phenomena.
Ocean parameterizations can be broadly grouped into lateral and vertical processes. Lateral parameterizations are relatively complicated to implement in general vertical coordinates because of the need to calculate neutral directions and the orientation of tracer fluxes. In contrast, the parameterization of vertical processes, that are mostly expressed in terms of mixing coefficients, are more readily ported to generalized vertical coordinates. The reason is that a general vertical coordinate grid merely appears like a variable resolution grid with different resolutions in each vertical column. Hence, many of the vertical mixing parameterizations were fairly rapidly ported from other codes or obtained through other projects. Notably, some of these vertical schemes impact on long-time scale phenomena, such as mixing by breaking of internal waves (Melet et al., 2012;section 2.2.3), and others on shorter time scales, such as surface boundary layer mixing (Large et al., 1994;Reichl & Hallberg, 2018; section 2.2.1). Even so, not all parameterizations used in earlier GFDL models are used in OM4 either because we did not develop code in time to meet our deadlines for CMIP6, or because they were obsolete.
A list of physical process parameterizations sorted by the time scales of their role in the ocean is loosely determined by the distance from the surface that the process acts. Surface-driven processes, such as mechanically driven mixing and buoyancy-driven convection, act rapidly to modify water masses. Restratification of the mixed layer by mixed layer instabilities (submesoscale baroclinic eddies; section 2.2.2) is a rapid process that dictates water mass properties. Thereafter, mesoscale eddies and horizontal circulation modify and govern redistribution on the large scale (section 2.2.6). Thus, when building OM4, we focused on optimization of the mixed layer representation using relatively short simulations of a few years and gradually added longer time scale processes evaluated in systematically longer simulations of decades to centuries. 1.3.2. Aspiring to Remove Dimensional Parameters Many parameterizations rely on uncertain or unconstrained coefficients. Interactions between the parameterizations can lead to nonlinear behaviors in parameter sensitivity experiments. Where possible, we avoided the use of dimensional coefficients, preferring nondimensional parameters that can more readily be applied universally, and often determined in large eddy simulations (LES) and laboratory studies. For example, coefficients in the surface boundary layer parameterization (section 2.2.1) are all nondimensional and determined by fitting the single-column parameterization to LES and highly accurate (yet expensive) higher moment closures (Reichl & Hallberg, 2018). That is, we did not optimize the surface boundary layer parameterization in the context of the OM4 model configurations. This approach greatly reduced the parameter space dimensionality when building OM4. It also forced us to seek explanations outside of the boundary layer scheme when investigating causes for biases.
However, there are notable exceptions to our goal of using nondimensional parameters found in other OM4 parameterizations. One concerns the use of a dimensional background diffusivity to represent unresolved nonlocal breaking internal gravity waves (section 2.2.3). Another is the parameterization of mixed layer baroclinic eddies (Fox-Kemper et al., 2011;section 2.2.2), with this scheme containing a dimensional frontal length for which there is no theory. Both the background diffusivity and the frontal length have sizable impacts on the integrity of OM4, thus requiring specific optimization exercises within the OM4 configurations.

Tools and Protocols for Code and Model Development
The exercise of developing a product such as a coupled model or ocean-sea ice model is complicated and requires significant resources in the form of computing hardware/software and human brainware. To help organize our development and enhance reproducibility, we employed tools and protocols that could be valuable for other groups. Notably, we applied software version control using the git tool (https://git-scm.com/) to the entirety of the model configuration, including the processing of input data, management of the computational code, input parameters, platform environment and analysis scripts. We also employed an open development paradigm in which all versioning activity was publicly conducted on GitHub at their website (https://github.com/NOAA-GFDL/MOM6-examples). This approach provided a detailed record of decisions and exposed the code and model configurations to the scrutiny of a broad community of users and developers. Another major benefit is that the solutions we describe here should be bitwise reproducible by anyone with access to the same hardware. For those without the same hardware, the results should be broadly reproducible scientifically on any system capable of computations of this scale. In our experience, the above process for code and model development greatly enhances the reproducibility of numerical models used for ocean and climate science.

Documenting the Development Pathway
The purpose of this paper is to document the developments we made in constructing OM4 as part of the CM4 coupled simulations. As CM4 formed our primary goal for model development, development focused almost exclusively on coupled simulations, with the ocean/sea ice OM4 configuration run intermittently. To quantify this focus in the development of CM4, model development groups at GFDL collectively ran a total of 45,000 model years of coupled simulations and only 3,000 model years of the forced ice-ocean configurations.
The key reason we focused on CM4 rather than OM4 is that forced ocean-sea ice models are missing key coupled feedbacks that can hide or distort model improvements (e.g., Griffies et al., 2009). Furthermore, improvements in a forced ocean-sea ice model may not necessarily translate to improvements in the climate model with an interactive atmosphere.
The results detailed in this paper use the configuration of the ocean-sea ice components also used by the climate model CM4. We are rather pleased that this configuration yields respectable simulations under an ocean-sea ice experimental protocol. This success could be ascribed to good fortune. However, we suggest that it is also due to a compatibility between the state of the coupled climate model and the atmospheric data used in the forced ocean-sea ice models. That is, the success of the forced ocean-sea ice models documented in this paper is a testament to the closely coordinated efforts of the ocean, sea ice, atmosphere, land, and coupled model development at GFDL. We hypothesize that the development of a climate model is optimized only with close coordination across component model development.
The remainder of this paper consists of the following sections. In section 2 we detail the model configuration, including the grid and the physical parameterizations. Further details of the model code and its algorithms are provided in Appendix A. In section 3 we offer an assessment of a suite of model diagnostics from OM4 when forced using the interannual Coordinated Ocean-sea ice Reference Experiment (CORE) protocol from Griffies et al. (2009) and Danabasoglu et al. (2014). The numerous published papers making use of this protocol provided us with important benchmarks during the development of OM4. We note that any assessment will be incomplete and subject to the variety of interests of those actively participating in the analysis. Nonetheless, our aim is to expose both the strengths and the weaknesses of the OM4 configuration, thus providing guidance for further developments. We close the main portion of the paper in section 4, offering summary remarks and points where our ongoing development is targeted.

Model Configuration
In this section we present the variety of details comprising the OM4 ocean and sea ice model configurations constructed using the MOM6 ocean code coupled to the Sea Ice Simulator version 2 (SIS2) sea ice code. The codes make use of a common infrastructure, the Flexible Modeling System (Balaji, 2012), which provides parallelism and couples the component models together. In addition to the GitHub repository (https:// github.com/NOAA-GFDL/MOM6-examples), the configurations are summarized for ES-DOC at this site (https://furtherinfo.es-doc.org/CMIP6.NOAA-GFDL.GFDL-CM4.omip1.none.r1i1p1f1).

Ocean Model Dynamical Core and Grid Resolution
In this subsection we describe the dynamical core and the grid resolution chosen for the MOM6 ocean component of OM4. Appendix A offers further details, including the model equations and data sets.

Vertical Lagrangian Dynamical Core
The dynamical core of MOM6 is based on the hydrostatic primitive equations formulated in their generalized vertical coordinate form (see section A.1). The treatment of the vertical direction uses a special limit of the ALE method. Our approach follows the pioneering efforts from Bleck (2002) for ocean models (see also Figure 1. Horizontal grids used in (a) the nominally 1/4°OM4p25 and (b) the nominally 1/2°OM4p5. Only every 20th/10th grid line is plotted, respectively. Black/blue/cyan lines indicate the Mercator grid, bipolar cap, spherical grid, respectively. In OM4p25, the red lines indicate a twisted south pole grid. In OM4p5, the magenta lines are the region of a cubic blend between a fine spherical grid (yellow) and the Mercator grid. Adcroft & Hallberg, 2006 for an overview), whereby methods used for isopycnal coordinate models are combined with the three-dimensional ALE method of Hirt et al. (1997). In this approach, the ocean equations are integrated forward for the duration of the thermodynamic time step using a strictly Lagrangian algorithm in the vertical; that is, each layer retains its volume although tracer can move across layers through diffusion. The thermodynamic time step is used for updating thermodynamic variables and is an integer multiple of the baroclinic-dynamics time step used for subcycling the momentum and continuity equations. After time stepping the ocean state forward in the Lagrangian limit, a new vertical grid is generated and the ocean state is mapped from its current vertical grid to the new vertical grid. The vertical grids can be arbitrary but here we use a hybrid of z * (Adcroft & Campin, 2004;Stacey et al., 1995) near the surface and a modified potential density (referenced to 2,000 dbar) elsewhere, again as motivated by Bleck (2002). Details about our choice for hybrid vertical coordinates are given in section 2.1.4, and further details on the MOM6 model equations and algorithm are provided in section A.1.

Horizontal Resolution and Grids
We make use of two horizontal grid resolutions for OM4, with the ocean and sea ice using the same grids. The nominal 0.25°configuration in OM4p25 has 1,440×1,080 cells in the zonal/meridional directions, and the 0.5°configuration in OMp5 has 720×576 cells. Both grids were constructed by stitching together patches of simpler grids (Figure 1). The final composite grids are described using the Mosaics grid specification of Balaji et al. (2007), which requires generating the grid at twice the nominal resolution to unambiguously define all the finite volume metrics (i.e., grid distances and areas) for the Arakawa staggered grids. The parameters determining the grid generation are recorded in Table 1. A Mercator grid with unit cell aspect ratio spans the equator and extratropics with meridional extents determined by the zonal degrees of freedom (indicated by black grid lines in Figure 1). The ocean model uses the Arakawa C-grid (Arakawa & Lamb, 1981) and we choose to place the zonal velocity component exactly along the equator, and to center the tracer cells there as well. Meridians in the Mercator grid are equally spaced longitudinally. The positions of the Mercator parallels are determined solely by the zonal resolution and the number of meridional rows. The meridional extents of the Mercator grid are resultant properties and used as constraints on the other grid patches.
A bipolar cap following Murray (1996) covers the northern end of the sphere exactly matching the northern limit of the Mercator grid at roughly 65°N (indicated by blue grid lines in Figure 1. See Table 1 for precise latitudes). The coordinate singularities are placed at 120°W and 60°E so that converged coordinate lines fall over land, namely in Siberia and North America, rather than over the ocean. Note. All grids were generated at twice the given resolution and the metric factors (i.e., grid distances and areas) calculated and stored in a preprocessing step. The ocean and sea ice models read, and then aggregate, the metric factors for the different staggered locations on their grid.
In OM4p5, a spherical grid with uniform latitudinal resolution with 54 rows is fitted between the southernmost latitude of 78°S and the southern edge of the Mercator grid that sits to the north (cyan grid in Figure 1). In OM4p25 the spherical grid patch has 110 rows also ending at 78°S. Additionally, OM4p25 has a 70-row displaced pole grid in the south designed to facilitate future use with ice shelf cavities (red region in Figure 1). Since this model does not yet have ice shelf cavities, the first 40 rows are all masked as land cells and thus discarded from the grid. The last 30 rows of the remaining grid thus have variable latitudes.
OM4p5 also has a band of enhanced meridional resolution along the equator (yellow region in Figure 1). Between ±5°the meridional resolution is set to a uniform 0.2564°, and beyond that the resolution is interpolated with a cubic function to match the Mercator grid at ±30°(magenta region in Figure 1). This refined equatorial resolution enhances intrinsic equatorial ocean variability associated with tropical instability waves.

Time Steps
MOM6 uses subcycling to time integrate the governing equations. The time step for baroclinic dynamics is limited by the fastest internal gravity wave speed and the horizontal grid spacing. We use 15 min for OM4p25 and 30 min for OM4p5. The barotropic time step is automatically determined every 2 hr and was of order 19 and 38 s for OM4p25 and OM4p5, respectively. Tracer advection is implemented iteratively to allow longer time steps. Hence, both models can use a 2-hr time step for time stepping the thermodynamic and tracer equations, which is notable for being longer than the coupling time step. The coupling between sea ice and ocean is 1 hr, which is sufficient to stabilize the coupled dynamic modes between the sea ice and ocean. The sea ice model uses a rheology time step of 50 s and a dynamics time step of 20 min for both OM4p25 and OM4p5.

Vertical Grid and its Resolution
Both the OM4p5 and OM4p5 model configurations use 75 degrees of freedom in the vertical direction with the same hybrid isopycnal-z * vertical coordinate prescription inspired by Bleck (2002). A target potential density referenced to 2,000 dbar, ρ 2 , is assigned to each interface (between layers) and chosen so as to resolve the dominant subsurface water masses in the World Ocean. When generating a vertical grid at each time step, independently in each column the position of each target potential density is found by interpolation from a profile of density, ρ p which has 1% of the compressibility retained but centered on p r =2,000 dbar; ρ p ¼ ρ S; θ; p r þ 0:01ðp−p r Þ ð Þ . This mild compressibility introduces an artificial stratification in the otherwise potential density and mildly refines resolution in deep weakly stratified regions, which can appear statically unstable in a pure potential density variable referenced to a distant pressure. Each interface is also assigned a minimum depth, in z * -space, which is used if the interpolated position of the isopycnal interface is shallower than this depth. The end result is that interfaces approximately follow isopycnals in the interior, but near the surface they suddenly become parallel to the sea surface. The depth of transition between isopycnal to z * is broadly shallowest in the tropics and deepens toward high latitudes ( Figure 14 provides an example of the vertical grid for high-latitude sections, shown as part of our discussion of dense gravity current overflows). When a new vertical grid is generated at the end of each thermodynamic step, all state variables are conservatively remapped between old and new grids (more details on the vertical remapping algorithm are provided in section A.3). For legacy reasons, the newly generated grid has a minimum thickness of 1 mm applied to layers. However, layers can subsequently become much thinner during the dynamics steps.

Vertical Location of State Variables
The prognostic state variables, namely horizontal velocity components u and v, salinity, S, and potential temperature, θ, are layer averaged in the vertical and governed by layer-integrated equations-see section A.1 for the layer-integrated governing equations. Therefore, the variables should be interpreted as centered in layers. The layer thickness, h, is the other model field of equal prominence to the state variables. The layer thickness describes the spatially and temporally variable vertical grid. Vertical diasurface velocity within a particular vertical coordinate can be interpreted as located at the interfaces between layers. Pressure and density enter the layer-integrated governing equations through line integrals around a finite volume , and so identifying a single location for these variables is misleading. Nevertheless, the most appropriate label for vertical staggering of variables is the Lorenz grid (Lorenz, 1960).

10.1029/2019MS001726
Journal of Advances in Modeling Earth Systems

Ocean Model Subgrid-Scale Parameterizations
In this subsection we detail the various subgrid-scale parameterizations used for the ocean model component of OM4.

Planetary Boundary Layer
Early in the development of OM4 we used the K-profile parameterization (KPP) surface boundary layer parameterization (Large et al., 1994) as obtained through the CVmix project (https://github.com/CVMix), in conjunction with a z * vertical coordinate throughout the water column. We struggled to obtain satisfactory solutions, which in retrospect was due to implementation choices when employing KPP and to then missing features in the package (Van Roekel et al., 2018). However, in the interim, Reichl and Hallberg (2018) developed a new energetically constrained parameterization of the surface boundary layer (EPBL) specifically for climate models. The parameterization provides vertical diffusivity and viscosity as well as the depth of active mixing (i.e., the boundary layer thickness).
The EPBL parameterization was evaluated in single-column mode against LES. It was also matched to reproduce k−ϵ solutions using the General Ocean Turbulence Model (Rodi, 1987;Schumann & Gerz, 1995;Umlauf et al., 2005) over a wide range of forcing and external parameters including rotation. Although EPBL potentially has many parameters, these were highly constrained in the process studies of Reichl and Hallberg (2018) using LES and k−ϵ. We use EPBL for OM4 and furthermore retained essentially the same parameter values as recommended by Reichl and Hallberg (2018) including their suggested use of a Langmuir turbulence parameterization.

Mixed Layer Baroclinic Eddies
Restratification of the mixed layer by mixed layer baroclinic instabilities and mixed layer eddies (MLE) are relatively fast restratification processes (faster than interior baroclinic instabilities; Boccaletti et al., 2007). Furthermore, these upper ocean eddies have order 1 Rossby number and are thus commonly referred to as submesoscale eddies. In OM4, the effects from these eddies are parameterized by a bolus transport acting on scalars following the method of Fox-Kemper et al. (2011). In particular, we prescribe this transport following the main text of Fox-Kemper et al. (2011) yet without incorporating the extra caveats and filters described in their appendix.
The Fox-Kemper et al. (2011) MLE scheme requires a "mixed layer depth" over which the scheme is active. In early implementations within MOM6 we used the instantaneous boundary layer depth (as reported by KPP), and then the mixed layer depth (MLD) given by the depth of a chosen potential density difference with the surface. However, we found the scheme to be poorly behaved and noisier in regions with a weak pycnocline. We also noted that using the instantaneous MLD led to an overly strong diurnal cycle in the restratification (not shown).
We developed the following time-filtering method to handle the diurnal cycle. Here, the depth used as input for the MLE parameterization, H r , is chosen to be the deeper of the instantaneous mixing depth, H m , and a modified running mean of H r . The instantaneous depth, H m , is chosen to be the depth of active mixing as reported by the EPBL surface boundary layer scheme (section 2.2.1). The running time mean requires a relaxation time scale, τ r , which controls how slowly the depth of restratification shoals following a retreating mixed layer. When the mixed layer is deepening, we assume the instabilities are growing rapidly and thus begin to restratify immediately. In contrast, when the mixed layer is retreating (shoaling), the running mean keeps H r deeper, which we use to represent the finite lifetime of the mixed layer eddies. The time stepping for H r as seen by the MLE parameterization is thus given by This asymmetric filtering approach avoids an overly strong diurnal cycle in restratification and quite closely reproduces the overturning circulation described in Fox-Kemper et al. (2011;not shown). In all experiments τ r =30 days which allows the parameterized eddies to be seasonally modulated but ensures that restratification of the deepest mixed layer is persistent for several weeks. This time scale of 30 days is longer than the time scales of mixed-layer instabilities but was empirically found to help counter a cold bias in the coupled model for which this ocean model was developed.
The Fox-Kemper et al. (2011) parameterization requires the specification of a frontal length scale, which represents the dominant lateral length scale of the unstable mixed layer baroclinic eddies (L f in their equation (6)). The ratio of the model grid size to the frontal scale is used to scale resolved density gradients within the parameterization. It is also used to scale away the parameterized effects (i.e., turn the scheme off) when the frontal scale is resolved by the model grid. In the absence of a theory for the frontal length, we use a single value globally. Ideally, this length is the same at all grid resolutions. However, we found a strong sensitivity of SST biases to the MLE parameterization parameters, with the frontal length the most efficient parameter found to optimize for reducing biases. Hence, the OM4p25 and OM4p5 models use different frontal lengths (see Table 2). We test the sensitivity of mixed layer depths to this frontal length by considering OM4p5e, which uses the same frontal length as in OM4p25 (see section 3.2).

Interior Diapycnal Mixing
Diapycnal mixing due to shear instability is parameterized as a function of the resolved shear and the flux Richardson number following Jackson et al. (2008). The parameterization was designed for regions of dense gravity current overflows (section 3.7) as well as within the equatorial upper ocean (section 3.4). The authors provide a table of recommended parameter values based on comparisons to direct numerical simulations and LES, and we made use of their recommended values for OM4.
Mixing generated by breaking internal gravity waves generated by tidal forcing (i.e., internal tides), is parameterized using St Laurent et al. (2002) for the local dissipation at the internal tide generation sites. A map of average tidal flow amplitude is provided as external input to the model and is invariant over the course of the simulations. Tidal amplitudes, used in conjunction with a bottom roughness (see section A.4.1) in the internal tide dissipation scheme, were derived from a gridded model/data inverse estimate (Egbert & Erofeeva, 2002) using the eight primary tidal harmonic constituents summed in quadrature and linearly interpolated to the OM4 horizontal grid. The energy input to internal tides is calculated using the model bottom stratification and we tuned this energy in OM4p25 to be equivalent to 1 TW of dissipation of the global barotropic tides. We neglected to tune the energy input for OM4p5, using the same parameters for both models, so the dissipation in OM4p5 is 1.5 TW. One third of the dissipated tidal energy is used for mixing in the column above using Polzin (2009) to vertically distribute the mixing. All these aspects of the internal tide-driven mixing followed the implementation of Melet et al. (2012).
We parameterize the far-field breaking of propagating internal tides and internal waves generated by winds and other sources in an ad hoc manner by prescribing a static background tracer diffusivity. The chosen background has the latitude dependence proposed by Henyey et al. (1986) as implemented and tested by Harrison and Hallberg (2008). Here, we use a pivot value of 1.5×10 −5 m 2 /s at ±30°latitude and an equatorial value of 2×10 −6 m 2 /s.

Bottom Boundary Layer
Mixing in the bottom boundary layer is parameterized by the energetic model of Legg et al. (2006), in which 20% of the kinetic energy lost to bottom drag acting on the resolved flow is converted to potential energy via mixing. The stress due to the bottom drag is quadratic and takes the form ρ o c d u where u is the bottommost model velocity and c d =0.003 is a constant bottom drag coefficient. We include a constant velocity scale, u b , to represent the contribution of unresolved motions to the quadratic drag. In principle, u b should be spatially variable and would mostly reflect tidal amplitudes. But here we use a constant u b ¼ 0:1ms −1 .

Lateral Friction
Both OM4p25 and OM4p5 use a biharmonic viscosity taken as the maximum of a dynamic biharmonic viscosity with a coefficient of 0.06  and a static biharmonic viscosity of the form u 4 Δ 3 where Δ is local grid spacing and u 4 =0.01 m/s. In addition, OM4p5 has a Laplacian viscosity taken as the maximum of a Smagorinsky dynamic viscosity with coefficient 0.15 and a static viscosity of the form u 2 Δ with u 2 =0.01 m/s (Megann & New, 2001). The Laplacian viscosity is scaled to zero wherever the first baroclinic deformation is resolved, namely in the tropics, using the resolution function of Hallberg (2013).

Mesoscale Eddy Parameterizations in OM4p5
The coarser model, OM4p5, parameterizes mesoscale eddies through terms active in the tracer equation.
There are two distinct unresolved processes parameterized separately in OM4. The first is the relaxation of baroclinicity to remove available potential energy from the large scale, for which we follow the ideas of Gent et al. (1995) as implemented using the stream function formulation of Ferrari et al. (2010). Operationally, this scheme is implemented as an interface height diffusion to avoid problems with layer thickness diffusion described by Holloway (1997). Furthermore, we apply the associated eddy induced transport via a bolus velocity as commonly used in layer models (e.g., Bleck, 2002). The second parameterization provides a diffusive mixing of tracers along neutral directions, proposed by Solomon (1971) and Redi (1982), implemented using a finite-volume general-coordinate methodology.
The diffusivity applied to interface heights is prescribed using mixing length theory so that the nominal diffusivity for a column is given by is the eddy velocity scale, E MEKE is the kinetic energy per unit mass in the mesoscale eddies, and L e is the mesoscale eddy length scale. The eddy length and velocity scales are predicted using a two-dimensional model of unresolved mesoscale eddy kinetic energy (MEKE) following Jansen et al. (2015). We deviate slightly from Jansen et al. (2015) in that we (a) use a harmonic mean of the Rhines and "LH" scales for the mixing length rather than the minimum of the two, and (b) retuned the nondimensional parameters multiplying the length scales to produce realistic levels of MEKE. This scheme monitors the rate at which mean available potential energy is removed by the interface diffusion. Vertical structure of this MEKE-dependent diffusivity is imposed by multiplying the horizontally varying diffusivity by the normalized vertical structure of an equivalent barotropic mode. This modal function is calculated from the model stratification in each column by solving for the first baroclinic mode but with a no-slip bottom boundary condition (Hallberg, 1997). The resulting function is zero at the bottom, increases monotonically away from the bottom with no zero crossings, and is generally surface intensified.
The diffusivity used for the neutral diffusion of tracers has two contributions. One contribution is equal to the two-dimensional MEKE-dependent diffusivity as given above. The other is a constant prescribed background value of 50 m 2 /s added to ensure some lateral mixing of biogeochemical tracers on the equator and in the tropics. OM4p5 also includes a MEKE-dependent viscosity as well as a latitude-dependent background viscosity with a maximum value of 2,000 m 2 /s at the poles.

Sea Ice and Iceberg Models
In this subsection we describe the sea ice and iceberg models used in OM4 and CM4. Although the iceberg module is enabled and part of OM4, there is no prescribed calving in the forcing data set so the icebergs do not play a role in these simulations.

Sea Ice Model
For OM4 we make use of the SIS2.0 sea ice model, which is an updated version of the SIS model described by Delworth et al. (2006). SIS2.0 uses thermodynamics similar to the Community Sea Ice model CICE4.1 (Hunke & Lipscomb, 2010). SIS2.0 also uses a C-grid horizontal stencil as in MOM6 and an elastic-viscous-plastic rheology following Bouillon et al. (2013). Here we briefly describe the model formulations for column thermodynamics, ice thickness distribution, and dynamics.
SIS2.0 has four sea ice layers and one snow layer with thermodynamics similar to Bitz and Lipscomb (1999), except that brine has the heat capacity of seawater rather than that of sea ice. Constant conductivities of 0.31 and 2.03 W·m −2 ·K −1 are used for snow and ice, respectively. Salinities in the ice layers are prescribed as in CICE4.1 (Hunke & Lipscomb, 2010). The thermodynamic solver couples the ice temperatures implicitly to the atmosphere through a skin temperature as in SIS (Winton, 2000). Enthalpy (heat) and water are conserved to machine roundoff. SIS2.0 uses the sea ice shortwave radiative transfer method of Briegleb and Light (2007), employing the same code used in CICE4.1 (Hunke & Lipscomb, 2010). The optical property tunings of snow, ice, and pond are set one observational standard deviation high (R_snw = R_ice = R_pnd = 1) to adjust CM4.0's sea ice cover toward observations. SIS2.0 uses a constant coefficient for ice ocean temperature coupling of 240 W·m −2 ·K −1 .
Five sea ice thickness categories are managed according to a Lagrangian scheme (Bitz et al., 2001). SIS2.0 does not use a subgrid ice ridging scheme, so the thickness categories are concentrated on low thicknesses in order to resolve thermodynamic processes. The separating thickness boundaries for the five thickness levels are 0.1, 0.3, 0.7, and 1.1 m. The thinnest category extends to zero thickness and thickest category is unbounded. An up/down pass through the categories each time step moves ice between thickness categories when bounds are exceeded, which follows that in SIS (Delworth et al., 2006). SIS2.0 does not include a side melting scheme. Frazil is added as an areal extension of the thinnest layer of thickness category containing ice.
Ice and snow thicknesses are advected with a modified upwind scheme. Tracer advection uses a third-order piecewise parabolic method. The elastic-viscous-plastic rheology is solved on the C-grid following Bouillon et al. (2009). Ice strength is calculated based on cell-mean thickness and concentration as in Hibler (1979). The dependence on concentration is exponential in this formulation and can lead to instabilities. To avoid this problem, an estimated post time step concentration is used for calculating the ice strength. A quadratic drag formulation is used at the bottom of the ice with a dimensionless ice-ocean drag coefficients of 3.24×10 −3 . The surface roughness length is set to 0.1 mm.

Iceberg Model
A subcomponent of the sea ice model represents icebergs. In the CM4 and ESM4 coupled models, frozen water from the land model component is handled by an iceberg model which is part of the sea ice component as described in Martin and Adcroft (2010). We use updates to rolling and melting parameterizations described in Stern et al. (2016). We also introduced a stochastic parameterization of unresolved eddy and tidal motions that is active only in shallow water and near coasts. This scheme helps avoid a persistent build up of icebergs where climatological model winds tend to blow icebergs into embayments. In the experiments described here, there is no calving in the CORE forcing data set so this component does not affect the solutions. Table 2 provides a summary of the model configurations assessed in this paper. The two baseline models are the OM4p25 0.25°configuration, which has no parameterization of mesoscale eddies in the tracer equation, and the OM4p5 0.5°configuration, in which unresolved mesoscales are parameterized (section 2.2.6). The OM4p5 model also has a different setting for the Fox-Kemper et al. (2011) submesoscale mixed layer restratification parameterization (section 2.2.2). To illustrate the difference in mixed layer depth arising from changing the restratification front length, we show some results for the OM4p5e configuration that sets the parameter to that used in OM4p25 (section 3.2). Additionally, to illustrate the role of the mesoscale eddy parameterization in OM4p5 (and OM4p5e) we turn off the mesoscale parameterizations in OM4p5n ("no mesoscale eddy parameterization") configuration, with a variety of diagnostics studying the importance of this change.

CORE Forced Simulations
We make use of the interannual CORE protocol from Griffies et al. (2009) and Danabasoglu et al. (2014) to force OM4 over the years 1948-2007. The CORE protocol is specified to be used in Phase I of Coupled Model Intercomparison Project version 6 Ocean Model Intercomparison Project (CMIP6/OMIP; Griffies et al., 2016). This approach prescribes a common atmospheric state and runoff as compiled by Large and Yeager (2009), and uses bulk formulae to compute the associated boundary fluxes as a function of the prognostic SST, surface ocean velocity, and ice cover. Notably, the bulk formula approach does not yield identical forcing across the suite of OM4 models, given that SST, surface ocean velocity, and sea ice cover are a function of the changing ocean state. Even so, CORE simulations are far more constrained than in a coupled climate model that uses interactive atmosphere and land models. Following the CORE interannual forcing (IAF) protocol, we ran the simulations for five consecutive cycles over years . The model clock starts in 1708 so that the fifth and final cycle corresponds to 1948 to 2007. When showing climatologies from the simulations, we focus on the most recent 20-year period of 1988-2007. The CORE protocol does not specify the restoring rate for surface salinity. In the experiments reported here, surface salinity is restored everywhere to climatology with a piston velocity of 1/6 m/day, which is the value used in earlier MOMbased simulations (e.g., Table 3 of Griffies et al., 2009 and Table 2 of Danabasoglu et al., 2014).

Surface Climate
In comparison to a fully coupled model, the SST in an ice-ocean model driven by a prescribed atmosphere is strongly controlled by the effective restoring to data. Thus the SST errors are somewhat smaller than in a coupled model but nevertheless indicative of model biases. Figure 2 shows the SST bias of the 20-year average at the end of the fifth cycle (years 1988-2007) relative to World Ocean Atlas (WOA) climatology. There is a general warm bias in the tropics of order 0.5°C which is very similar for all four models shown suggesting the bias is not a function of resolution (compare OM4p25 and OM4p5), eddy-parameterization (compare OM4p5n with OM4p5) or vertical coordinate (compare OM4p25 with OM4p25-z*). This mild tropical warm bias is consistent with the shallow mixed layers discussed in section 3.2.
In contrast to the lower latitudes, the choice of vertical coordinate leads to significant SST differences in the Southern Ocean with relative surface cooling in the Pacific and Indian sectors of the Southern Ocean.
The role of resolution and the need to parameterize mesoscale eddies are manifest in changes in the Western Boundary jets and Southern Ocean. The bias in the eddy-permitting 1/4°OM4p25 and eddy-parameterized 1/2°OM4p5 are very similar, while the 1/2°OM4p5n without parameterized eddies has significant warm biases in the Southern Ocean, a warm bias in the Labrador Sea, and very different bias in the Kuroshio extension. These results are consistent with the general view of the role of eddies (e.g., Farneti et al., 2010;Griffies et al., 2015;Hallberg & Gnanadesikan, 2006). One difference between OM4p25 and OM4p5 is a stronger cold bias in the northwest corner of the North Atlantic in OM4p5. This "blue spot" is often associated with weak overflows from the Greenland-Icelandic-Norwegian Seas, or a weak deep western boundary current (Zhang et al., 2011), and the same could be true here with a slightly weaker Atlantic meridional overturning circulation as discussed in section 3.6.

Mixed Layer Depth
The ocean surface mixed layer is directly influenced by air-sea exchange of heat, momentum, and gases from above, vertical mixing and entrainment due to various ocean turbulence processes, and subject to horizontal transport and overturning processes on a variety of space and time scales. The MLD and its variability indicates the depth over which atmospheric and oceanic properties are exchanged on seasonal, interannual, and longer time scales. There is a growing suite of in situ ocean measurements, particularly from Argo floats, that enhance our confidence in the observed MLD both temporally and globally. We compare the MLD computed from the model with that computed from the WOA climatology, providing a valuable means to assess model skill.
When the mixed layer deepens, heat and gaseous content of the seawater that have been recently exchanged with the local atmosphere penetrate vertically toward the ocean interior. In this analysis we consider the seasonal cycling of the deep MLDs (forced by buoyancy loss and strong winds), as the water masses influenced by these events are critical for the formation of the deep and intermediate water masses that fill the ocean interior. The globally deepest MLDs occur in high-latitude winter regions subject to strong buoyancy loss and strong winds, where convectively driven MLDs can exceed 1 km and variability in their composition is ultimately reflected in the properties of interior ocean seawater. Accurately predicting these deep mixed layers is critical for predicting carbon dioxide (CO 2 ) uptake and the ocean's overturning circulation.
The MLDs shoal in summers, where abundant levels of solar radiation maintain vertical stratification that limits maximum MLDs to Oð10Þ m or even less in the tropics. The biases in these regions contribute less directly to long-term ocean properties. Even so, capturing these shallow MLDs is critical for ocean climate simulations due to intraseasonal ocean-atmosphere feedbacks that affect large-scale atmosphere-ocean coupling (e.g., El Niño-Southern Oscillation), as well as their importance in biogeochemical cycling of compounds such as nutrients, oxygen, and carbon dioxide.

Observation-Based Measures
We use the WOA monthly temperature and salinity climatology to derive the observational MLD benchmark for our assessment of OM4. For this analysis, the MLD for each month is computed for 1°boxes. The minimum and maximum of the monthly MLD are then calculated and used to create a climatology with global coverage. The constant density threshold of de Boyer Montégut et al. (2004) is employed here to compute MLD in an identical manner for both the observation and the model. The MLD is therefore formally defined here as the thickness of the layer over which the vertical density variation is less then Δρ=0.03 kg/m 3 from the surface layer (of order 1-m depth), where the density is computed relative to the surface pressure (i.e., potential density σ 0 ). The model three-dimensional monthly averaged density over the final 20 years of the fifth CORE cycle is used to compute the model MLD for the data shown below. This averaging period represents the fraction of the WOA record with the most observations (including the ARGO period).

Minimum of Monthly Mean Mixed Layer Depths (Summer)
The minimum monthly mean observed MLDs and the various OM4 configurations are shown in Figure 3, which corresponds to summer MLDs in the respective hemispheres. The observed minimum MLDs are shallower than 30 m globally except for a deep band in the Southern Ocean, which is a pattern generally represented in both OM4p25 and OM4p5. The primary differences between observations and OM4p25 are a shift in spatial pattern and a slight net deep bias in the Southern Ocean and a slight band of deep bias relative to the observation in the tropics. The OM4p5 subtropical bias contrasts to the OM4p25 subtropical bias, with a slight shallow bias occurring throughout the majority of the region due to the strengthened submesoscale restratification parameterization (see Table 2). The Southern Ocean MLDs in OM4p5 also show a spatial shift in the occurrence of the deepest MLDs, though the net bias is shifted relative to OM4p25 (again due to the strengthened submesoscale restratification).
The influence of changing the horizontal resolution on minimum MLDs is isolated by comparing OM4p25 to OM4p5n (Figure 3d). The finer resolution of OM4p25 leads to slightly deeper minimum MLDs in the Southern Ocean and the subtropical regions. By introducing the mesoscale eddy parameterization, the differences between OM4p5e and OM4p25 are reduced in the Southern Ocean (compare Figure 3d and Figure 3e). The mesoscale eddy parameterization is not active in the tropical region (following Hallberg, 2013), so it is unknown if the eddy parameterization would reduce the differences found in the tropics. As expected, the effect of the strengthened submesoscale restratification parameterization is to shoal the MLDs in OM4p5 everywhere relative to OM4p5e due to the increased restratification effect (Figure 3f).

Maximum of Monthly Mean Mixed Layer Depths (Winter)
The maximum monthly mean observed MLDs and the various OM4 configurations are shown in Figure 4, which corresponds to winter MLDs in the respective hemispheres. The Southern Hemisphere maximum MLDs in both OM4p25 and OM4p5 tend to yield a deep MLD bias in the Southern Ocean. In the subtropics there is a shift in bias in OM4p5 relative to OM4p25, similar to that discussed in the minimum MLDs. The patterns of MLD bias are complex in the Northern Hemisphere high latitudes, with both OM4p25 and OM4p5 yielding deep biases in the Labrador Sea contrasting with the shallow biases in the central north Atlantic Ocean.
The influence of decreasing the resolution (Figure 4d) reveals that OM4p5n has significantly deeper MLDs than OM4p25 at high latitudes. Deeper MLDs for OM4p5n follows from the absence of a mesoscale eddy parameterization to restratify deep MLD regions, whereas OM4p25 explicitly permits eddies (though not fully resolved) that act to restratify (Danabasoglu et al., 1994). This overly deep MLD bias in OM4p5n is significantly improved by introducing the mesoscale eddy parameterization, as revealed by comparing OM4p5n to OM4p5e (Figure 4e). The difference due to reducing the resolution (Figure 4d) is nearly opposite to the difference due to introducing the eddy parameterization (Figure 4e). The effect of the strengthened submesoscale restratification parameterization is examined by comparing OM4p5 to OM4p5e (Figure 4f), with the difference between observations and OM4p25, (c) the difference between observations and OM4p5, (d) the difference between OM4p25 and OM4p5n, indicating the effect of reduced resolution, (e) the difference between OM4p5n and OM4p5e, indicating the effect due to the addition of the mesoscale eddy parameterization, and (f) the difference between OM4p5 and OM4p5e, indicating the effect due to the changed submesoscale parameterization. Note that differences and biases are plotted with the same color scale which saturates for some difference plots. Red indicates a shallow bias, often associated with a warm bias, whereas blue is a deep bias. submesoscale scheme shoaling the MLD everywhere except the deepwater production region in the Weddell Sea and Labrador Sea.

Discussion
Overall OM4 is skillful at simulating the MLD, particularly for the yearly minimum values (generally the summer) when the biases are very small over much of the ocean. The skill at predicting wind and wave driven MLDs supports the approach of optimizing the boundary layer parameterization from one-dimensional (Reichl & Hallberg, 2018) and LES (section 2.2.1). However, there are still significant biases in yearly maximum (winter) MLDs that are critical to understand and rectify in subsequent model development efforts. Furthermore, the minimal deep bias in subtropical latitudes contrasts with the warm subtropical SST bias and suggests a difference in the stratification within the seasonal thermocline.

Ventilation Measured by Ideal Age
The analysis of MLD (section 3.2) provides a partial measure of ventilation with direct comparisons to observational measures. We also find it useful to consider the ideal age as introduced by Thiele and Sarmiento (1990). The maximum ideal age is 300 years, which is the length of the simulation. In Figure 5, we show the ratio of the difference between the maximum age and model age at the end of the simulation, compared to the maximum age. All simulations show young, recently ventilated waters near the surface and older waters in the interior. The OM4p25 and OM4p5 simulations are remarkably similar both spatially (r 2 =0.98, Figures 5a and 5b), and in total ventilation below 3,000 m (4.45×10 16 m 3 vs. 3.9×10 16 m 3 , Figure 5d). The OM4p5n simulation displays very strong mixing to depth especially in subpolar regions (Figure 5c), with total ventilation below 3,000 m of 1.3×10 17 m 3 . The addition of mesoscale and submesoscale eddy parameterizations in OM4p5 allows this model to reproduce much of the features of ocean ventilation present in the finer resolution OM4p25. We conclude, as have numerous other studies, that these eddy parameterizations are critical for the representation of circulation and ocean heat and CO 2 ventilation at coarse resolution. Figure 6 shows the annual mean upper ocean temperature difference between the OM4 simulations and Tropical Atmosphere Ocean/Triangle Trans-Ocean Buoy Network (TAO/TRITON) mooring observations (McPhaden et al., 1998), at four mooring longitudes in the equatorial Pacific. All three OM4 simulations show an equatorial Pacific warm SST bias, in contrast to the cold SST bias that emerges in the fully coupled CM4 and ESM4 simulations (not shown). The OM4 simulations also show a cold bias below 200 m, which together with the warm surface bias leads to excessive thermal stratification (overly large ∂T/∂z) of the equatorial Pacific upper ocean. The stratification bias is particularly strong just above the thermocline, at 170°W (between 80 and 120 m), 140°W (20-50 m), and 110°W (0-30 m), suggesting a shallow thermocline and too little downward mixing of surface heat at these locations. Excessive near-surface thermal stratification also predominates in the off-equatorial tropical Pacific, as well as in the tropical Atlantic and Indian Oceans (not shown). Figure 6 shows that OM4p25's equatorial temperature biases are slightly reduced compared to OM4p5, although their vertical structure generally looks quite similar. OM4p5n actually has the smallest biases, with the least intense thermocline along the equator. Thus OM4p5's eddy parameterization worsens the Figure 5. Ocean ventilation inferred from ideal age tracer, where we plot the normalized age f=(A max −A)/A max , with A max the maximum age of the simulation (300 years). The normalized age has values between zero and unity, with zero indicating recently ventilated (young) waters and unity indicating unventilated (old) waters. We show zonal mean latitude-depth plots for (a) OM4p25, (b) OM4p5, and (c) OM4p5n. Panel (d) shows area-integrated ventilation below 3,000 m for the three simulations: OM4p25 (blue), OM4p5 (green), and OM4p5n (red). OM4p5n is clearly an outlier, with a significant amount more ventilation than OM4p5 and OM4p25. temperature biases relative to observations, although the parameterization does make OM4p5's equatorial Pacific upper-ocean temperature structure closer to that in OM4p25. Figure 7 shows the annual mean upper ocean zonal velocity difference between the OM4 simulations and the TAO/TRITON mooring observations in the equatorial Pacific. Consistent with the shallower thermocline in the OM4 simulations (Figure 6), the simulations also show a shallower equatorial undercurrent (EUC) than observed. The simulated EUC is also too strong, with excessive time mean vertical shear (∂u/∂z) in the western and central Pacific in the top 100 m of the ocean. Eastward velocity biases predominate at nearly all depths at all four mooring sites, indicating a strengthening of the eastward EUC and a weakening of the westward South Equatorial Current in the simulations. The resolution refinement and eddy parameterization have little effect on the time mean equatorial Pacific velocity profiles, though OM4p25 shows a slight weakening of the velocity biases.

Potential Temperature Trends
The 6-hourly prescribed atmospheric state tightly constrains the surface ocean properties under the CORE IAF forcing. The 60-year cyclic forcing is evident in the time series of global mean SST (Figure 8a) but the spread among the ensemble of configurations (OM4p25, OM4p5, OM4p5e, and OM4p5n) only varies between 0.01 and 0.05°C. Initially, OM4p5n is the coldest ensemble member but becomes the warmest by the end of the fifth forcing cycle. Despite some subtle differences, OM4p25, OM4p5, and OM4p5e are also well-balanced in terms of their net downward surface heat fluxes, with averages of less than 0.1 W/m 2 over the 300 years of simulation (Figure 8b). In the OM4p5n configuration, decadal-average heat fluxes into the ocean are initially 1.5 W/m 2 and remain positive over the first two forcing cycles before becoming negative over cycles three through five. Figure 7. As in Figure 6, but for upper-ocean zonal velocity differences (cm/s).  Wittenberg et al. (2018), the TAO/TRITON data are processed using daily means of measurements from both fixed depth mechanical current meters and acoustic Doppler current profilers, with a simple average of the two where both are available. The merged daily means from the moorings are vertically interpolated and then averaged onto a monthly grid. These monthlies are then used to compute a 12-month climatology, whose annual mean is subtracted from the simulations; 1988-2007 corresponds to the period of overlap between the OM4 simulations and the available observations at the four equatorial mooring sites.
The OM4p25 simulation is remarkably stable with a volume mean ocean potential temperature drift of 0.0006°C/century (Figure 8c). In contrast, the OM4p25-z* simulation is a clear outlier by manifesting a monotonic warming in excess of 0.25°C/century and a surface heat uptake in excess of 1 W/m 2 . Recall that the sole difference between OM4p25-z* and OM4p25 is the vertical coordinate, with OM4p25-z* using z* throughout the vertical whereas OM4p25 uses a hybrid isopycnal-z* (see Table 2). The depth versus time evolution of potential temperature (Figure 9a) shows slight cooling above 500-m depth, a maximum of warming near 2,000-m depth, and cooling in the abyssal waters below 4,000 m. The volume mean drift is slightly larger in OM4p5 (0.008°C/century) and OM4p5e (0.02°C/century). The depth versus time evolution of potential temperature is similar in the OM4p5 and OM4p5e configurations (Figures 9b and 9c), with much of the warm drift concentrated between 500 and 1,000-m depth. Both configurations also show cooling in the subthermocline waters below 1,000 m and there is enhanced warming relative to the 0.25°configuration between 3,000 m and 4,000 m. The abyssal cooling is similar in magnitude between the 0.5°and 0.25°c onfigurations. While OM4p5 and OM4p5e look very similar throughout the water column, they differ most in the upper 500 m where OM4p5 tends to cool slightly whereas OM4p5e warms slightly. This difference highlights the role of the submesoscale parameterization in model drift in the 0.5°c onfigurations, with the parameter setting used in OM4p5 resulting in stronger restratification than the other model configurations (see Table 2).
Unlike the other three configurations previously discussed, OM4p5n has centennial-scale structure in the volume mean potential temperature drift. The ocean warms at a rate of 0.2°C/century over the first two forcing cycles, driven primarily by warming in the upper 1,000 m. Beginning in the third cycle, this upper-ocean warming ceases and is replaced by cooling the abyssal waters. Over CORE-II cycles three through five, the volume mean ocean temperature cools at a rate of −0.1°C/century. The switch from upper ocean warming to abyssal cooling is related to changes in high-latitude ventilation and deep water formation. The OM4p5n experiment ran for an additional five cycles and the volume mean temperature continued its slight cooling trend of −0.02°C/century (Figure 8c, inset).

Poleward Heat Transport and Atlantic Overturning Circulation
In Figure 10 we show the directly calculated northward global ocean heat transport from OM4p25, OM4p5, and OM4p5n simulations. The OM4p25 and OM4p5 simulations exhibit very similar northward global ocean heat transport at all latitudes. The simulated maximum northward global ocean heat transport (∼1.5 PW around 20°N) is weaker than the implied maximum northward global ocean heat transport (∼1.8 PW around 20°N) derived from air-sea surface heat fluxes using National Centers for Environmental Predictions (NCEP) reanalysis  Also shown are the implied mean northward global ocean heat transport derived from air-sea surface heat fluxes using NCEP reanalysis data (Trenberth & Caron, 2001) and the observation-based in situ estimates (Ganachaud & Wunsch, 2003, G&W). Note the relatively stronger poleward heat transport for OM4p5n in the Southern Hemisphere whereas it has a weaker poleward transport in the midlatitude Northern Hemisphere. CORE = Coordinated Ocean-sea ice Reference Experiments; NCEP = National Centers for Environmental Prediction.
data (Trenberth & Caron, 2001). The result is similar to that found previously in some CORE simulations using different coarse-resolution ocean models with eddy parameterization (Griffies et al., 2009).
Compared with OM4p25 and OM4p5 simulations, the simulation with no eddy parameterization (OM4p5n) exhibits much weaker northward ocean heat transport in the Northern Hemisphere south of 40°N, and much stronger southward ocean heat transport in the Southern Hemisphere especially over the Southern Ocean. The difference is mainly due to the absence of parameterized subgrid-scale eddies in OM4p5n, which leads to unrealistic open ocean deep convection in the Southern Ocean (see more discussions in section 3.8), much colder subsurface temperature near Antarctic, and stronger southward volume transport in the western boundary currents in the Indo-Pacific Ocean. The above results indicate the important role of eddy parameterization in OM4p5, especially over the Southern Ocean.
The northward Atlantic heat transport from OM4p25, OM4p5, and OM4p5n simulations ( Figure 11) are similar to each other and slightly stronger in the tropical Atlantic and the midlatitude North Atlantic in the OM4p5n simulation. The simulated maximum Atlantic heat transport (∼1 PW near 26°N) is comparable to many previous CORE simulations (Danabasoglu et al., 2014) and weaker than the direct estimate (∼1.2 PW for the period of 2004-2013) from the RAPID measurements (Johns et al., 2010) and the implied value derived from air-sea surface heat fluxes using NCEP reanalysis data (Trenberth & Caron, 2001). The underestimation of the northward Atlantic heat transport around 26°N has been documented for many previous CORE simulations (Danabasoglu et al., 2014) and is likely due to the shallow meridional overturning.
The slightly stronger northward Atlantic heat transport in OM4p5n is related to the stronger Atlantic meridional overturning circulation (AMOC) strength in this simulation with a maximum volume transport of ∼18 Sv at 26°N, compared with OM4p25 and OM4p5 simulations that have a maximum volume transport of ∼15 and ∼14 Sv at 26°N, respectively ( Figure 12). The directly measured maximum AMOC transport at 26°N from the RAPID data is ∼17 Sv for the period of 2004-2015 (Cunningham et al., 2007;Smeed et al., 2018). Without the parameterized mesoscale eddy restratification effect, deep convection at northern high latitudes is very strong and contributes to the relatively stronger AMOC in OM4p5n. In all three simulations (OM4p25, OM4p5, and OM4p5n), the simulated AMOC penetration depth is shallower (∼3,000 m) than observational estimates (∼4,300 m) at 26°N from the RAPID array ( Figure 13; McCarthy et al., 2015), as also found in previous CORE simulations using many other ocean models (Danabasoglu et al., 2014). As discussed in section 3.7, the simulated shallower AMOC reflects a deficiency in modeling the deep penetrating Nordic Sea overflow and contributes to the simulated weaker-than-observed northward Atlantic heat transport (Danabasoglu et al., 2010,Wang et al., 2015Zhang et al., 2011).

Atlantic Overflows
In both OM4p25 and OM4p5, dense water moves from the Nordic seas and into the North Atlantic through the Denmark Straits at (27°W, 66°N) and the Faroe Bank Channel at (8°W, 62°N; Figure 14, upper panels). However, the density signal is rapidly eroded and does not penetrate below 1,000 m, across the grid Figure 11. Directly calculated 20-year (1988-2007) mean northward Atlantic heat transport (PW) from OM4p25, OM4p5, and OM4p5n CORE-II simulations. Also shown are the implied mean northward Atlantic heat transport derived from air-sea surface heat fluxes using NCEP reanalysis data (Trenberth & Caron, 2001) and the observation-based in situ estimates (Ganachaud & Wunsch, 2003, G&W). NCEP = National Centers for Environmental Prediction.
resolutions. Since the dense signal at the sills does not connect with the dense water at the bottom of the North Atlantic, the dense bottom water is not supplied by the Nordic overflows in these simulations. This rapid dilution of the overflow signal contrasts with observations (Girton & Sanford, 2003;Mauritzen et al., 2005) and simulations with an isopyncal coordinate ocean model (Wang et al., 2015), which all show dense water from the overflows descending well below 1,000 m.
In both of the Nordic overflow regions, the vertical coordinate system transitions from z * coordinates in the dense water upstream of the sill to isopycnal coordinates downstream of the sill (Figure 14, lower panels). The AMOC in the models is integrated from the bottom, so the simulated AMOC strength relative to the surface is defined as the difference between the maximum value (near 1,000 m) and the surface value (on the order of 1 Sv due to freshwater transport over the North Atlantic). Zonal integrals are calculated along model grid lines that are not truly zonal in the Arctic, north of 65°N, due to the bipolar grid north of 65°N (see Figure 1). The model overturning is the residual mean, which includes effects from parameterized mesoscale (in OM4p5) and submesoscale eddies.

Journal of Advances in Modeling Earth Systems
Some of the dilution in the overflow waters may therefore be due to numerical mixing during downslope descent in the z * coordinate regime (Winton et al., 1998). However, explicitly parameterized mixing is also significant (Figure 15), particularly the bottom boundary layer component of mixing. Observations (e.g., Mauritzen et al., 2005) suggest diffusivities as high as κ=10 −2 −10 −1 m 2 /s in the Faroe Bank Channel overflow, of a similar order of magnitude to our parameterized diffusivities. However, careful optimization of the bottom boundary layer mixing was not part of our model development process, with future development motivated by these results.
In contrast with the Nordic Sea overflows, OM4p25 generates a Mediterranean outflow with a structure that matches observational measures (Baringer & Price, 1997), such as an intermediate water mass with a salinity maximum at 1,200-m depth (Figure 16, top panels). The bottom boundary mixing scheme, which provides a significant fraction of the diffusion in the Mediterranean outflow ( Figure 16, lower panels) was specifically implemented to prevent formation of a split Mediterranean outflow plume (Legg et al., 2006(Legg et al., , 2009. By inducing the correct amount of mixing in the Mediterranean outflow, where intense mixing leads to a middepth outflow, too much bottom mixing can result in the Nordic outflows ( Figure 15). Future research is needed to determine the relative roles of parameterized mixing and vertical coordinate choice in causing excessive mixing in the Nordic overflows. We suspect that this excessive dilution of the Nordic overflows contributes to the shallow overturning in OM4p25, a connection demonstrated in Wang et al. (2015).

Southern Ocean
The Southern Ocean is home to the Antarctic Circumpolar Current (ACC), whose depth-integrated mass transport is the world's largest, and the Southern Ocean meridional overturning circulation, which plays a key role in the earth's climate system through the ventilation/sequestration of heat and carbon (e.g., Frölicher et al., 2014;Morrison et al., 2014;Rintoul & Naveira Garabato, 2013;Rintoul et al., 2001;Sabine et al., 2004). The ACC transport and the meridional overturning circulations are tightly woven into the fabric of Southern Ocean dynamics, with mechanisms affecting one also affecting the other. Given the importance of the Southern Ocean for climate, the scientific utility of a climate model is influenced by the quality of its Southern Ocean simulation. Consequently, there have been many studies of the Southern Ocean in CMIP5-  Figure 14. Sheardriven diffusivity (Jackson et al., 2008;top left), internal-tide driven diffusivity (Melet et al., 2012;top right), bottom boundary layer mixing (Legg et al., 2006; bottom left), surface boundary layer (Reichl & Hallberg, 2018;bottom right). Note in particular the rather large diffusivities in the bottom boundary layer, thus contributing to strong mixing of dense bottom overflow waters as they move downslope.

10.1029/2019MS001726
Journal of Advances in Modeling Earth Systems class coupled climate models. For example, Downes and Hogg (2013) studied the Southern Ocean circulation and mesoscale eddy compensation (defined below), Sallée et al. (2013) focused on the mixed layer depths and associated ventilation, and Heuzé et al. (2013) characterized the features of Antarctic Bottom Water. Downes et al. (2015) and Farneti et al. (2015) complemented these CMIP5 climate model studies by considering a variety of Southern Ocean features, including sea ice, in CMIP5-era CORE-II simulations. These and other studies motivate elements of the following assessment of the Southern Ocean as simulated by the suite of OM4 models.
We here compare three OM4 configurations: the standard OM4p25 and OM4p5 configurations as well as OM4p5n. The modification made to the mixed layer restratification parameterization in going from OM4p25 to OM4p5 makes a negligible difference in this region, thus allowing us to not include OM4p5e in this analysis. Differences between OM4p5n and OM4p5 are predominantly due to the mesoscale closure used in OM4p5.

Mass Transport Through the Drake Passage
Over the ACC latitudes, the maximum in westerly wind stress increases more than 50% during the 60-year CORE-II cycle and it trends southward by a few degrees latitude (see Figures 1 and 2 from Farneti et al., 2015), in agreement with the observed positive trend in the Southern Annular Mode (Marshall, 2003;Thompson et al., 2000). We expect the depth-integrated zonal mass transport to increase if it is directly Figure 16. Sections along the Mediterranean overflow at 35.8°N, from 5°W to 9°W, from OM4p25. Top left shows potential density referenced to 2,000 dbar. Top right shows salinity. Lower four panels show components of the vertical diffusivity. Shear-driven diffusivity (Jackson et al., 2008;top left), internal-tide-driven diffusivity (Melet et al., 2012;top right), bottom boundary layer mixing (Legg et al., 2006;bottom left), surface boundary layer (Reichl & Hallberg, 2018;bottom right). All quantities are averages of years 1988-2007.

10.1029/2019MS001726
Journal of Advances in Modeling Earth Systems related to the wind stress. Eddy saturation in the ACC is the process whereby an increase in zonal winds manifests as an increase in eddy energy yet without a corresponding increase in net zonal transport (Farneti et al., 2010Hallberg & Gnanadesikan, 2006;Hogg et al., 2015;Meredith et al., 2011;Straub, 1993). We here exhibit diagnostics that measure the degree to which OM4 simulations exist in an eddy saturated regime. For this purpose, it is sufficient to focus on mass transport through the Drake Passage since trends in transport are quasiuniform longitudinally around the ACC. 3.8.1.1. Time Series of the Transport In Figure 17 we exhibit the time series of depth-integrated mass transport through the Drake Passage for the five CORE-II cycles. OM4p25 and OM4p5 are notable for reaching a quasi-repeating behavior after the third cycle, whereas OM4p5n continues to drift. We ran OM4p5n for another five CORE-II cycles (10 cycles total) with continued drift throughout. The transports in OM4p25 and OM4p5 are within the central portion of the transport found in other CORE-II models analyzed by Farneti et al. (2015, see their Figures 5 and 6). Donohue et al. (2016) estimate the observed transport at 173±9Sv, thus indicating that OM4p5 is well within that range, OM4p25 is slightly below, and OM4p5n is well above. Donohue et al. (2016) estimate a sizable contribution to zonal mass transport from a relatively strong bottom flow. Assuming a constant velocity throughout the water column equal to the bottommost measured velocity, they estimate that this flow transports 46 ± 9 Sv through the Drake Passage, which corresponds to an averaged bottom zonal velocity of 1.3 cm/s. The remainder of the transport (roughly 127Sv) arises from thermal wind balanced flow diagnosed relative to the bottom. Figure 18 shows the time-averaged residual mean zonal velocity (Eulerian mean plus parameterized eddy induced) through the Drake Passage. Each simulation is dominated by eastward flow; however, the vertical and latitudinal structures differ. OM4p25 shows relatively fine latitudinally spaced currents, thus reflecting the presence of zonal density fronts and corresponding currents characteristic of the observed ACC (e.g., Dufour et al., 2015;Rintoul & Naveira Garabato, 2013). These currents extend into the abyss, although in some latitudes they reverse sign thus reflecting a first mode baroclinic structure (Peña-Molino et al., 2014). There is a particularly strong westward abyssal current on the northside of the Drake Passage, below and slightly southward of the strong eastward surface flow abutting South America. OM4p5n simulation exhibits a relatively large drift over the first 100 years, whereas OM4p5 and OM4p25 repeat their transport after roughly three cycles, with the start of each new cycle signaled by a reduction in transport. The transport in OM4p5 and OM4p25 broadly increases over the course of a 60-year cycle though with OM4p5 showing a more pronounced increase than OM4p25 (see Figure 19). Observational measurements analyzed by Donohue et al. (2016) place the net transport at 173±9 Sv. The latitudinal average bottom zonal velocity in OM4p25 (not shown) is roughly an order of magnitude less than the Donohue et al. (2016) estimates. Hence, only about 5 Sv can be attributed to barotropic transport so that the corresponding zonal transport in OM4p25 is dominated by thermal wind flow. The other two models also exhibit relatively weak bottom transport, along with broad latitudinal structures reflective of coarsely resolved ACC jets. Notably, OM4p5 exhibits eastward flow throughout the depth with a weak flow near the bottom (equivalent-barotropic structure as in Killworth, 1992). OM4p5n shows stronger flow into the abyss; however, it exhibits a slight reversal of flow (to the west) at some latitudes. The rather large transport in OM4p5n arises from both the relatively large baroclinicity in this model, with isopycnals generally sloping much more steeply than in OM4p25 and OM4p5 (see Figure 23), as well as the relatively strong flow throughout the depth. Mesoscale eddies, either represented or parameterized, generally reduce the available potential energy and thus reduce density slopes (Gent et al., 1995). The absence of a mesoscale closure in OM4p5n thus allows for its density surfaces to steeply slope more than in OM4p5 (with mesoscale closure) and OM4p25 (with explicit eddies, albeit incompletely resolved).

Signatures of Eddy Saturation
The very large transport found in OM4p5n, as well as its large amplitude unsteady behavior, motivate us to focus just on OM4p5 and OM4p25 for a characterization of eddy saturation. For that purpose, we display in Figure 19 the Drake Passage time series for OM4p25 and OM4p5 during the fifth CORE-II cycle, excluding the first 10 years of each cycle to avoid contamination from 2007 forcing transitioned back to 1948. A linear fit indicates that OM4p5 increases its transport at a rate more than twice that of OM4p25 during years 1958-2007. We conclude that OM4p5, with its parameterized mesoscale eddy transport, exhibits less eddy saturation than OM4p25, with its explicit (though incompletely represented) mesoscale eddies. We compare this result to the Farneti et al. (2015) study of 17 CORE-II simulations. They found that models such as OM4p25, with 0.25°horizontal grid spacing and no mesoscale eddy parameterization, generally exhibit more eddy saturation than models making use of a mesoscale eddy parameterization where the eddy diffusivity is depth independent. OM4p5, however, makes use of a depth-dependent diffusivity as determined by an equivalent barotropic structure (see section 2.2.6). Further examination of the reasons for differing eddy saturation behavior between OM4p25 and OM4p5 is the topic of ongoing research.

Southern Ocean Meridional Overturning Circulation
We now consider the Southern Ocean meridional overturning circulation and the response of this circulation to the positive trend in the Southern Annular Mode corresponding to an increase in the eastward zonal wind stress. In Figure 20 we show the Southern Ocean meridional overturning circulation as a function of depth and as a function of potential density referenced to 2,000 dbar. For all simulations we include contributions from the submesoscale eddy parameterization (active just in the mixed layer), and for the OM4p5 simulation  we also include contributions from the mesoscale eddy parameterization. Note the qualitatively similar behavior found in OM4p25 and OM4p5 whereas OM4p5n is a clear outlier.

Description of the Overturning Cells
As a function of depth, the OM4p5n overturning possesses a stronger upper cell (Deacon cell) than the other configurations. We suggest that this behavior arises from the deeper penetration of the wind driven circulation due to the weaker vertical stratification in OM4p5n. Note that part of the upper cell strength in depth coordinates results from meridional transport along zonally sloped isopycnals that cross depth surfaces. Hence, when mapping this "Deacon cell" circulation in density coordinates, its strength greatly reduces (Döös & Webb, 1994;Karoly et al., 1997).
There is a relatively strong abyssal and deep overturning in OM4p5n (the blue counter-clockwise circulation centered at 60°S) when viewed in density space. This cell is affected by deep reaching convection as revealed by extremely deep mixed layers in OM4p5n ( Figure 21). Furthermore, the deep mixed layers give rise to the large north-south density difference in turn giving rise to the strong thermal wind transport in the ACC. It is striking how the mesoscale eddy parameterization in OM4p5 greatly reduces the mixed layer depth relative to OM4p5n, and thus affects the overall structure of the meridional overturning (Danabasoglu et al., 1994 first noted this point). Note that this strong abyssal overturning is largely absent from the depth space overturning since the overturning is associated with weak net vertical motion when integrated around a latitude circle. Furthermore, as seen by the mixed layer depths in Figure 21, this deep overturning in OM4p5n is largely contained within the horizontal Weddell Gyre circulation.
Following Karoly et al. (1997) and McDougall and McIntosh (2001), we exhibit in Figure 20 the density overturning circulation mapped as a function of the zonal mean and time mean depth of the density surface, thus providing a means to determine the depth range over which the density-space overturning occurs. For each simulation, the bottom cell (which includes the Antarctic Bottom Water (AABW) cell) in the far south appears over a large depth range, reflecting the predominantly open ocean source for this part of the overturning circulation. The OM4p5n simulation shows a particularly strong bottom cell that corresponds to the widespread open ocean deep convection as revealed in the mixed layer (Figure 21), which in turn appears to suppress the upper cell strength. In contrast, OM4p25 and OM4p5 show a more modest amount of AABW formation.

Signatures of Eddy Compensation
Following Marshall & Radko (2003; see also Appendix A of Farneti et al., 2015 for a succinct review), we conceive of the upper cell of the meridional overturning circulation in the Southern Ocean (i.e., the red clockwise cell in the density overturning in Figure 20) as the sum of an Eulerian mean circulation, Ψ , and eddy-induced circulation, Ψ * At latitudes and depth range of the Drake Passage, the Eulerian mean is directly driven by the zonal mean zonal wind stress, τ (3) where f is the Coriolis parameter. The Eulerian mean thus represents the wind driven upwelling of deep Southern Ocean waters into the ACC. As the wind stress increases so too does the Eulerian circulation, acting to increase the store of available potential energy that in turn enhances the mesoscale eddy activity feeding off the available potential energy reservoir. This eddy activity in turn leads to the counterclockwise eddy-induced circulation that in part counteracts the Eulerian overturning. This process is known as eddy compensation.
As for the depth integrated ACC transport, we wish to characterize the response of the overturning circulation to increases in the zonal winds found in the CORE-II forcing. If the models are "fully eddy compensated," the increases in winds will increase the Eulerian overturning, but that increased Eulerian overturning will in turn lead to stronger eddies thus leaving the net overturning roughly unchanged. Following Farneti et al. (2015, see their Figure 28), we determine the degree of eddy compensation upon changes to the wind stress by comparing the time series of the maximum Southern Ocean overturning in the upper cell (measured in density space) south of 30°S. Figure 22 shows the resulting time series for OM4p25 and OM4p5 along with a linear fit. As for the Drake Passage transport shown in Figure 19, we here see that OM4p5 has roughly twice the response as OM4p25, thus indicating that meridional overturning in OM4p25 exhibits a more eddy compensated response than OM4p5. These results agree with earlier work showing that, in response to westerly wind intensification, the degree of eddy compensation increases with the horizontal resolution of the model (e.g., Dufour et al., 2012;Morrison & Hogg, 2012).

Potential Temperature and Salinity Stratification
In the vicinity of Antarctica, large amounts of heat stored in the ocean subsurface can potentially access the atmosphere and the ice shelves, with important implications for sea ice cover and ice shelf melting (e.g., Jenkins et al., 2016;Rignot & Jacobs, 2002). Model representation of stratification and water mass properties in this region is thus key to climate projections. However, ocean models generally do poorly in the weakly stratified polar regions partly due to the complex coupling between atmosphere, ocean, sea ice and ice sheet. In particular, many models are subject to spurious deep convection sometimes associated with the opening of large open-ocean polynyas (e.g., de Lavergne et al., 2014;Heuzé et al., 2013). Some recent studies have suggested that increasing oceanic resolution could lead to improvement in model performance, in particular through a better representation of the mesoscale eddy field and/or dense water overflows (Dufour et al., 2017;Newsom et al., 2016). In this subsection, we briefly assess the OM4 models' skills in representing stratification in the vicinity of Antarctica. Figure 23 displays potential temperature and salinity cross sections in the Weddell Sea sector along with potential density surfaces, together highlighting the stratification in this region known to be one of the  [1948][1949][1950][1951][1952][1953][1954][1955][1956][1957] during the fifth CORE-II cycle (fourth cycle has similar relative trends). The red line is a linear fit to the transport, with OM4p5 having roughly twice the slope than OM4p25. We conclude that OM4p25 exhibits more eddy compensation behavior than OM4p5.
most weakly stratified of the Southern Ocean. Both OM4p25 and OM4p5 show a similar density structure with a large-scale doming of isopycnals centered about 65°S characteristic of the Weddell Sea gyre, transitioning north of 60°S to sloping isopycnals associated with the ACC region. The surfaces are slightly more tilted in OM4p5 (more baroclinicity) than in OM4p25, consistent with enhanced eddy saturation in OM4p25 as discussed in section 3.8.1.
OM4p5n is a clear outlier in Figure 23, exhibiting a relatively unstratified water column at 65°S with a chimney of intense deep convection consistent with deep mixed layers shown in Figure 21. Deep convection leads to the formation of AABW resulting in an intensification of the bottom cell of the overturning circulation (see Figure 20). Though deep convection in the Weddell Sea should lead to injecting the relatively warm subsurface waters toward the surface, this upward heat flux does not lead to the opening of a polynya in OM4p5n. Still, we note a particularly strong bias in sea ice concentration and extent in the Antarctic sea ice of OM4p5n, with anomalies located in the Weddell Sea that might be partly due to the representation of stratification (see section 3.9, Figures 25 and 26).
Both OM4p25 and OM4p5 show a significant drift in the density structure of the Weddell Sea sector from the initial state, revealed by a downward displacement of isopycnal contours. This drift is likely due to a combination of the absence of dense water resupply from the continental shelf, either due to a lack of dense water formation and/or representation of dense water overflow, and of spurious diapycnal mixing in the interior ocean leading to the consumption of AABW (e.g., Dufour et al., 2017;Lee et al., 2002). This drift is not apparent in OM4p5n only because intense deep convection occurs throughout the water column supplying the abyssal ocean with dense waters. OM4p25 does not show significant improvement on the representation of stratification in the deep ocean relative to OM4p5. This result suggests that, despite the use of hybrid coordinates, the OM4 models still lack an adequate representation of overflows (see also section 3.7). This conclusion motivates further work on the representation of shelf processes and dense water overflows around Antarctica.

Arctic and Antarctic Sea Ice
We here evaluate the simulated climatology and interannual variability of Northern and Southern Hemisphere sea ice from the three OM4 simulations OM4p25, OM4p5, and OM4p5n. The simulated sea ice concentration (SIC) and sea ice extent (SIE) are assessed against monthly averaged passive microwave satellite SIC observations from the National Snow and Ice Data Center (NSIDC) processed using the NASA Team Algorithm (data set ID: NSIDC-0051, Cavalieri et al., 1996). Pan-Arctic and Pan-Antarctic SIE are defined respectively as the areal sum of all grid points whose SIC exceeds 15% in the Northern and Southern Hemispheres. All climatologies are computed over the years 1979-2007, which is the common period in which both satellite SIC observations and CORE IAF exist. The analyses shown here are based on forcing cycle five of the CORE IAF.

Sea Ice Climatology
The magnitude and timing of the seasonal cycle of Pan-Arctic SIE is reproduced with quite good fidelity in each of the OM4 CORE runs (see Figure 24a). Each of these runs has a modest low SIE bias in boreal winter and a more notable high SIE bias in boreal summer. These SIE biases are relatively small compared with analogous CORE IAF simulations performed with other ice-ocean models (compare with Figure 5 of Wang et al., 2016). The climatological spatial patterns of simulated March SIC generally display positive SIC biases (too much sea ice) in the Greenland and Barents Seas and negative biases in the Labrador Sea and Sea of Okhotsk (see Figures 25a-25c).
The Arctic winter SIC biases are similar across the three experiments with the exception of the OM4p25 run, which has a reduced bias in the Sea of Okhotsk. This result indicates that both resolved and parameterized eddy heat and momentum fuxes have a relatively minimal influence on the Arctic winter sea ice edge in this model. The three experiments also share common Arctic summer SIC biases, with a notable positive bias in the Canadian Arctic Archipelago and positive biases along the climatological ice edge within the Central Arctic (see Figures 25d-25f). This positive bias in the Canadian Arctic Archipelago is a common feature in ice-ocean models forced by CORE IAF (see Figure 4 of Wang et al., 2016), possibly related to a negative bias in downwelling shortwave radiation over this region in the CORE forcing (Wang et al., 2012).
The model biases and their sensitivity to ocean configuration are more significant in the Southern Hemisphere. Overall, we find a high degree of consistency between the OM4p25 and OM4p5 simulations. These experiments both have positive Antarctic winter SIE biases, owing primarily to positive SIC biases in the Western Pacific and Indian Ocean sectors and the Amundsen and Bellingshausen Seas (see Figure 24b and Figures 25g-25i). The Antarctic summer SIE biases in these simulations are smaller, but this is due to canceling SIC biases, as these experiments have too much sea ice in the Eastern Weddell Sea and too little ice in the Western Weddell Sea and around most of the Antarctic coastline (Figures 25k-25l). The timing of the seasonal cycle of Antarctic SIE is also biased in these runs: The observed SIE reaches clear minima and maxima in February and September, respectively, whereas the simulated SIE has minimum values in March and maximum values in September/October. The winter and summer Antarctic SIE biases in the OM4p25 and OM4p5 runs have similar magnitude to other ice-ocean models driven by CORE IAF . Also, the negative summer SIC bias in the Western Weddell Sea is a ubiquitous feature across CORE IAF-forced simulations (seen in each of 15 models studied by Downes et al., 2015) and the  positive winter SIC bias in the Australian sector (90°E to 150°E) is also a relatively common feature (present in 11 of the 15 models).
The OM4p5n run with no eddy parameterization is significantly different from the OM4p25 and OM4p5 simulations, indicating that mesoscale eddies have a substantial impact on Antarctic sea ice. The OM4p5n run has negative Antarctic SIE biases in all months of the year and SIC biases that are generally negative with the exception of positive winter anomalies in the West Pacific sector. These sea ice differences are consistent with the lack of restratification due to eddies and the corresponding deeper Southern Ocean mixed layers in the OM4p5n run, which tends to reduce winter sea ice growth (see Figure 21). The similarity between the OM4p25 and OM4p5 experiments suggests that the mesoscale eddy parameterization used in OM4p5 is capturing the effects of eddies on Southern Ocean stratification and sea ice with good fidelity.

Interannual Variability and Trends of Sea Ice
We next consider the ability of the simulations to capture the observed trends and interannual variability in Arctic and Antarctic sea ice. In Figure 26, we plot Arctic and Antarctic SIE time series for March and September from the three CORE experiments and NSIDC observations. All three experiments capture the observed trends and interannual variability in Arctic SIE with reasonably good fidelity. The detrended correlations for September SIE range from 0.77-0.81 and the detrended correlations for March SIE range from 0.71-0.83 (see summary statistics in Table 3). The September and March correlations from the OM4p25 run (0.81 and 0.83, respectively) are modestly higher than the detrended correlations from the 14 ice-ocean models presented in Wang et al. (2016), whose correlation values range from 0.60-0.80 and 0.37-0.80, respectively. These experiments also capture the magnitude of interannual variations well, as the detrended standard deviations agree closely with observed values. The negative Arctic SIE trends in each of these experiments are underestimated in both September and March by roughly 30%. The underestimation of the March SIE trend was a common feature across 13 of the 14 ice-ocean models examined by Wang et al. (2016), whereas the September SIE trends exhibited larger intermodel spread.
The interannual variability of Antarctic sea ice is not captured well by these simulations, with detrended correlation values ranging from 0.31-0.38 and −0.11-0.45 for March and September SIE, respectively. These low correlation values are potentially attributable to lower-quality forcing data in the Southern Ocean due to sparse observations in this region, but this would need to be examined using a larger suite of ice-ocean models. The OM4p25 and OM4p5 experiments have March sea ice variability that is too strong and September sea ice variability that agrees more closely with observations. The OM4p5n experiment displays qualitatively different September sea ice variability, characterized by large oscillations with a roughly 20-year period. The OM4p5n September SIE has more variability on interannual time scales, with a magnitude similar to observations. The OM4p25 and OM4p5 experiments have negative September SIE trends, failing to capture the slightly positive observed trend. Overall, these results show that Arctic sea ice trends and variability are simulated with much higher fidelity than Antarctic trends and variability.

Closing Comments
We have described two main configurations of a new global ocean/sea ice model at an eddy-permitting 1/4°resolution (OM4p25) and a noneddying 1/2°-resolution (OM4p5), along with some variants to illustrate key sensitivities that influence the design of the main models. One of the earlier decisions during the development was to adopt a hybrid isopycnal-z* vertical coordinate instead of the simpler z* coordinate. While attempting to build the z* coordinate model, we were unable to reduce the explicit vertical mixing enough to compensate for a large oceanic heat uptake (not shown). Here, we demonstrated a significant reduction in heat uptake by the ocean, of order 0.25°C/century, simply by switching from the z* coordinate to the hybrid isopycnal-z* coordinate. This result is in full agreement with other studies about the advantages of using Lagrangian vertical coordinates (e.g., Chassignet et al., 2003;Megann et al., 2010Megann et al., , 2018. We compared the role of resolution in setting the physical fidelity of the solutions. The mean state of the coarser OM4p5 model is quite comparable to the finer OM4p25 model due to particular choices of the mesoscale eddy parameterizations. We also adjusted the submesoscale parameterization of mixed-layer instabilities between the two resolutions which had a mild effect on MLDs. Although the mean states of the two models are comparable, we demonstrate that the eddy-permitting model exhibits more eddy saturation/compensation behavior in the Southern Ocean than when the eddies are parameterized in the coarser model (see Farneti et al., 2015 for a model intercomparison of this behavior).
There remain several large-scale biases in these models that require further work. For example, there is a systematic shallow mixed-layer bias, and a shallow meridional overturning resulting from weak Nordic Sea overflows. For some of these biases, such as dilution of the Nordic Sea overflows by excessive mixing, we have hypotheses for how to remedy them in future development. However, other biases remain elusive. Whether these biases are unique to this model configuration or are more generic can be answered through model comparison studies such as the ongoing CMIP6/OMIP (Griffies et al., 2016).
There is a tendency for model developers to focus on the problematic aspects of a model simulation since it is the problems that developers aim to remedy. Despite the highlighted biases detailed in this paper, the two main models described here are remarkably faithful to the observations relative to earlier generation models. For example, the interior and volume-integrated drifts are small compared to other models using the same atmospheric state as evaluated by Griffies et al. (2009) and Danabasoglu et al. (2014). Development of our new models started from a new code base, thus motivating us to reevaluate every decision such as physical parameterizations and representation of the vertical coordinate. While some of the remaining biases may be due to missing parameterizations or nonoptimal choice for the hybrid vertical coordinate, we are nonetheless pleased to have produced a respectable model using fewer parameterizations than previous generations of GFDL ocean models.
of planetary vorticity f ¼ 2Ωsinϕ, the constant Boussinesq reference density ρ 0 , and the gravitational potential Φ=gz, In the horizontal momentum equations F represents the accelerations due to the divergence of stresses including those provided through boundary interactions.
The prognostic thermodynamic variables are potential temperature, θ, and salinity S, which are related to in situ density ρ through the Wright (1997) equation of state. In the potential temperature and salinity equations, fluxes due to diabatic, vertically oriented processes are indicated by J (z) . The tendency due to the convergence of fluxes oriented along neutral directions is indicated by N γ . Our implementation of this neutral diffusion parameterization is detailed in Shao et al. (personal comm.).
Transforming these equations to general vertical coordinate r=r(x,y,z,t) (Kasahara, 1974;Starr, 1945) yields The time derivatives are now computed with the generalized vertical coordinate fixed rather than the geopotential. We introduced the specific thickness, z r =∂z/∂r, which measures the inverse vertical stratification of the vertical coordinate surfaces.
Similar to Bleck (2002), MOM6 is discretized in the vertical by integrating between surfaces of r to yield layer equations where the layer thickness is h ¼ ∫z r dr and variables are treated as finite volume averages over each layer: where δ r =dr (∂/∂r) is the discrete vertical difference operator. The pressure gradient accelerations in the momentum equation A13 are written in continuous-in-the-vertical form for brevity; the exact discretization is detailed in Adcroft et al. (2008). The MOM6 time-stepping algorithm integrates the above layer-averaged equations forward allowing the vertical grid to follow the motion, that is, _ r ¼ 0, so that the red terms are dropped. This approach is generally known as the Lagrangian method, but here the Lagrangian method is used only in the vertical direction. After each Lagrangian step, a remap step is applied that generates a new vertical grid of the user's choosing. The ocean state is then mapped from the old to the new grid. The physical state is not meant to change during the remap step, yet truncation errors make remapping imperfect. We employ high-order accurate reconstructions to minimize errors introduced during the remap step (White & Adcroft, 2008;White et al., 2009). The connection between time stepping and remapping are described in section A.3.

A.2 Specifics of the Ocean Model Equations
We here provide more details for the terms appearing in the ocean model equations described in section A.1.

A.2.1 Horizontal Momentum Equation
Equation A13 is the horizontal momentum equation written in its vector-invariant advective form with u ¼ b x u þ b y v the horizontal velocity, p the hydrostatic pressure, f the Coriolis parameter, and the vertical component of the vorticity using ∇ r for the curl operator. The discretization of the Coriolis term is the enstrophy conserving scheme of Sadourny (1975). The geopotential coordinate, z, has a value z=0 at a resting ocean surface, z=η(x,y,t) at the ocean free surface, and z=−H(x,y) at the ocean bottom. We use the Boussinesq approximation (volume conserving kinematics) with ρ 0 = 1,035 kg/m 3 the reference density. MOM6 has an option for compressible non-Boussinesq flow [mass conserving kinematics]. However, we chose the Boussinesq option largely based on legacy. Time and horizontal derivatives are computed holding the generalized vertical coordinate fixed rather than the geopotential The transport of seawater crossing surfaces of constant r is measured by the dia-surface velocity component (see section 6.7 of Griffies, 2004) ∂z ∂r with z r the specific thickness that is assumed one-signed throughout the ocean, and D/Dt the material time derivative operator. In the ocean interior where r is aligned with isopycnals, the dia-surface velocity becomes the diapycnal velocity whose value is directly related to irreversible processes such as mixing that act on potential temperature and salinity. In the unstratified mixed layers, r=z * so that z r _ r ¼ ð∂z=∂z * ÞDz * =Dt , which is close to the familiar vertical velocity component Dz/Dt. Viscous dissipation (Laplacian and biharmonic friction following  and mechanical boundary forces (winds, bottom stress) contribute to the divergence of the deviatoric (symmetric and tracefree) stress tensor, F ¼ ∇·τ. MOM6 and the real ocean have no vertical sidewalls, and MOM6 treats all solidearth boundaries with bottom stress parameterized as a quadratic drag. A.2.2 Hydrostatic Balance Equation (A14) is the discrete version of the hydrostatic balance. The horizontal pressure gradient force is implemented as a contact force following the method of Adcroft et al. (2008). These equations differ from Bleck (2002) who uses the Montgomery potential to calculate pressure gradient accelerations.

A.2.3 Thickness and Tracer Equations
Volume conservation appears in the form of a prognostic flux-form layer thickness equation A15, with the nonnegative layer thickness given by where dr is the thickness of a layer in r−space (e.g., the density difference between target density classes or the thickness between target depths). The layer thickness increases where horizontal thickness fluxes converge, ∇ r · h k u ð Þ<0, and where dia-surface flow converges, δ r ðz r _ rÞ<0. The volume flux h k u is computed using the quasi-third-order PPM scheme (Colella & Woodward, 1984) using a positive-definite limiter rather than the monotonic limiter. This last choice avoids limiting of positive extrema and thus retains third-order accuracy everywhere except near vanishing layers. Transport in the thickness equation is discretized compatibly with that in the flux-form potential temperature and salinity equations A16 and (A17). Compatibility is required to maintain global and layer integrated conservation properties for volume, heat, and salt. Tracer reconstruction for transport uses PPM with monotonic limiters but using third-order interpolation for edge values. This reduces the size of the stencil which helps the computational efficiency of the transport scheme. The flux convergences, N γ θ and N γ S , provide subgrid-scale neutral diffusion for the potential temperature and salinity, whereas δ r J ðzÞ θ and δ r J ðzÞ S provide subgrid-scale vertical diffusion as well as boundary fluxes. In the interior, both subgrid fluxes vanish when their respective tracers are spatially uniform, thus ensuring that the tracer equation reduces to the thickness equation when the tracer is uniform.
Parameterized subgrid-scale advection from the submesoscale (Fox-Kemper et al., 2011) and mesoscale (Gent et al., 1995) parameterizations are combined with the lateral advection of thickness and tracer, thus providing a residual mean advective transport for the scalar fields. Furthermore, we implement subgrid advective terms solely as lateral transports, thus interpreting them as layer bolus transport as appropriate for vertical Lagrangian models rather than a three-dimensional eddy-induced advection as appropriate for vertical Eulerian models (see McDougall & McIntosh, 2001 for details).

A.2.4 Equation of State
The equation of state, (A18), determines in situ density as a function of potential temperature, salinity, and pressure. We evaluate the pressure in the equation of state according to −g ρ 0 z. Doing so maintains energetic consistency for the Boussinesq fluid according to section 2.4.3 of Vallis (2017). We make use of the Wright (1997) equation of state so that θ is potential temperature and S is the practical salinity. Although MOM6 has the more updated equation of state from IOC et al. (2010), the required changes for thermodynamic variables were implemented only after the basic model configuration was developed. Time constraints on model development prompted us to retain usage of Wright (1997) for OM4.
The freezing point of seawater is approximated as T f ¼ −0:054S−7:75×10 −08 p; where p is in units of Pascals and S is in units of 1×10 −3 . When the local temperature anywhere in the ocean column falls below the freezing point, the water-equivalent volume of ice is calculated and the fusion heat locally added back to the ocean to raise the liquid seawater temperature back to the freezing point. The frozen water and salt are sent to the sea ice model.

A.3 Basics of the Vertical Lagrangian-Remap Method in MOM6
As discussed by Adcroft and Hallberg (2006), there are two general classes of algorithms that frame how hydrostatic ocean models are formulated. The two classes differ in how they treat the vertical direction. Quasi-Eulerian methods follow the approach traditionally used in geopotential coordinate models, whereby vertical motion is diagnosed via the continuity equation. Quasi-Lagrangian methods are traditionally used by layered isopycnal models, with the Lagrangian approach specifying motion that crosses coordinate surfaces. Indeed, such dia-surface flow can be set to zero using Lagrangian methods for studies of adiabatic dynamics. MOM6 makes use of the vertical Lagrangian remap method, as pioneered for ocean modeling by Bleck (2002), which is a limit case of the Arbitrary-Lagrangian-Eulerian method (Hirt et al., 1997). Dia-surface transport is implemented via a remapping so that the method can be summarized as the Lagrangian plus remap approach and is essentially a one-dimensional version of the incremental remapping of Dukowicz and Baumgardner (2000).
The MOM6 implementation of the vertical Lagrangian-remap method makes use of two general steps. The first evolves the ocean state forward in time according to a vertical Lagrangian limit with _ r ¼ 0. Hence, the horizontal momentum, thickness, and tracers are time stepped with the red terms removed in equations A13, (A15), (A16), and (A17). All advective transport thus occurs within a layer as defined by constant r-surfaces so that the volume within each layer is fixed. All other terms are retained in their full form, including subgrid scale terms that contribute to the transfer of tracer and momentum into distinct r layers (e.g., dia-surface diffusion of tracer and velocity). Maintaining constant volume within a layer yet allowing for tracers to move between layers engenders no inconsistency between tracer and thickness evolution. The reason is that tracer diffusion, even dia-surface diffusion, does not transfer volume.
The second step in the algorithm comprises the generation of a new vertical grid following a prescription, such as whether the grid should align with isopcynals or constant z * or a combination. The ocean state is then vertically remapped to the newly generated vertical grid. The remapping step incorporates dia-surface transfer of properties, with such transfer depending on the prescription given for the vertical grid generation. To minimize discretization errors and the associated spurious mixing, the remapping step makes use of the high-order accurate methods developed by White and Adcroft (2008) and White et al. (2009).
The underlying algorithm for treatment of the vertical can be related to operator splitting of the red terms in equations A15 and A16. If we consider, for simplicity, an Euler-forward update for a time-step Δt, the time stepping for the continuity and temperature equation can be summarized as  (A26) is essentially used to define the dia-surface transport z r _ r by prescribing h (n+1) . For example, to recover a z coordinate model, h (n+1) =Δz, and z r _ r becomes the Eulerian vertical velocity, w.
Within the above framework for evolving the ocean state, we make use of a standard split-explicit time stepping method by decomposing the horizontal momentum equation into its fast (depth integrated) and slow (deviation from depth integrated) components. Furthermore, we follow the methods of Hallberg and Adcroft (2009) to ensure that the free surface resulting from time stepping the depth integrated thickness equation (i.e., the free surface equation) is consistent with the sum of the thicknesses that result from time stepping the layer thickness equations for each of the discretized layers; that is, ∑ k h ¼ H þ η.
A.4 Data Sets OM4 requires the following external data sets as inputs: (i) bathymetry, (ii) geothermal heat flux through the sea floor, (iii) ocean color to regulate short-wave radiation penetration, (iv) tidal forcing fields and topographic roughness for internal wave mixing parameterization, (v) atmospheric state and river runoff to calculate the driving fluxes of heat, fresh water, and momentum, and (vi) initial conditions for potential temperature and salinity. Each of these data sets is preprocessed offline and read by the model at run-time. The details of pre-processing are no less important than those of the numerical model in determining the solutions. Hence, the scripts and steps used for each data set are version controlled and published with the model code. A.4.1 Bathymetry Model bathymetry was derived from a blend of the GEBCO 30-s gridded topography (Weatherall et al., 2015) and International Bathymetric Chart of the Arctic Ocean (Jakobsson et al., 2012) using simple averaging to the model grid. Additional hand edits were applied to reopen closed straits or replace missing barriers and islands. The list of all modifications are recorded in the repository and stored as auxiliary data with the model input files. Topographic roughness estimates were computed by first performing a local planar fit of the fine resolution points within each model grid cell and calculating the root-mean-square deviation about this surface. Roughness was used as part of the bottom-enhanced internal tide mixing scheme from section 2.2.3.

A.4.2 Geothermal Heat Flux
A geothermal heat flux is applied at the sea floor using a map from Davies (2013). The data are interpolated to the nodes of the super grid (quarter cells) of each model using bilinear interpolation and then integrated to the model grid using the trapezoidal rule. This interpolation is not exactly conservative but the coarse resolution of the source data and relatively fine resolution of the model supergrids makes the errors small compared to the uncertainty in the data.

A.4.3 Ocean Color
A monthly climatology of chlorophyll from the NASA ocean color product SeaWIFS (doi:10.5067/ ORBVIEW-2/SEAWIFS_OC.2014.0), years 1997-2010, is used for calculating short-wave penetration. A.4.4 Atmosphere and River Data from CORE When we began development of OM4, the interannual CORE data set of Large and Yeager (2009) was the only standard approach used in the community for developing ocean/ice models. We thus chose to focus our assessment on CORE since it is more mature and there are a number of benchmark papers that examine 15-20 other ocean/sea ice models using CORE.

A.4.5 Initial Conditions and Hydrographic Data
The initial hydrographic state is based on the World Ocean Atlas January climatology Locarnini et al., 2006), which is interpolated to the model grid at runtime. The same data source is used for many analyses of biases and comparison.