Incorporation of the C-GOLDSTEIN efficient climate model into the GENIE framework: “eb go gs ” configurations of GENIE

Abstract. A computationally efficient, intermediate complexity ocean-atmosphere-sea ice model (C-GOLDSTEIN) has been incorporated into the Grid ENabled Integrated Earth system modelling (GENIE) framework. This involved decoupling of the three component modules that were re-coupled in a modular way, to allow replacement with alternatives and coupling of further components within the framework. The climate model described here (referred to as "eb_go_gs" for short) is the most basic version of GENIE in which atmosphere, ocean and sea ice all play an active role. Among improvements on the original C-GOLDSTEIN model, latitudinal grid resolution is generalized to allow a wider range of surface grids to be used. The ocean, atmosphere and sea-ice components of the "eb_go_gs" configuration of GENIE are individually described, along with details of their coupling. The setup and results from simulations using four different meshes are presented. The four alternative meshes comprise the widely-used 36 × 36 equal-area-partitioning of the Earth surface with 16 depth layers in the ocean, a version in which horizontal and vertical resolution are doubled, a setup matching the horizontal resolution of the dynamic atmospheric component available in the GENIE framework, and a setup with enhanced resolution in high-latitude areas. Results are presented for a spin-up experiment with a baseline parameter set and wind forcing typically used for current studies in which "eb_go_gs" is coupled with the ocean biogeochemistry module of GENIE, as well as for an experiment with a modified parameter set, revised wind forcing, and additional cross-basin transport pathways (Indonesian and Bering Strait throughflows). The latter experiment is repeated with the four mesh variants, with common parameter settings throughout, except for time-step length. Selected state variables and diagnostics are compared in two regards: (i) between simulations at lowest resolution that are obtained with the baseline and modified configurations, predominantly in order to evaluate the revision of the wind forcing, the modification of some key parameters, and the effect of additional transport pathways across the Arctic Ocean and the Indonesian Archipelago; (ii) between simulations with the four meshes, in order to explore various effects of mesh choice.


Introduction
To comprehensively model climate change on multimillennial timescales, to jointly explore and statistically analyse large subsets of model parameter space of combined climate-system components, or to run a multitude of long-term scenarios for anthropogenic climate forcing, we presently need to use "Earth System Models" (ESMs) which include all the necessary components and processes, and which are of high computational efficiency. Such ESMs must specifically represent the atmosphere, ocean, sea-ice, and, for some purposes ice sheets and biogeochemical cycles. For long-timescale palaeoclimate studies of Earth system dynamics, over several millions of years, changes in geography, orography and bathymetry must also be accommodated. These models must be computationally much faster than the general circulation models that are used to simulate climate change on shorter (centennial) timescales, to enable model integration over long timescales, and/or to accomplish large ensembles. A new generation of such ESMs of Intermediate Complexity (EMICs) has emerged in recent years (Claussen et al., 2002). In EMICs, computational speed is achieved through a combination of simplified physics, low resolution and efficient numerics (e.g., asynchronous coupling, implicit time-stepping). EMICs have been used operationally to explore long-term climate change and commitment (IPCC, 2007a), and for a variety of palaeoclimate studies (e.g., IPCC, 2007b). In one such application, Plattner et al. (2008) found no systematic difference in predictions between EMICs and higher-complexity models.
As a contribution to the EMIC class of models, the Grid-ENabled Integrated Earth system modelling (GENIE) framework is both flexible and scaleable in complexity (Lenton et al. ( , 2007; http://www.genie.ac.uk). Further advantages of the GENIE structure include the availability of grid-based tools for launching and analysing experiments (Lenton et al., 2007), and the use of a version-control infrastructure and the automatic routine testing of the code and any modifications on a daily basis. Further information on the general GENIE infrastructure can be found elsewhere (see http://source.ggy.bris.ac.uk/wiki/GENIE).
At the present time, the GENIE framework supports two main levels of complexity (dimensionality and physics) in the atmosphere, GENIE-1 and GENIE-2 (see Lenton et al., 2007). The present paper describes a basic version of GENIE-1, originally known as C-GOLDSTEIN (Edwards and Marsh (2005), henceforth EM05). This configuration of GENIE-1 is commonly referred to by the shorthand "eb go gs", where "eb" refers to the Energy and Moisture Balance Model (EMBM) of the atmosphere (by analogy, "ig" is commonly used to refer to the IGCM, the fullydynamical atmospheric model component available in the GENIE framework), "go" refers to the Global Ocean-Linear Drag Salt and Temperature Equation Integrator (GOLD-STEIN) ocean model, and "gs" refers to the sea ice scheme used in conjuction with GOLDSTEIN.
The purpose of this paper is to provide a basic description of eb go gs, emphasizing recent updates on earlier versions (e.g., Lenton et al., 2007). The standard version of eb go gs described in this paper is based on GENIE version 2.7.4 (tagged as "rel-2-7-4" in the source code repository). Adjustments to the code introduced after the release of GENIE version 2.7.4 were required for one of the alternative mesh options. GENIE version 2.7.7 has been verified to reproduce the results when set up using the configurations presented herein and the parameter and topography-related configuration files from the present study are planned (at the time of this writing) to be included in an upcoming GENIE release version (planned to be GENIE version 2.7.8, and to be tagged as "rel-2-7-8"). A further purpose is the description and basic evaluation of a range of meshes which can now be accomodated in the eb go gs model, motivated by a range of scientific foci and requirements, which include very high computational efficiency of the model, previously established model configurations, increased zonal or meridional resolution to better accomodate regional geographic or topographic features or localised processes, or compatibility with other model components.
We outline, and present results from, a baseline setup of the eb go gs model, which is currently in wide use as a standard configuration in conjunction with the simulation of biogeochemical cycles in GENIE, employing the BIOGEM and ATCHEM (Ridgwell et al., 2007a) (and often also the SEDGEM (Ridgwell and Hargreaves, 2007)) model components (henceforth referred to as the "biogem" configuration). This particular configuration of GENIE-1 was included as model "GENIE-16" (emphasizing 16 levels in the ocean) in the multi-model inter-comparison study of Cao et al. (2009). We further outline three horizontal meshes as alternatives to the widely-used horizontal resolution of 36 by 36 equalarea grid cells for GENIE-1, and present selected results obtained in spin-up experiments with each. All experiments are carried out with code on the GENIE development branch "Marsh et al-2009-GMDD", which contains specific differences not available in the release version "rel-2-7-4" (which, except in the case of one of the alternative meshes, do not affect the baseline performance of the model) and the specific parameter and topography-related configuration files for the set-ups shown here (unless they have already been present). Some aspects of, and results from, the "biogem" configuration used here are reported elsewhere (e.g., Cao et al., 2009). Variants of eb go gs in the "rel-2-7-4" release version have been subject to tuning exercises, results of which will be reported elsewhere.
The paper is organized as follows. In Sect. 2, we describe the model, with four sub-sections that describe the individual model components, coupling, time-stepping and meshes, and model parameters and inputs. In Sect. 3, we present the results obtained from default spin-up model runs, using a standard mesh (in two configurations), and selected results obtained with the three alternative meshes. In Sect. 4, we discuss the wider use of GENIE, emphasizing that the eb go gs setup underpins more complete model representations of the Earth System, already used to investigate phenomena such as abrupt climate change and biogeochemical cycles, in studies of past or future Earth System dynamics.

Model components
The three components of eb go gs -the frictionalgeostrophic GOLDSTEIN ocean model, an EMBM of the atmosphere, and a thermodynamic-dynamic model of sea ice -are largely based on the respective model components of the C-GOLDSTEIN climate model of EM05, and are each described in turn below. In addition to the three modules, run-time interpolation of the wind forcing fields has been implemented as an independent model component ("geniewind"), with potential future extension of this module in mind. The current functionality of this module is limited to the interpolation of externally provided wind and wind-stress datasets during model initialization, and hence no specific description of this model component is provided in this section (see Sect. 2.4 for further details).

Ocean
GOLDSTEIN is a 3-D frictional-geostrophic ocean model based on the reduced physics of the thermocline equations, described for a single-basin configuration in Edwards et al. (1998). A global ocean version which constitutes a precursor to the ocean component of C-GOLDSTEIN was introduced in Edwards and Shepherd (2002) with idealized global geography and resolution intermediate between box models and models with an approximation of realistic sub-basinscale topographic features. The presently-decribed version of GOLDSTEIN is essentially as described in EM05, with a few modifications, outlined below. Various further additions or improvements to GOLDSTEIN include: a surface mixed layer scheme (based on the scheme by Kraus and Turner (1967)); a bottom-boundary layer scheme (based on the scheme by Killworth and Edwards (1999)); an overflow scheme using the parameterisation of Danabasoglu et al., 2010; revised surface boundary conditions replacing virtual salinity and tracer fluxes with volume fluxes (following Huang (1993)); a more general equation of state through the parameterisation of thermobaricity (K. Oliver, personel communication, 2008); an enhanced diapycnal mixing scheme (Oliver and Edwards, 2008); the additional convection step used by Müller et al. (2006) preceding the conventional convection scheme. Incorporation of these improvements into the GENIE framework are already available as options or are currently in progress, and have been described elsewhere or are anticipated to be documented in future publications.
The prognostic variables at each model grid point are temperature and salinity (actually the difference relative to a reference salinity, default value 34.9), as well as the threedimensional velocity field. In summary, the following physical processes are explicitly represented: -Horizontal and vertical transport of heat, salinity, and other specified tracers (for example, tracers of biogeochemical cycles, which are transported as passive tracers within GOLDSTEIN), through advection, convection and mixing (comprising the combined parameterisation for isoneutral mixing and eddy-induced advection (Griffies, 1998), which can optionally be reverted to basic horizontal and vertical mixing); -Surface exchange of heat with the atmosphere and sea ice, provided by the EMBM coupling scheme (see Sect. 2.2), or through coupling with an alternative atmospheric component, is applied as a surface boundary condition; -Surface exchange of moisture with the atmosphere, sea ice and land, provided (as for heat exchange) by the EMBM coupling scheme (see Sect. 2.2); -Surface input of momentum (zonal and meridional wind  stress at locations where velocities are defined), provided by the new model component "genie-wind" (although technically a model component, its current sole  purpose is to interpolate prescribed wind speed and  wind stress fields onto the configured model grid during model initialization). Alternatively, if 'genie-wind' is deactivated, the wind stress fields are obtained from the EMBM, where they can be read in pre-interpolated form (this option is used in the simulation "biogem" described below).
We have incorporated three separate modifications relative to the original GOLDSTEIN model of EM05, two of which are options which generalize the application to different grids and geographies, while one represents a minor adjustment of the dynamics: -We have refined the calculation of barotropic flow to permit an arbitrary number of islands (disconnected land-masses). As noted in EM05, barotropic flow around islands, and hence through straits, can be calculated from the solution of a set of linear constraints arising from the integration of the depth-averaged momentum equations around each island. The standard C-GOLDSTEIN model setup of EM05 evaluates barotropic flow around Antarctica as the single contribution arising from such constraints. While EM05 explored the effect of additionally-enabled barotropic transport along the Bering Strait in a sensitivity simulation, no transport across the Indonesian Archipelago was enabled. Barotropic flows through Bering Strait and the Indonesian Throughflow (ITF) have significant effects on the modern global circulation (Hirst and Godfrey, 1993;Hu and Meehl, 2005), and the modification of the scheme for the solution of the barotropic streamfunction in GOLDSTEIN enables the representation of both flows in the coarse model grids used for the present study. Two of the model geographies introduced here feature the additional inclusion of barotropic flow around the North Pole (the North Pole is implicitly represented by a meridional no-flow condition at the northern model grid boundary), and an additional contribution to the barotropic streamfunction around Greenland is included, enabling net transport through the Davis Strait.
-A wider variety of horizontal mesh structure and resolution can now be accommodated, by generalizing the model equations to allow an arbitrary latitudinal distribution of grid point rows. The original code assumed equal spacing in sine of latitude, which, in conjuction with equal spacing in the longitude, results in equalarea cells. Three different mesh options are currently available: (i) constant spacing in the sine of latitude (parameter "igrid" set to 0); (ii) constant spacing in degrees 960 R. Marsh et al.: Incorporation of the C-GOLDSTEIN efficient climate model into the GENIE framework ("igrid" set to 1); (iii) a copy of the latitudinal axis of the IGCM atmospheric model component ("igrid" set to 2).
-An adjustment has been made to the variable upstreamweighting coefficients in the advection scheme that originally did not consistently take account of the grid curvature. The upstream weighting is an approximation of the Fiadeiro and Veronis (1977) scheme -see Müller et al. (2006) for further details.
In addition to these modifications, a scaling error in the non-dimensionalization of coefficients in the equation of state, which remained at the same "hard-wired" values when the default depth scale was updated from 4000 m to 5000 m, was discovered and corrected in June 2006 (Ridgwell et al., 2007a). The overall effect of this correction was a 25 % increase in the dimensional values of wind, currents, overturning, diffusivities and inverse timescales.
Except for the tracers transported by GOLDSTEIN, all quantities appearing in the model equations are in nondimensionalized form. With the exception of the depth scale, dimensional scale values and constants specific to the ocean component are as previously reported (Edwards et al., 1998;Edwards and Shepherd, 2001), and are listed here in Table A1.

Atmosphere
The Energy Moisture Balance Model (EMBM) of the atmosphere is based closely on that of the UVic Earth System Model (Weaver et al., 2001), essentially as described in EM05, with two exceptions (see below). Further enhancements of the EMBM include features that become available when the atmospheric component is coupled with the Efficient Numerical Terrestrial Scheme (ENTS) land model component, as well as forcing mechanisms related to prescribed ice sheets, both of which are beyond the scope of the present manuscript and have been described by Williamson et al. (2006) and Holden et al. (2009) respectively. The prognostic variables are air temperature and specific humidity in the atmospheric boundary layer. The heat and moisture balances of the atmosphere are outlined in Eqs. (1) and (2) of EM05.
In the EMBM, vertical fluxes are in SI units, while horizontal velocities and fluxes are non-dimensionalized, as for GOLDSTEIN. In contrast to the length of the time-step for updating the state of the ocean and sea-ice, as well as for evaluating heat and moisture fluxes between the three model components, the length of the time-step used in the EMBM transport equations is shorter (typically by a factor of five) and the transport equation is solved using an implicit scheme.
In summary, the following physical processes are explicitly represented in the EMBM module: -Horizontal transport of heat and moisture in the atmosphere by winds and mixing; -Precipitation, instantaneously removing all moisture corresponding to the excess above a relative humidity threshold (the threshold is configurable, but typically the default value of 85 % is used); -Surface exchange of heat with ocean, sea ice and land, according to balance of incoming shortwave radiation (a function of planetary and sea ice albedo), sensible and latent heat fluxes, re-emitted long-wave radiation, and outgoing planetary long-wave radiation (a function of radiative gases: CO 2 ; water vapour; CH 4 ; N 2 O); -Surface exchange of moisture with the ocean, sea ice and land, according to the balance of precipitation, evaporation/sublimation and runoff, the latter through instantaneous re-allocation of precipitation over land to coastal grid-boxes according to simple representation of catchments (see Sect. 2.3 and Fig. 1).
Two important changes have been incorporated since we originally described the EMBM in EM05. Firstly, the seasonal cycle of insolation is included by default; incoming shortwave radiation is computed for each latitude and point through the annual cycle at which exchange fluxes are evaluated, using modern day parameters for eccentricity, obliquity and precession. Secondly, the same flexibility of the horizontal mesh structure as allowed for the ocean module has been implemented. Scale values, constants, and configurable parameters specific to the atmosphere are listed in Table A2. In one further modification, the meridional profile of heat diffusivity at high southern latitudes may be scaled (reduced) southward of a specified latitude, according to a tuneable parameter, to improve the reprsentation of regional air temperatures and air-sea heat exchange.

Sea ice
We implement the sea ice model described in Weaver et al. (2001), details of which are provided in the Appendix of EM05. The model comprises thermodynamics plus advection with ocean current at the uppermost level of the ocean model. The prognostic variables are ice thickness, and ice areal fraction, or concentration. The transport of sea ice includes sources and sinks for these variables, computed by the coupling scheme that is technically part of the EMBM, hence this is where the thermodynamic sea-ice calculations are undertaken. As an update from the version used by EM05, the same horizontal mesh types can now be selected as for the EMBM and GOLDSTEIN ocean model components. Not specified in Weaver et al. (2001) is the formula of Millero (1978) (Eq. 9 of EM05), which we use for T f (S o ), R. Marsh et al.: Incorporation of the C-GOLDSTEIN efficient climate model into the GENIE framework 961 the freezing point for seawater as a function of surface salinity: where S o is in practical salinity units. Also, in addition to the climatological albedo everywhere, we further apply a surface albedo over sea ice, α ice , allowed to vary as a function of air temperature T a , following Holland et al. (1993), and expanding on Eq. (7) in EM05: The surface temperature of the sea ice, T i , is also diagnosed for use in bulk flux formulae as the temperature "seen" by the atmosphere over sea ice. As outlined by Weaver et al. (2001), T i is found as an iterative solution. The sea-ice transport scheme has been modified to optionally use the same implicit time-step formulation as is used for the atmospheric transport. This option can be necessary with relatively fine horizontal resolution at high latitudes.
In the sea-ice model, vertical fluxes are in SI units, while horizontal fluxes are non-dimensionalized. In summary, the following physical processes are explicitly represented: -Horizontal transport of sea-ice concentration and thickness; -Surface exchange of heat with the atmosphere, according to the balance of incoming shortwave radiation, sensible and latent heat fluxes, and re-emitted long-wave radiation; -Exchange of heat with the ocean, according to the balance of freezing and melting; -Surface exchange of freshwater between the atmosphere and ocean; precipitation is assumed to pass straight to the ocean, ignoring ice, and sublimed water from the ocean surface is added directly to the atmosphere; -Exchange of freshwater with the ocean, according to the balance of freezing and melting.
Following each sea-ice transport time-step, sea-ice concentration is constrained to the range from 0 to 1, and sea ice with height less than a minimum of 0.01 m is removed, both as described in the appendix of EM05. Constants and scale values specific to the sea-ice component are listed in Table A3.

Coupling strategy
Apart from advective drift of sea ice with ocean surface currents, and given prescribed wind forcing, the three eb go gs model components are coupled through vertical heat and freshwater fluxes. As the details of each flux component are fully outlined in Weaver et al. (2001), we only outline here how these components enter the coupling, recapitulating and enlarging on the C-GOLDSTEIN model description of EM05. Positive fluxes are defined as directed into the corresponding model component.
While the coupling scheme forms part of the EMBM module, it is called directly by the main GENIE programme. We employ asynchronous coupling: one "ocean" and "sea ice" time-step per set number of "atmosphere" time-steps, the ocean, atmosphere and sea ice being forced with heat and freshwater fluxes evaluated by the coupling scheme for each ocean and sea ice time-step. For the simulations peresented herein, the number of atmospheric timesteps (per ocean/seaice timestep) is set to 5, to ensure numerical stability in updating the atmosphere.

Heat fluxes
Following Weaver et al. (2001), heat is exchanged between the three components as follows. The corresponding total heat flux into the atmosphere, Q ta , is computed as: Q SW is the net shortwave radiation at the top of the atmosphere, accounting for planetary albedo (see Fig. 5e). C A is an absorption coefficient, parameterizing heat absorption by water vapour, dust, ozone and clouds. Over the oceans, C A = 0.3, i.e., the ocean absorbs 70 % of Q SW . Over land, C A = 1.0; with the most rudimentary terrestrial representation, land has no heat capacity, and all net shortwave radiation must be absorbed in the overlying atmosphere. Q LH is the latent heat released through condensation and subsequent precipitation. Q LW is net (upward minus downward) re-emitted long-wave radiation (the upward part from the ocean or sea ice). Q SH is the sensible heat flux from either ocean or sea ice. Q PLW is the outgoing planetary long-wave radiation.
The total heat flux from the atmosphere into an ocean gridbox (including sea ice), Q t , is partitioned according to fractional sea ice cover, A i : where Q to is the total heat flux into the ice-free ocean from the atmosphere, and Q ti is the total heat flux into sea ice from the atmosphere:  over open ocean and sea ice, respectively (ρ o is the density of seawater, L v and L s are the latent heat of vaporization and sublimation, and E o and E i are the evaporation and sublimation rates). Note that both Q SW,o and Q SW,i are calculated according to planetary albedo (see Fig. 5e), but Q SW,i is further modified by sea ice albedo (see Eq. 2), while sensible and latent heat fluxes are computed according to property differences between the atmosphere and either ocean or sea ice.
The heat flux from sea ice into the ocean, in the sea ice fraction of an ocean grid-box, Q b , is defined by a corrected form of Eq. (10) of EM05: where T o is the temperature of the uppermost ocean level, z the thickness of the uppermost ocean layer, C po is the specific heat of seawater beneath sea ice, and τ i is a relaxation timescale, held constant at 17.5 days. This timescale Interpolated surface wind speed (m s −1 )

Fig. 2.
Annual-mean wind stress (arrows) and derived surface wind speed (colours) for each configuration.
imposes an upper limit of that order of magnitude for the (ocean) time-step in the coupling scheme. Given Eqs. (5) and (7), the net surface heat flux into the ocean, Q H , is given by: From Eqs.
(3), (5) and (6), the atmosphere, ocean and sea ice are directly coupled through Q SH and Q LW . There is also "indirect" coupling through latent heat fluxes. Evaporative cooling of the ocean is followed by latent heat release in the atmosphere upon condensation of the evaporated water vapour and subsequent precipitation. This coupling may be "non-local": evaporated water vapour may be transported away from an oceanic source, to condense out elsewhere. Equations (7) and (8) describe the thermal influences of sea ice and the ocean on each other.

Freshwater fluxes
The total freshwater flux into the atmosphere, through exchanges with other components, is simply the evaporation Interpolated wind speed (m s −1 ) rate over the open ocean fraction of a grid-box plus the sublimation rate over the sea ice fraction, minus the precipitation rate. The freshwater flux into the ocean is applied as F S , an implied salinity flux: where S 0 is a representative salinity, P is the precipitation rate, R is the runoff rate, and F w,i is the net freshwater flux (into the ocean) associated with freezing and melting of sea ice (see EM05 for further details). Runoff is routed according to prescribed runoff directions for each land grid-box (see Sect. 2.3 and Fig. 1). With no snow scheme, we do not store precipitation on sea ice. Instead, any such precipitation is added to the underlying ocean, as if the sea ice is permeable.
The atmosphere and ocean are locally coupled through E and P , while the ocean is further forced non-locally by the atmosphere through R. Sea ice and the ocean are coupled through F w,i : brine rejection on freezing increases salinity; fresh water released on melting decreases salinity.

Meshes, time-stepping and numerics
In this paper, we present the results obtained with four distinct meshes:  -The default mesh, identified as "3636s16l" (also used in the "biogem" configuration); -Doubling 3636s16l resolution in all spatial dimensions, we obtain a mesh identified as "7272s32l"; -For preferential resolution of latitude (in equal increments of 3 • ), but the same resolution of longitude as 3636s16l, is a mesh identified as "3660l16l"; -For near-equal increments of latitude and longitude (equal increments in longitude of 5.625 • and similar near-equal increments of latitude), is a mesh identified as "6432i16l", exactly matching the surface grid resolution of the dynamical atmospheric component available in GENIE (Lenton et al., 2007).
In referring to different meshes, the convention used here is to specify the number of grid cells for longitude first, followed by that for latitude, then appending "s" where we specify constant increments in sine of latitude, "l" where we specify constant increments in latitude, and "i" in a special case where we match the mesh of the IGCM dynamical atmospheric model. The spacing in the longitudinal direction is always in equal increments of longitude, which results in equal-area partitioning of the mesh for the "s" variant. Finally, we specify the number of ocean depth levels followed by an "l". The depth levels are exactly logarithmicallyspaced, with a maximum depth of 5000 m, and we now routinely specify a lowest resolution of 16 depth levels. This default number of levels is doubled over earlier implementations (e.g., EM05, Lenton et al., 2006, Ridgwell et al., 2007a, as the benefits in terms of resolved physics and biogeochemistry are considered to out-weigh the extra computational cost (roughly a doubling of computation time over corresponding 8-level configurations). Depths and thicknesses associated with ocean layers in the 16-level meshes are listed in Table A4. The motivation for these four meshes is according to scientific focus. The 3636s-type meshes afford great efficiency and have been most widely used, for very large ensembles (e.g., , multiple biogeochemical simulations (e.g., Cameron et al., 2005, Ridgwell et al., 2007a, and very long simulations, of the order of 1 Ma, to investigate rock weathering and the carbon cycle (Colbourn, 2011); notably the 16-layer version has now been established as a standard model setup (e.g., Cao et al., 2009). The 7272s32l mesh tests the effect of doubled spatial resolution in isolation, by exactly replicating the land mask and bottom depths of the 3636s16l configuration. The 3660l16l mesh provides notable enhancement of the high-latitude meridional resolution compared to the other configurations, intended to improve representation of water mass formation, deep sinking and cryospheric processes. The 6432i16l mesh exactly matches the surface exchange mesh of the IGCM, eliminat-ing the need for the interpolation that can degrade simulations, and has higher zonal and polar resolution than the 3636s16l mesh but lower meridional resolution in low latitudes (Lenton et al., 2007). GENIE bathymetries are generated by first averaging 5-min (3636s16l and 7272s32l) or 2-min (3660l16l and 6432i16l) gridded elevation data, originated as ETOPO5 (National Geophysical Data Center, 1988) and ETOPO2 (US Department of Commerce, National Oceanic and Atmospheric Administration, 2006) respectively, in model gridcells. Manual intervention is currently necessary to handle narrow land bridges (e.g., Central America) and straits (e.g., Bering Strait, Indonesian Archipelago) that are typically "averaged out" by automatic procedures. The land-sea mask of the 3636s16l, "biogem", and 7272s32l setups are essentially identical to that used by EM05. The bathymetry of the 3636s16l setup corresponds to the bathymetry of "biogem" Annual-mean air temperature ( • C) (a) and specific humidity (g kg −1 ) (b) for 3636s16l and "biogem" configurations; "3636s16l" minus "biogem" differences for annually-averaged air temperature (c) and specific humidity (d).
(which is identical to the bathymetry of the "GENIE-16" setup described by Cao et al., 2009), except for one former land grid cell which now represents the location of the ITF. The bathymetry of the corresponding grid cell was set to a depth (728.8 m) approximating that of the Makassar Strait. The land/sea masks for 6432i16l and 3660l16l are from previous 8-layer model setups (the former outlined in Lenton et al. (2007)), again with the exception of modified grid cells to represent the Makassar Strait, and (in the former) the Davis Strait. The 16-layer bathymetries for the 3660l16l and 6432i16l setups were generated by averaging depths from the ETOPO2 bathymetry onto the model grid first, then iteratively reducing/elevating averaged ocean depths in small decrements/increments in order to reduce large absolute values of the horizontal second derivative below a value of 2 × 10 −8 m −2 . Further shallow depths at and around the Drake Passage were lowered to a minimal depth of 557.8 m. Figure 1 shows ocean topography discretized on the four different meshes. Also required, for a given horizontal mesh, is a runoff direction for land grid-cells (north, south, east, west), and these directions are also indicated in Fig. 1. Topography and catchment area representations are combined in single mesh-specific files. The solution of the barotropic streamfunction requires two further files, identifying disconnected land masses and specifying closed paths in model coordinates around these land masses. One further file specifies grid-cells representative of the Atlantic, the Pacific and the Southern Ocean (except in the "biogem" configuration, where the model auto-detects the basin mask). The protocol for these filenames, and the respective original data sources, are provided in Table 1.
In principle, the model code would allow arbitrary variation of resolution with latitude. However, horizontal resolutions higher than 7272s32l are liable to lead to a degradation of efficiency owing to numerical boundary instabilities that require an implicit scheme for solution (Edwards and Shepherd, 2001). At sufficiently high resolution, an instability related to baroclinic eddy formation is likely to occur, although the model does not include sufficient physics to represent eddies correctly. In practice, the atmosphere and sea ice share the same horizontal resolution as the ocean. Although the GENIE framework provides a facility for the coupling of modules with different horizontal resolution, in the case of the eb go gs configuration, the use of non-matching meshes across the model components would lead to unnecessary interpolation with associated errors and degradation in performance.
In eb go gs, we use a mixture of explicit and implicit numerics, with the latter chosen for atmospheric and sea ice transport (except in the "biogem" configuration, where the sea-ice time-step uses the explicit option). According to the horizontal mesh in use, time-steps are specified as follows:  Edwards and Shepherd, 2001). The convection depth level is a grid-specific diagnostic model output averaged over the whole spin-up simulation.
-number of atmosphere time-steps: asynchronous with default of 5 atmospheric time-steps per ocean time-step (i.e., 480 yr −1 or 960 yr −1 , depending on the horizontal resolution).
The default 3636s16l time-step of 96 yr −1 is chosen to conveniently resolve each month with 8 time-steps.

Prescribed inputs and key parameters
Input/output fields to/from the three components of eb go gs are listed in Table 2. Present values for key parameters of the three components are listed in Table 3. As EM05 originally defined these parameters in the context of four governing equations, we provide only a brief summary here. Advective and diffusive transports of heat (moisture) in the EMBM are parameterized by six (three) key parameters. Ocean mixing in GOLDSTEIN is parameterized by isopycnal and diapycnal diffusivities. Two further key parameters control the ocean circulation: momentum drag is proportional to an inverse timescale, and EM05 further introduced "a constant scaling factor, W ... which multiplies the observed wind stresses in order to obtain stronger and more realistic winddriven gyres". Horizontal diffusion of sea ice thickness and concentration is in proportion to sea ice diffusivity. The choice of parameter values in Table 3 has been guided by experimental experience. Optimally-tuned parameter values for similar model setups of GENIE or closely-related models have previously been obtained by Hargreaves et al. (2004), Beltran et al. (2006), ), and Matsumoto et al. (2008, using a variety of methods. The  study includes the results of the two preceding tuning exercises that used different methodologies, whereas the latter two tunings follow the methodology of . Note that resulting parameter values from these studies apply to the default horizontal mesh (3636s16l or 3636s08l setups), and are not likely to be optimal for other resolutions. One specific tuned parameter set, the "GENIE-16" setup of Cao et al. (2009), is presently favoured when coupling eb go gs with the biogeochemical module of GE-NIE (A. Ridgwell, personel communication, 2011), herein referred to as the "biogem" configuration (Table 3, where we identify the unmodified values in parentheses). This configuration is based on a parameter set originally tuned before the model was subsequently modified by extending the cross section of the Drake Passage, through local modifications of the bathymetry, and by reducing the meridional diffusivity for heat over and around Antarctica as described in the Appendix A of Cao et al. (2009) (A. Ridgwell, personel communication, 2011, and outlined here in Sect. 2.1.2. The value of −56 • N (note sign convention) for the northward extent of the heat diffusion adjustment has been introduced to replace the mesh-specific setting and, for the "3636s16l" and "biogem" cases, both settings have identical effects.
Three parameters from the "biogem" set -the reduction of meridional heat diffusivity at southern high latitudes, the meridional slope for atmospheric heat diffusivity and hori-zontal sea ice diffusivity -have been modified (in the case of sea ice diffusivity the implicit numerical scheme option was also additionally enabled). The corresponding parameter set (see Table 3), subsequently referred to as the "default" parameter set, ensures numerical stability across a wider range of meshes. This parameter set is compatible, and used, with all four mesh options presented in this study, to facilitate comparison of results. The reduction of meridional heat diffusivity at southern high latitudes (relative to "biogem") further acts to reduce errors in ocean and atmospheric temperature (see "Results" below).
In addition to these key parameters, the EMBM also requires several prescribed input fields (see Table 1). In the "biogem" configuration, these fields comprise: (i) Wind stress components (Fig. 2c), interpolated onto the model grid from the wind stress climatology of Josey et al. (1998), read in by and provided via the EMBM module to the ocean component; (ii) Zonal and meridional wind speed fields (Fig. 3c)     wind speed components in the two highest-latitude grid rows are zonally-averaged (10-m winds are used, although they represent vertically-averaged zonal advection of heat, as well as zonal and meridional vertically-averaged zonal and meridional advection of moisture).
For new configurations, the reading of wind-related forcing fields from pre-interpolated and mesh-labelled files has been replaced with a run-time interpolation scheme contained in "genie-wind". This approach facilitates the flexible configuration of alternative meshes in GENIE-1, as no mesh-specific preparation of the wind-related forcing files is required. The surface winds and wind stresses here are longterm averages of 1000 mb wind fields and momentum flux fields, respectively. We use monthly means for 1979-2009, from the NCEP/DOE 2 Reanalysis dataset (Kanamitsu et al., 2002), from the NCEP/DOE 2 Reanalysis data provided by the NOAA/OAR/ESRL PSD, Boulder, Colorado USA, from their website http://www.esrl.noaa.gov/psd/. Interpolations of these fields onto the model meshes are shown in Figs. 2a, c-e and 3a, c-e. As the reanalysis wind stress fields are glob-ally defined, this facility avoids inconsistencies in the windstress forcing fields, previously introduced through limited spatial coverage of the observation-based wind stress climatology of Josey et al. (1998). This problem typically occurred close to continental boundaries, where the model land mask deviates from the location of continents as assumed for the parent dataset. While the new approach produces defined values for the whole model ocean domain, it misinterprets some orographically-induced momentum fluxes as wind-stress forcing of the ocean (particularly in the 7272s32l setup, for which the land-sea mask is coarsely resolved in relation to the model mesh). Also, the previous forcing fields featured no wind stress for some months or the full year over sea ice-covered areas, resulting in zero or reduced windstress on the combined ocean and sea ice system in such regions, which may or may not coincide with sea ice cover in the model. While sea ice in eb go gs is not directly forced by wind stress, it is advected with ocean surface currents that are thus forced. With this model update in mind, we undertake here a comparison of the new and legacy wind forcing (Fig. 4). Amplitudes of wind stress are typically bigger in the reanalysisbased wind stress forcing, notably in the Southern Ocean (Fig. 4a). Differences in the surface winds (Fig. 4b) reflect the different structure of the 10 m and 1000 mb fields (in reanalysis products). Figure 2 also shows the surface wind speed, derived by the model from the wind-stress climatology, which is used by the coupling scheme to parameterise air-sea (and air-sea ice) exchange of heat and moisture. Note that the model zonally averages the surface winds of the two highest-latitude grid rows, which represents a much broader area in the "3636s"type and "7272s"-type meshes compared to the other mesh types. Differences of the derived surface wind stress fields between 3636s16l and "biogem" (Fig. 4c) are prominent near seasonally or perennially sea-ice covered areas in the Southern Ocean and near continental boundaries where the wind stress field is influenced by orography and/or was set to zero in the forcing fields used by "biogem". Differences between Ekman pumping diagnosed for the two different wind stress climatologies (Fig. 4d) include enhanced upwelling along the western continental boundaries of northern South America and Africa, but a weaker upwelling off the Brazilian coast. Figure 5 shows meridional profiles of the heat and moisture diffusivities for the different experiments ( Fig. 5a-d), and the prescribed climatological albedo (Fig. 5e). The meridional profile of heat diffusivity in the atmosphere is controlled by five tunable parameters (the first five 'eb' parameters in Table 3). Note the sharp inflexion in the meridional component of heat diffusivity, from high to low values moving southwards, at around 56 • S. This "feature" of our diffusivity profile enforces the thermal isolation of high southern latitudes, designed to improve the realism of simulated Antarctic climate, Southern Ocean sea ice, and the formation of Antarctic Bottom Water (see Sect. 3.2).

Model evaluation against observations
For basic model evaluation, selected annual-mean observations (see Table 1) are interpolated onto each mesh (using simple bi-linear interpolation), and to each depth level in the ocean. Points of the model grid not fully enclosed by eight points within the ocean domain of the data-based products have been extrapolated by horizontally searching for the nearest oceanic grid point in the data-based field (vertically interpolated to the respective model depth layer). Atmospheric observations of 2 m and 1000 mb air temperature, 2 m specific humidity, and 1000 mb relative humidity are from the NCEP/DOE 2 Reanalysis data provided by the NOAA/OAR/ESRL PSD, Boulder, Colorado, USA, from their website http://www.esrl.noaa.gov/psd/. As for the wind data, these datasets were long-term averaged over 1979-2009. Ocean observations of in-situ temperature and salinity are from the Locarnini et al. (2006) and Antonov et al. (2006) climatologies. Model-observation differences, considered as errors, are evaluated for various fields of the atmosphere and the ocean components as root mean square (RMS) errors, with weighting to account for both the number of model points in the respective field and for differences in the observed variance of the data-based fields (analogous to the error measure used by EM05), in a variant also weighting the individual model points according to variable grid-cell areas or volumes. Table 5 lists these weighted model-observation RMS errors, for the atmosphere and ocean, and for each mesh, where the model data are annual-mean values for the last year of each experiment. The simulated humidity fields used in the RMS error computation are adjusted by removing the diagnosed amount of precipitable water (if any) at the time of recording model output. We return to these errors when discussing global evaluation (Sect. 3.6).

Results
We show results from spin-up experiments with the five configurations. Except for spatial resolution and the time-step length, experiments on the four different meshes are configured with an identical parameter set, which is somewhat modified with respect to the "biogem" configuration (see below). We focus on results from 3636s16l, illustrate and explain differences with respect to "biogem", and selectively consider differences attributed to the choice of mesh.
In all experiments, eb go gs is integrated from an initially uniform state (ocean temperature of 0 • C and salinity of 34.9 psu, and dry atmosphere at 0 • C) for 5000 simulated 976 R. Marsh et al.  years to approximate a near steady, seasonally-oscillating state (see Fig. 6). After a few hundred years, the radiation balance and the atmospheric hydrological cycle are largely equilibrated, evident in the rapid equilibration of globalmean air temperature (Fig. 6c) and specific humidity (not shown). During the model spin-up, the ocean heat content approaches equilibrium, as shown by the evolution of deep ocean temperature and salinity (not shown). Except for the "biogem" setup, after a few millennia, the global ocean circulation and seasonal sea ice distributions are very close to equilibrium, as illustrated by the long-term trends in maxima and minima of Atlantic and Pacific overturning rates (Figs. 6a, b) and sea ice volume (Fig. 6d) . In the "biogem" setup, somewhat stronger drifts occur even after 5000 yr (particularly in the Pacific overturning streamfunction), which are slowly reduced over subsequent millennia of continued simulation (not shown). Broadly speaking, timescales of equilibration are similar across the four meshes, with no systematic differences that can be ascribed to grid resolution. The seasonality of the insolation forcing is evident in the simulation, although subdued in the ocean circulation due to the use of annual-mean wind-stress forcing.
We henceforth consider annual means in the oscillatory near-equilibrium state, as realised after each 5000-yr spin-up, of the atmosphere, ocean, and sea ice (Sects. 3.1-3.3  ocean circulation (Sect. 3.5), and evaluation of global modelobservation differences (Sect. 3.6). Figure 7 shows the annual-mean air temperature (Fig. 7a) and specific humidity (Fig. 7b) for 3636s16l. The air temperature fields are generally smoother than the specific humidity fields. As noted by EM05, the rapid relaxation of specific humidity (total removal of above-threshold relative humidity via precipitation at every time-step of the coupling scheme) means that this field largely reflects the underlying surface temperature, while transport and source terms modify the pattern in a way consistent with unrepresented physical processes. In this experiment, atmospheric diffusivity parameters promote thermal isolation of high southern latitudes to a lesser extent than in "biogem". The resulting differences in air temperature and specific humidity, "3636s16l" minus "biogem" are shown in Fig. 7c,d. With enhanced thermal isolation in the "biogem" setup, Antarctic air temperature is reduced by up to around 15 • C. while Arctic air temperatures are correspondingly higher by around 5 • C. As specific humidity varies with air temperature, it is correspondingly lower in "biogem" by up to 2 g kg −1 around Antarctica, and higher by up to 0.5 g kg −1 in the Arctic. Larger shifts of around ± 2 g kg −1 occur in the tropics, possibly due to differences in the spatial pattern of the applied advective transport fields (see Fig. 4b). For further evaluation of each experiment, Fig. 8 shows annual-mean "model minus observations" differences in surface air temperature (compared to the 1000 mb air temperature from the reanalysis data -see Sect. 2.5). The atmosphere is mostly colder than observed. Some of the negative differences over the continents, in particular over western North America, may be due to the neglect in eb go gs of the orographic effects on air mass trajectories, associated condensation and latent heating. The positive differences around eastern Eurasia and eastern North America may be associated with unrealistic diffusive 978 R. Marsh et al.  transport into the continental interior of heat released over the warm western boundary currents. The largest negative differences in "biogem" and 3636s16l are a too-cold Antarctica in "biogem" and a too-cold Arctic in 3636s16l, for reasons discussed above. Between the four meshes for which parameter settings are fixed, differences are again relatively invariant. Figure 9 shows the annual-mean sea surface temperature (Fig. 9a), salinity (Fig. 9b) and density (Fig. 9c), along with a diagnostic of convection averaged over the full spin-up duration (Fig. 9d). The latter illustrates localized centres of deep convection in the North Atlantic and notably around Antarctica. '3636s16l minus "biogem" differences of these properties and diagnostics are shown in Fig. 10. South of 55-60 • S, surface temperatures are mostly lower in "biogem" by up to 1 • C (Fig. 10a), although in many locations in this region, sea ice formation restricts minimum temperature to the freezing point for seawater, and differences are consequently small. North of this region, with the exception of a few locations, temperatures are higher, by up to 3 • C, in "biogem", consistent with reduced southward heat transport. In 3636s16l, temperature at coastal grid-cells is generally reduced by around 1-3 • C, most notably in upwelling regions and this is likely due to the revised wind forcing (see Fig. 4c, d). The largest positive surface salinity differences are located in the Arctic, associated with differences in sea ice and associated freeze/melt fluxes (see Sects. 3.3 and 3.4). The largest negative salinity differences are caused by the opening of the ITF in the 3636s16l setup (Fig. 10b). The reduced salinity downstream the ITF corresponds to a negative anomaly in density in the uppermost ocean layer (Fig. 10c). Despite subtle differences in surface density in the Southern Ocean, there is considerably more convective mixing in 3636s16l (Fig. 10d). Convection associated with dense water formation is locally strongest immediately adjacent to Antarctica, on average spanning about 8 depth levels (which,  Fresh water flux into the ocean (m yr −1 ) Fig. 18. Annual-mean components of freshwater flux (m year −1 ) into the ocean for 3636s16l (maps) and "3636s16l" minus "biogem" differences (contours). Bold-outlined areas in (c) represent the three pairs of regions for Atlantic-to-Pacific freshwater flux adjustment as shown by .

Ocean properties
assuming convection always occurs at the surface ocean, extends to a depth of 1158 m).
For the Atlantic and Pacific basins, we also show temperature and salinity in meridional-vertical sections for 3636s16l in Fig. 11. Several features are clear from these sections. Temperature and salinity values below around 1000 m are influenced by the inter-hemispheric overturning circulation, in particular the export of North Atlantic Deep Water (NADW). At the low horizontal and vertical resolutions used here, representation of Antarctic Bottom Water (AABW) is aided by the thermal isolation of Antarctica (see above). '3636s16l minus "biogem" differences are shown in Fig. 12. In "biogem", the abyssal ocean is considerably warmer, most notably the Atlantic sector where deep-water temperature is higher by up to 2 • C (Fig. 12a). Deep Atlantic salinity is correspondingly higher in 'biogem" by up to around 0.2 psu. The implication is that warmer and more saline NADW forms, circulates and mixes in "biogem". For sub-surface water masses in the upper 1000 m of the Atlantic, south (north) of the Equator, temperatures and salinities are lower (higher) in "biogem". The dominant differences in the upper 1000 m of the Southern Hemisphere are indicative of the formation and northward spreading of cooler and fresher intermediate water in "biogem". In the North Pacific, negative temperature and salinity differences at around 1000 m (Figs. 12c, d)  For model-observations evaluation, of each experiment, Figs. 13 and 14 show annual-mean "model minus observation" differences for zonally-averaged ocean temperature in the Atlantic and Pacific, respectively, including their extensions into the respective sectors of the Southern Ocean. Errors are negative throughout much of the ocean below around 1000 m depth (with the exception of "biogem", for which warmer NADW is a notable influence), and for the Atlantic sector of the Arctic Ocean in all but the 3660l16l mesh. These anomalies are associated with cold biases in surface climate (see below). In all five experiments, the upper layers in northern high latitudes of the Atlantic are dominated by negative errors reaching in excess of −5 • C, consistent with limited northward extent of the Atlantic overturning (see Fig. 21 below) and weak associated heat transport.

Sea ice
Sea ice is particularly sensitive to parameters and the different high-latitude resolutions provided by the various mesh types (see Fig. 6d). For this reason, in Fig. 15 we show annual-mean sea ice height and concentration, for all five experiments. The extent of seasonality (in months of ice-cover) can be approximated by multiplying annual-mean concentration by 12. In all experiments, thickest and most perennial sea ice develops in the Arctic, in reasonable agreement with observations (e.g., National Snow and Ice Data Center (2010)). Arctic sea ice is substantially thicker in 3636s16l than in "biogem", while with the higher resolution in the high latitudes of the alternative meshes, the model produces even thicker sea-ice cover, reaching up to 3 m (in 3660l16l). In the Southern Ocean, more extensive sea ice develops in "biogem", compared to 3636s16l, with perennial sea ice adjacent to Antarctica (albeit very little) only in the former experiment (see Fig. 6d). This is attributed to the larger reduction of atmospheric heat diffusivity in the former experiment. Sea ice in both hemispheres is broadly similar in 3636s16l and 7272s32l (Fig. 15a, c). More realistic sea ice distribution (in particular around Antarctica) is obtained with the higher polar resolution in 6432i16l and 3660l16l (Fig. 15d, e), by comparison with observations. The reasons for differences between meshes may not be directly related to spatial resolution, but may rather be due to feedbacks involving the atmosphere and the ocean.

Vertical and internal fluxes of heat and moisture
For 3636s16l only, Figs. 16-18 show annual-mean vertical and internal components of heat and freshwater flux, for both atmosphere and ocean. While these year-5000 annual means may be subject to slight inter-annual variability, they are still representative for typical flux fields of the near-equilibrium state at the end of the simulations. The short-wave heat flux into the atmosphere (Fig. 16a) declines polewards due to the angle of incidence of insolation (accounting for orbital parameters), and is further reduced with latitude by planetary albedo (see Fig. 5e). Further demarkation at land-sea boundaries, with stronger fluxes over land, is a consequence of zero heat capacity of land in the absence of a more sophisticated land scheme (the heat capacity of land is neglected). Consequently, the atmosphere absorbs all incident shortwave radiation over land, while the atmosphere absorbs 30 % of short-wave radiation over the ocean. The atmosphere loses heat through net long-wave heat fluxes (Fig. 16b), dominated by the outgoing planetary longwave radiation. The sensible heat flux into the atmosphere (Fig. 16c) reaches high values where the surface ocean loses heat as relatively warm water advects towards deep sinking sites, notably the North Atlantic subpolar gyre, and in a zonal strip around the high-latitude Southern Ocean. Such heat exchange is strongest where the annual sea-ice coverage is lowest (see Fig. 15), and the ocean is in direct contact with a cold atmosphere. The latent heat flux into the atmosphere (Fig. 16d) is directly proportional to the precipitation rate through the associated condensation of water vapour (see Fig. 17b below). The highest fluxes are located in the tropics, but secondary maxima coincide with strong evaporation in the vicinity of warm ocean currents. Combined together, these four components yield a familiar climatological pattern of net heat flux into the atmosphere (Fig. 16e), with heat gain (loss) in low (high) latitudes that are balanced by meridional heat transports (not shown). Notable differences between 3636s16l and "biogem" are more extreme maxima in latent heat flux around the Equator and a shift of sensible heat flux closer to Antarctica due to reduced sea-ice cover in 3636s16l. Table 1. Mesh-specific files required for the eb go gs experiments shown: <world> is "wv3jh2" for 3636s16l, "worjh2" for "biogem", "w32jh2" for 7272s32l, "igcm16" for 6432i16l, and "3660l2" for 3660l16l; filenames with the prefix "+" are read in by all configurations but are only actively used when the run-time interpolation of wind forcing fields is not enabled (herein only by configuration "biogem"); the prefix "*" denotes files not required for "biogem"; the prefix "%" denotes alternative fields used for the error estimation.  (Josey et al., 1998) "eb" +taux v.interp Ocean wind-stress component SOC climatology (Josey et al., 1998) "eb" +tauy u.interp Ocean wind-stress component SOC climatology (Josey et al., 1998) "eb" +tauy v.interp Ocean wind-stress component SOC climatology (Josey et al., 1998) "eb" * NCEP-DOE Reanalysis 2 ltaa 1000 mb wind speed.nc Zonal and meridional 1000 mb wind speed fields NCEP-DOE Reanalysis 2 "eb" * NCEP-DOE Reanalysis 2 ltaa wind stress.nc Zonal and meridional momentum flux into the atmosphere (negative wind stress) NCEP/DOE Reanalysis 2 (Kanamitsu et al., 2002) "eb" Model-data error computation %NCEP-DOE Reanalysis 2 ltaa 1000 mb T r.nc Observational 1000 mb air temperature and relative humidity fields NCEP/DOE Reanalysis 2 (Kanamitsu et al., 2002) "eb"

%NCEP-DOE Reanalysis 2 ltaa surface T q.nc
Observational 10 m air temperature and specific humidity fields NCEP/DOE Reanalysis 2 (Kanamitsu et al., 2002) "eb" WOA05 an TS.nc Observational ocean temperature (in situ) and salinity fields World Ocean Atlas 2005 Antonov et al., 2006) "go" The short-wave heat flux into the ocean (Fig. 17a) is a constant 70 % of net shortwave radiation at the top of the atmosphere (30 % being absorbed in the overlying atmosphere) over most of the ocean surface. Varying with latitude as the atmospheric counterpart, short-wave heat flux into the ocean is further reduced at high latitudes due to the higher albedo associated with sea ice. The net long-wave heat flux (Fig. 17b) is small relative to the atmospheric counterpart, representing the balance of re-radiated long-wave heat fluxes by the ocean (upward) and the atmosphere (downward). The sensible heat flux into the ocean (Fig. 17c), negative everywhere, is identical and opposite to the counterpart flux into the atmosphere (over the ice-free ocean). The latent heat flux into the ocean (Fig. 17d), again negative everywhere, mirrors the pattern of evaporation (see Fig. 18a below). Combined together (Fig. 17e), these four components yield a pattern of net heat flux into the ocean at low latitudes, with strong localized heat loss over western boundary currents and at sites of high latitude convection (see Fig. 9d), again balanced by meridional heat transport (not shown). The largest differences between the net surface heat fluxes in 3636s16l and "biogem" are found around Antarctica, associated with reduced sea-ice extent in 3636s16l.
Evaporation is equivalently a negative freshwater flux into the ocean (Fig. 18a), with largest absolute values over the subtropical gyres. Precipitation (Fig. 18b) to an extent coincides with strong evaporation. As is typical of low complexity atmospheric models such as the EMBM, precipitation over land is poorly represented, even with moisture advection by specified winds. Substantial areas of continental interior Table 2. Input/output fields to/from the individual modules. Note, while for eb go gs configurations (as for C-GOLDSTEIN, see EM05), we define neither land temperature nor standing water on land, differences between air and land temperatures, as well as land moisture, are important coupling quantities when other land options are used, and so are made available when eb go gs is coupled to the land surface scheme ("ENTS") of Williamson et al. (2006). Also, when the atmospheric chemistry module ("ATCHEM") is included, further input to the atmosphere include greenhouse gas concentrations. When the ocean module is coupled to the ocean biogeochemistry module ("BIOGEM") a range of fields from the three modules are made available to BIOGEM, including a number of passive oceanic tracers which are subjected to ocean transport and are modified by BIOGEM.

Field
"eb" "go" "gs" "genie-wind"  Williamson et al. (2006), but such enhancement requires coupling with the land scheme -an elaboration on the eb go gs configurations that are documented here. The precipitation that does fall over continents instantaneously reaches the ocean as runoff (shown in combination with the Atlantic-to-Pacific freshwater flux adjustment in Fig. 18c). The consequence of seasonal freezing and melting of sea ice is a net freshwater flux that effectively removes freshwater in regions of sea ice formation and adds freshwater in regions of sea ice melt (Fig. 18d). Combined together (Fig. 18e), these four components yield net freshwater gain at low and middle latitudes, with net freshwater loss in the subtropics. Locally, strongest net freshwater gain is associated with major rivers. In comparision with "biogem", the most prominent differences in freshwater input to the ocean are due to enhanced precipitation in the tropics in 3636s16l and consequent higher runoff at low latitudes. The enhanced precipitation is in turn asso-ciated with warmer low-latitude climate in biogem (compare model-data differences in Fig. 8a, b). Figure 19 shows the barotropic streamfunction for each of the five experiments and Table 4 lists the transport strengths through selected straits between isolated land masses. The horizontal circulation is dominated by the circumpolar flow, which is sensitive to wind and buoyancy forcing, in combination with friction. Barotropic flow around disconnected land masses is computed following the evaluation of closed-path integrals around the islands, individually specified for each of the five configurations (see Fig. 1). While in the analytical limit of the integration, the barotropic values over the islands are independent of the exact location of these paths, truncation errors due to the discretisation of the integral are introduced. A multiplicity of potential paths exists for the discretised model geographies, and so the realised steadystate solution depends on the choice of the path. For the strength of circumpolar flow, we obtain a range of values by specifying the path integral around Antarctica for different permissible completely-straight paths in some of the configurations (Table 4). This model artefact was not previously documented.

Ocean circulation
The extent of friction in the Drake Passage in particular depends on meridional resolution -with fewer grid-cells between Antarctica and South America, frictional effects are more severe due to the proximity of lateral boundaries and (if any) shallow topographic features (see EM05). While there is considerable variation between experiments, horizontal transports through the Drake Passage are in general weak compared to observation-based estimates (Table 4). This is in part a consequence of the non-viscous representation of frictional effects in GOLDSTEIN (see Killworth, 2003).
For the same wind forcing, differences between the meshes are largest in the Southern Ocean, where transport through the Drake Passage ranges widely, from 46 Sv in 6432i16l to 86 Sv in 7272s32l. Following the above reasoning, these differences are due in part to the number of grid-cells, and resolution of shallow topography, across the Drake Passage -with fewer grid-cells (as for 3636s16l and 6432i16l), the frictional influence of adjacent coasts and shallow topographic features is stronger and this weakens transport. The substantially weaker circumpolar transport in 6432i16l, compared to the experiments on different meshes with the same parameter settings, may also be associated with differences in buoyancy forcing over the Southern Ocean. This is further evidenced by the contrast between Drake Passage transport in "biogem" (37 Sv) and in 3636s16l (59 Sv). A major difference between these two experiments is in the deep Southern Ocean, where temperatures in 3636s16l are substantially lower (see Fig. 12a). Table 4. Annual-mean net transport through straits (positive transport is directed northward or westward) simulated by the five configurations and data-based estimates. The last row lists the net transport through the Drake Passage simulated with configurations where the location for the path integral around Antarctica was varied among all permissible straight paths (the value in parentheses indicates the grid row for this path on the respective mesh). With relatively smaller corresponding salinity differences, deep waters adjacent to Antarctica are consequently denser, isopycnals slope more steeply upwards across the Southern Ocean, and a substantially stronger circumpolar current is supported by the lateral pressure gradient. Table 4 lists net volume transport along pathways into and out of, as well as within, the Arctic Ocean, which are facilitated in the grids with higher resolution in this area. Net flow through the Bering Strait is into the Arctic Ocean and ranges from 0.2 to 2.3 Sv across the four mesh variants with identical wind forcing, bracketing the data-based estimate of Roach et al. (1995). In 3636s16l and 7272s32l, this is balanced by simulated outflow of Arctic waters across the Greenland, Iceland, and Norwegian (GIN) Seas (contrary to data-based estimates which, if assumed not to have changed substantially between the respective dates of their measurements, suggest net inflow across this section (Roach et al., 1995;Curry et al., 2011)). The 3660l16l simulation features net outflow across the GIN Sea, partially offset by outflow through the Davis Strait. Net transport through the Davis Strait in the 6432i16l setup exceeds the net inflow through the Bering Strait, and consequently must be balanced by net inflow across the GIN Seas. Enhanced spatial resolution in the region of the GIN Seas seems to facilitate the development of cyclonic circulation in the high-latitude northern Atlantic and its northward extension representing the recirculation of Atlantic waters across the Barents Sea and along the East Greenland Current, which is best expressed in 6432i16l (Fig. 19d). Simulated ITF strength is within, or slightly below, the data-based estimates for the transport along Makassar Strait (Gordon et al., 1999(Gordon et al., , 2008, but systematically below the estimate of Ganachaud and Wunsch (2000) for the total transport across the Indonesian Archipelago.
In Sect. 3.2, we showed that there are large differences in deep and bottom water properties, attributed to differences in deep water formation at high latitudes in the Northern and Southern Hemispheres. As the overturning circulation is in general sensitive to the resolution at high latitudes, in Figs. 20 and 21 we show annual-mean meridional overturning streamfuctions (for the global ocean and the Atlantic, respectively) for all five experiments. Substantial differences are apparent in both hemispheres, notably in Table 5. Weighted and unweighted root mean square errors, for the atmosphere and ocean, for each experiment. Annual-mean observations are gridded onto each mesh (see Table 1 and Sect. 3.6 for the sources of the datasets and further details). The computation of the sums over individual grid points follows the method of EM05, the integrals over the respective model domain additionally weights the individual grid points according to the volume or area of the corresponding grid box or cell. Note that an error value of 1 would be achieved by hypothetical model output of a uniform field with the same average value as the data-based field. the Southern Ocean, although the Atlantic overturning is in "Conveyor" mode for each experiment. This is facilitated by the specification of an explicit additional Atlantic-Pacific moisture flux, effectively a flux correction (default value 0.727 × 0.32 Sv = 0.23 Sv -unmodified from the value used in "biogem"). Without this additional flux, the typical pattern of Atlantic overturning does not develop, in all configurations. In the Pacific, the overturning streamfunction (not shown) is characterized by the northward import and upwelling of deep water. For the five configurations, annuallyaveraged Atlantic overturning reaches 11-15 Sv. A clear difference between experiments is the presence or absence of an AABW cell in the Atlantic sector (anti-clockwise overturning as shown by negative values below 3000 m in Fig. 21): absent in "biogem", 6432i16l and 3660l16l; present in 3636s16l and 7272s32l. This abyssal cell is sensitive to differences in deep-water density, in turn dependent on processes at North Atlantic and Southern Ocean deep water formation sites.

Global evaluation
The global errors listed in Table 5 indicate relatively little variation, across the four experiments using the "default" parameter set, in the different measures of RMS error, indicating that global-mean quantities are relatively insensitive to the choice of mesh, given the identical parameter settings.
The most notable differences are a better simulation of relative humidity in the atmosphere and a worse sea-surface salinity distribution in 3660l16l and 6432i16l, respectively. Differences in errors between the "biogem" setup and the simulations with the "default" parameter set are slightly more substantial, largest difference being found for the air temperature field. When the dimensions of the grid cell areas (or volumes) representing the model points are taken into account, the different emphasis on spatial resolution for the different meshes can be seen in the RMS errors. Sea-surface temperature, for example, is closer to observations with enhanced resolution in the high latitudes, while sea-surface salinity is more realistic in the simulations with enhanced latitudinal resolution in the lower latitudes. These errors are presented here to illustrate how different versions of GENIE may be systematically co-evaluated with a new model option. However, since no attempt has been made here to tune key parameters for each of the new meshes, the corresponding simulations cannot be adequately evaluated on the basis of the errors listed in Table 5.

Summary and discussion
We have described the basic climate model, "eg go gs", that is facilitated by the GENIE framework. While some of the model description is retrospective and can be found  Edwards et al., 1998;Edwards and Shepherd, 2001  amenable for investigations of Earth System dynamics at selected time-slices in the geological past. Future papers will report a series of updates on the present results, as progress continues in several areas: parameter tuning; methodologies for palaeoclimate experiments; further developments of GENIE-1 components that enhance complexity and/or add new mechanisms. In summary, future developments of GENIE-1 are expected to lead to improvements on the results reported here, and in general to advance the model's applicability for future climate change and paleoclimate studies.

Appendix A
Variables in the ocean, atmosphere and sea ice components of eb go gs are non-dimensionalized according to dimensional scales. Tables A1-A3 list the respective values for dimensional scalings, along with constants, offsets, and selected parameters, for the three components. Table A4 lists depths and thicknesses associated with ocean layers in the "biogem", 3636s16l, 3660l16l and 6432i16l meshes.