Description of the Earth system model of intermediate complexity LOVECLIM version 1.2

The main characteristics of the new version 1.2 of the three-dimensional Earth system model of intermediate complexity LOVECLIM are briefly described. LOVECLIM 1.2 includes representations of the atmosphere, the ocean and sea ice, the land surface (including vegetation), the ice sheets, the icebergs and the carbon cycle. The atmospheric component is ECBilt2, a T21, 3-level quasigeostrophic model. The ocean component is CLIO3, which consists of an ocean general circulation model coupled to a comprehensive thermodynamic-dynamic sea-ice model. Its horizontal resolution is of 3 ◦ by 3, and there are 20 levels in the ocean. ECBilt-CLIO is coupled to VECODE, a vegetation model that simulates the dynamics of two main terrestrial plant functional types, trees and grasses, as well as desert. VECODE also simulates the evolution of the carbon cycle over land while the ocean carbon cycle is represented by LOCH, a comprehensive model that takes into acCorrespondence to: H. Goosse (hugues.goosse@uclouvain.be) count both the solubility and biological pumps. The ice sheet component AGISM is made up of a three-dimensional thermomechanical model of the ice sheet flow, a visco-elastic bedrock model and a model of the mass balance at the iceatmosphere and ice-ocean interfaces. For both the Greenland and Antarctic ice sheets, calculations are made on a 10 km by 10 km resolution grid with 31 sigma levels. LOVECLIM1.2 reproduces well the major characteristics of the observed climate both for present-day conditions and for key past periods such as the last millennium, the mid-Holocene and the Last Glacial Maximum. However, despite some improvements compared to earlier versions, some biases are still present in the model. The most serious ones are mainly located at low latitudes with an overestimation of the temperature there, a too symmetric distribution of precipitation between the two hemispheres, and an overestimation of precipitation and vegetation cover in the subtropics. In addition, the atmospheric circulation is too weak. The model also tends to underestimate the surface temperature changes (mainly at low latitudes) and to overestimate the ocean heat uptake observed over the last decades. Published by Copernicus Publications on behalf of the European Geosciences Union. 604 H. Goosse et al.: Description of LOVECLIM version 1.2 74

1 Introduction LOVECLIM (Fig. 1) is a three-dimensional Earth system model of intermediate complexity (EMIC, Claussen et al., 2002), i.e. its spatial resolution is coarser than that of stateof-the-art climate General Circulation Models (GCMs) and its representation of physical processes is simpler.In LOVE-CLIM, the most important simplifications are applied in the atmospheric component because it is usually the most demanding one in terms of computing time in GCMs.Thanks to those modelling choices, LOVECLIM is much faster than GCMs.On one single Xeon processor (2.5 Ghz), it is possible to run 100 years, with all the components activated, in about 4 h of CPU time.This is a key advantage as it is affordable to perform large ensembles of simulations (as required to test the influence of parameter choices or to analyse natural variability of the system) and the long simulations needed to study past climates and long-term future climate changes.Compared to some other EMICs, LOVECLIM includes a 3-D representation of the system, facilitating the description of some physical processes such as the formation and development of weather systems as well as the comparison with data coming from different regions.
The first two components of LOVECLIM, which were coupled at the end of the 1990's, are the atmospheric model ECBilt (Opsteegh et al., 1998) and the sea-ice-ocean model CLIO (Goosse and Fichefet, 1999), forming what has been later referred to as ECBilt-CLIO2 (e.g., Goosse et al., 2001Goosse et al., , 2002)).Those two components are still presently the core of LOVECLIM, but with significant improvements compared to the original versions.In particular, the radiative scheme and the parameterization of the surface fluxes in ECBilt have been completely revised (e.g., Schaeffer et al., 1998Schaeffer et al., , 2004, see http://www.knmi.nl/onderzk/CKO/differences.html).Initially, in ECBilt-CLIO2, ECBilt and CLIO were interacting through the OASIS software (Terray et al., 1998).This has been modified in later versions where new Fortran routines, specifically developed for the model, take care of the exchanges between all the model components.
ECBilt-CLIO was further coupled to the terrestrial biosphere model VECODE (Brovkin et al., 2002), leading to ECBilt-CLIO-VECODE (e.g., Renssen et al., 2003Renssen et al., , 2005)).More recently, two additional components were added (Driesschaert et al., 2007): the ocean carbon cycle model LOCH (Mouchet and Franc ¸ois, 1996) and the ice sheet model AGISM (Huybrechts, 2002).As the list of acronyms ECBilt-CLIO-VECODE-LOCH-AGISM was becoming too long, it has been decided to form a new acronym, based on the names of all model components: LOVECLIM which stands for LOch-Vecode-Ecbilt-CLio-agIsM.For simplicity, the new name LOVECLIM should be used even if only some components of the model are activated in a particular study.
ECBilt-CLIO and LOVECLIM 1.0 have been publicly released on the KNMI (Koninklijk Nederlands Meteorologisch Instituut) webiste (http://www.knmi.nl/onderzk/CKO/ecbilt. html) and UCL (Université catholique de Louvain) website (http://www.climate.be/modx/index.php?id=81), respectively.However, the public version of LOVECLIM does not include LOCH and AGISM, as the main developers of those two components wish that potential users contact them first to organize a collaboration before obtaining the permission to activate those parts of the code.
In contrast to LOVECLIM 1.0, version 1.1 of LOVE-CLIM (Goosse et al., 2007) has not been publicly released.However, the new LOVECLIM1.2,which is publicly available since December 2009 (http://www.astr.ucl.ac.be/ index.php?page=LOVECLIM%40Description), is very similar to LOVECLIM 1.1 regarding the physics of the model.Some minor modifications were included and some small bugs, which had limited impacts on model results, have been corrected (http://www.astr.ucl.ac.be/index.php?page= LOVECLIM@bugs).In addition, some technical updates have been performed before the official release.In particular, a standard set up for simulating the Last Glacial Maximum (LGM) climate is now available (Roche et al., 2007).
However, no full description of the model is currently available.For each new version, only the new components and the major differences compared to previous versions were described.As a consequence, in order to determine exactly which processes are represented in a version of LOVECLIM, a new user or a scientist interested in model results has to follow the full history of the code over the last 10 years.He/she will thus likely miss some elements because they are too briefly mentioned or only available in internal reports.In addition, he/she will not know for sure if some physical parameterizations or model parts described in early papers are still valid for the latest versions.
We take here the opportunity of the release of LOVE-CLIM1.2 to describe in more detail the present state of the model.We will not discuss extensively all the model equations and parameterizations as this would correspond to hundreds of pages.Nevertheless, the main characteristics of the model will be described and a short evaluation of model results performed.We consider that it is sufficient, in the large majority of cases, for new users and to estimate if the model is an adequate tool for performing a particular analysis (as well as to estimate the associated limitations).Scientists interested in a specific point are referred to the cited papers, the present manuscript providing an up-to-date list of useful references and web addresses where the important internal reports can be obtained.

ECBilt: the atmospheric component
The atmospheric model, developed at KNMI, was first coupled to a simple ocean model (which used a flat bottom) and a thermodynamic sea-ice model (e.g., Haarsma et al., 1996;Opsteegh et al., 1998;Selten et al., 1999, Weber andOerlemans, 2003).Those ocean and sea-ice components have been removed and replaced by CLIO, keeping only the atmospheric part in ECBilt-CLIO and in LOVECLIM.
ECBilt has a dynamical core derived from the work of Marshall and Molteni (1993).It is governed by the equation for q, the quasi-geostrophic potential vorticity, written in isobaric coordinate (Holton et al., 2004;Opsteegh et al., 1998): q is defined as V ψ is the rotational component of the horizontal velocity, f is the Coriolis parameter, f 0 is f at 45 • (north and south), k d and k r are diffusion and damping coefficients, R is the ideal gas constant, c p is the specific heat for constant pressure, σ is the static stability parameter, α is the specific volume, Q is the diabatic heating, F ζ contains the ageostrophic terms in the vorticity equation and F T is the advection of the temperature by the ageostrophic wind.Equation ( 1) is written with ψ, the streamfunction, as an independent variable.ψ is thus the main variable in the dynamical core of ECBilt.ψ is related to ζ , the vertical component of the relative vorticity vector, by Knowing ψ, it is then possible to compute the geopotential height φ, using the linear balance equation: The temperature T is computed from φ sing the hydrostatic equilibrium and the ideal gas law: The ageostrophic terms F ζ and F T are included in Eq. ( 1) in order to improve the representation of the circulation at low latitudes, in particular the Hadley cells.These terms are obtained by computing the vertical velocity and the horizontal divergence diagnostically (Opsteegh et al., 1998).
www.geosci-model-dev.net/3/603/2010/Geosci.Model Dev., 3, 603-633, 2010 is the quasi-geostrophic potential vorticity, T is the air temperature, Ts is the surface temperature and qa is the total water content.is the streamfunction, q is the quasi-geostrophic potential vorticity, T is the air temperature, T s is the surface temperature and q a is the total water content.
Equation ( 1) is solved using spectral methods with a horizontal T21 truncation and three vertical levels at 800 hPa, 500 hPa and 200 hPa (Fig. 2).This corresponds in the physical space to a grid resolution of about 5.6 • in latitude and in longitude.The radiative scheme and the thermodynamic exchanges between the layers and with the surface are computed in this physical space.Temperature is obtained at the surface and at the 650 hPa and the 350 hPa horizons.The model also contains a thermodynamic stratosphere.
The humidity in the atmosphere is represented in ECBilt by a single prognostic variable: the total precipitable water content between the surface and 500 hPa.This variable is transported horizontally using a fraction (60%) of the sum of geostrophic and ageostrophic winds at 800 hPa to take into account the fact that humidity is generally higher close to the surface where wind speeds are lower.Above 500 hPa, the atmosphere is assumed to be completely dry in the model.All the water that is transported by atmospheric flow above this 500 hPa level thus precipitates.Precipitation also occurs if the total precipitable water in the layer is above a relevant threshold (in the LOVECLIM1.2,this threshold is set equal to 0.83 times the vertically integrated saturation specific humidity below 500 hPa, assuming a constant relative humidity in the layer, see Table 1).The convection and associated precipitation are parameterized as in Held and Suarez (1978).
The longwave radiative scheme of ECBilt is based on a Green's function method (Chou and Neelin, 1996;Schaeffer et al., 1998).The following formula is applied for all the model levels: where Flw is the longwave flux, Fref is a reference value of the flux when temperature, humidity and the concentrations of greenhouse gases are equal to the reference values, FG is a function allowing one to compute the contribution associated with the anomalies compared to this reference in the vertical profile of temperature (T ) and in the concentrations of the various greenhouse gases in the atmosphere (GHG ).The last term represents the anomaly in the longwave flux due to the anomaly in humidity q a (see Schaeffer et al., 1998 for an explicit discussion of those terms).The coefficients Fref, G1 and those included in the function FG are spatially dependent.amplw and explw are adjustable coefficients to take into account the uncertainties in the model, in particular those related to its crude representation of the changes in the vertical profiles of temperature and humidity.In LOVECLIM1.2,explw is equal to 0.40; amplw is equal to 1, except between 15 • S and 15 • N, where it is equal to 1.8.All the reference states are derived from a climatology based on the NCEP-NCAR reanalysis (Kalnay et al., 1996).Equation ( 6) is applied for both clear sky and overcast conditions.The total upward and downward longwave fluxes are then the weighted average of the two contributions as a function of the cloud cover, using prescribed clouds (ISCCP D2 dataset, see Rossow et al., 1996).
The downward and upward shortwave fluxes in ECBilt are computed at the 3 levels in the atmosphere, at the surface and at the top of the atmosphere using also a linearised scheme.The transmissivity of the atmosphere (as the cloud cover, see above) depends on the location and the season but is not computed prognostically.The surface albedo is a function of the fraction of the grid box covered by ocean, sea ice, trees, desert and grass (see Sects.2.2, 2.3 and 2.7).The insolation at the top of the atmosphere is obtained using the orbital parameters computed following Berger (1978).
The surface fluxes of sensible and latent heat are computed from estimates of temperature, humidity and wind speed at 10 m and from the characteristics of the surface using standard bulk formulae.The wind speed at 10 m is supposed to be equal to 0.8 times the wind speed at 800 hPa.For the temperature and humidity, the extrapolation from the higher levels is based on anomalies compared to spatially dependent reference profiles derived from the NCEP-NCAR reanalysis (Kalnay et al., 1996), as in the longwave radiative scheme.
The land-surface model is part of the ECBilt code and has the same grid as the atmospheric model.The surface temperature and the development of the snow cover are computed by performing the heat budget over a single soil layer, which has a spatially homogenous heat capacity.For the moisture, a simple bucket model is used.The maximum water content of the bucket is a function of the vegetation cover.If, after evaporation, precipitation and snow melting, the water content is higher than this maximum, the water is transported immediately to an ocean grid point corresponding to the mouth of the river whose basin includes the model grid box.
The equations governing the ocean flows are deduced from the Navier-Stokes equations written in a rotating frame of reference with some classical approximations such as the Boussinesq approximation, the thin shell approximation, and the hydrostatic approximation.The effects of small-scale processes, not explicitly represented by the model, are included in the momentum equation using a simple harmonic operator along the horizontal.For the scalar quantities (in particular potential temperature and salinity), the model relies on both the isopycnal mixing formulation (Redi, 1982), using the approximation of small slopes (Cox, 1987), and the eddy-induced advection term, as proposed by Gent and McWilliams (1990) (see also Mathieu and Deleersnijder, 1999 and Table 2).The parameterization of vertical mixing (Goosse et al., 1999) is derived from Mellor and Yamada's level 2.5 model (Mellor and Yamada, 1982).The vertical viscosity and diffusivity are considered to be proportional to the characteristic velocity (q t ) and length (l) of the turbulent motions.The characteristic velocity q t is computed through a prognostic differential equation for the turbulent kinetic energy, while l is derived from a simple diagnostic equation.While applied over the whole water column, this turbulence closure is mainly active in the surface layer.At depth, the vertical viscosity and diffusivity is generally equal to a background value which follows a profile similar to the one proposed by Bryan and Lewis (1979).In addition, a convective adjustment scheme is applied when the water column is statically unstable on a vertical depth range greater than 100 m.This is achieved by increasing the vertical diffusivity to 10 m 2 s −1 .
In order to improve the representation of the dense water that flows out of the continental shelves and descends toward the bottom along the continental slope, CLIO includes Campin and Goosse's (1999) parameterization of downsloping currents.If the density of a grid box on the continental shelf (or on a sill) is higher than the density of the neighbouring box over the deep ocean at the same depth, shelf water flows along the slope until it reaches a depth of equal density.In order to verify volume conservation, this transport is compensated by a vertical and then horizontal return flow from the deep ocean to the shelf.
CLIO has a free surface.To avoid imposing for all the model equations the small time step needed to explicitly resolve fast external inertia-gravity waves, the split-explicit method is applied (Gadd, 1978).The numerical integration is carried out in two stages: the depth-integrated part (or barotropic mode) and the depth-dependent one with a zero vertical mean (baroclinic mode).The low numericalcost 2-D barotropic mode, which includes the surface gravity waves, is integrated with a small time step (5 min), while the more expensive 3-D baroclinic mode is solved using a much longer time step (3 h).
The various variables are staggered on a B-grid following the classification of Arakawa (Mesinger and Arakawa, 1976) (Fig. 3).The horizontal discretization is based on spherical coordinates, using a resolution of 3 • in longitude by 3 • in latitude and a realistic bathymetry (within the limits of the model resolution).Actually, two spherical subgrids (Deleersnijder et al., 1997) are associated to avoid the singularity at the North Pole (Fig. 4).The first one is based on classical longitude-latitude coordinates.It covers the Southern Ocean, the Pacific Ocean, the Indian Ocean and the South Atlantic.The second spherical subgrid has its poles located at the equator, the "north pole" in the Pacific (111 • W) and the "south pole" in the Indian Ocean (69 • E).The remaining parts of the ocean are represented on this "rotated" grid, i.e., the North Atlantic and the Arctic.The two subgrids are connected in the equatorial Atlantic where there is a correspondence between the meridians of the South Atlantic on one  V are the two components of the barotropic velocity, η the surface elevation, u, v, and w the three components of the velocity, S the salinity, θ the potential temperature, q 2 (two times) the turbulent kinetic energy, and Ks and K u the vertical diffusion and vertical viscosity.grid and the parallels of the other grid in the North Atlantic.Because of the grid system, the direct connection between the Pacific and the Arctic through the Bering Strait is not explicitly computed, but the transport there is parameterized by a linear function of the cross-strait sea-level difference in accordance with the geostrophic control theory (Goosse et al., 1997).The vertical discretization follows the simple so-called "z-coordinate", with 20 levels in the vertical in the standard version.
The sea-ice component of CLIO is an updated version of the sea-ice model of Fichefet andMorales Maqueda (1997, 1999).It uses the same horizontal grid as the ocean model.Sensible heat storage and vertical heat conduction within snow and ice are determined by a three-layer model (one layer for snow and two layers for ice).Each grid box is partly covered by sea ice of uniform thickness (i.e., the model includes only one sea-ice thickness category) and open water (leads).Vertical and lateral growth/decay rates of the ice are obtained from prognostic energy budgets at both the bottom and surface boundaries of the snow-ice cover and in leads.When the load of snow is large enough to depress the snowice interface under the water level, seawater infiltrates in the model the entirety of the snow layer below the ocean surface and freezes there.This snow and the frozen seawater form then a new layer (snow ice), implying in the model an increase in sea ice thickness (Fichefet and Morales Maqueda, 1997).The parameterization of the surface albedo is taken from Shine and Henderson-Sellers (1985), with corrections for clear and overcast conditions as recommended by Greenfell and Perovich (1984).This albedo formulation takes into consideration the state of the surface (frozen or melting) and the thickness of the snow and ice covers.
For the momentum balance, sea ice is considered as a twodimensional continuum in dynamical interaction with the atmosphere and the ocean.The viscous-plastic constitutive law proposed by Hibler (1979) is used for computing the internal ice force.The ice strength is taken as a function of the ice thickness and compactness (Hibler, 1979).The physical fields that are advected are the ice concentration, the snow volume per unit area, the ice volume per unit area, the snow enthalpy per unit area, the ice enthalpy per unit area, and the brine reservoir per unit area.
The model equations are solved numerically as an initial value-boundary value problem by using finite difference techniques.A staggered spatial grid of type B is utilized.The heat diffusion equation for snow and ice is solved by means of a fully implicit numerical scheme, which avoids the development of numerical instabilities when the snow or ice thickness becomes small.The ice momentum balance is treated basically as in Zhang and Hibler (1997).A no-slip condition is imposed on land boundaries.The contribution of advection to the continuity equations is determined by making use of the forward time marching scheme of Prather (1986).This method is based on the conservation of the second-order moments of the spatial distribution of the advected quantities within each grid cell.It preserves the positiveness of the transported variables and has very small numerical diffusion.The advantage of employing this elaborate scheme is that, for a coarse resolution grid such as the one used here, it determines the location of the ice edge with a higher accuracy than the more conventional upstream schemes do.
A standard quadratic law is applied for calculating the stress at the ice-ocean interface.The heat flux from the ocean to the ice is computed by the parameterization of McPhee (1992), while the salt and freshwater surface exchanges are based on mass conservation.As CLIO includes a free surface, the exchanges of freshwater are represented by a vertical velocity at surface equal to precipitation -evaporation + runoff.However, for relatively subtle reasons linked to the way the free surface is represented in the model, applying such a straightforward method is not possible at the ice-ocean interface (Tartinville et al., 2001).As a consequence, all the mass exchanges between the ocean and sea ice are implemented as negative and positive salt fluxes, the freshwater fluxes being then virtual salt fluxes that have the same dilution effect as the corresponding freshwater exchanges.
All the model equations, parameters, and numerical schemes are described in detail in the user's guide of the CLIO model (http://www.astr.ucl.ac.be/index.php?page= CLIO%40Description).

VECODE: the continental biosphere component
The model for the terrestrial biosphere VECODE (VEgetation COntinuous DEscription model) (Brovkin et al., 2002;Cramer et al., 2001) was specifically designed with the purpose of interactive coupling with a coarse resolution atmospheric model for long-term simulations.It is a reduced-form dynamic global vegetation model (DGVM), which simulates changes in vegetation structure and terrestrial carbon pools on timescales ranging from decades to millennia.VECODE consists of three sub-models: (1) a model of vegetation structure (bioclimatic classification) calculates plant functional type (PFT) fractions in equilibrium with climate; (2) a biogeochemical model estimates net primary productivity (NPP), allocation of NPP, and carbon pool dynamics; (3) a vegetation dynamics model.PFTs (see, e.g., Prentice et al., 1992;Chapin et al., 1996 for the PFT concept) are used to describe the vegetation cover.For any given climate, there is a unique stable composition of PFTs corresponding to the climate (in this context, climate is understood as a long-term average of atmospheric fields).If climate changes, the vegetation model simulates the transition from the equilibrium for the previous climate to a new equilibrium with the new climate.The time scale of this transition is determined from the carbon cycle model.
A fractional bioclimatic classification (Brovkin et al., 1997) is developed in order to adapt discrete bioclimatic classifications (e.g.Life Zones by Holdridge, 1947, or BIOME by Prentice et al., 1992) for coarse resolution climate models.Two basic PFTs are used: trees and grasses.The sum of the tree fraction, f , and the grass fraction, g, is equal to the vegetation fraction, v; the rest corresponds to the desert fraction, d = 1 − v.These transient fractions are different from equilibrium fractions (vegetation in equilibrium with climate), denoted by f , v. Semi-empirical parameterizations are used for f and v: where G 0 is the growing degree-days above 0 (GDD0, i.e., the sum of the surface air temperature for all the days with a mean daily temperature higher than 0 • C), P r is the annual mean precipitation, c, a, a for , b, a des , b 2 and P 0 r are bio-climatic parameters (Table 3), G min is the minimum GDD0 for trees, and P min r is the minimum precipitation for vegetation.
Those equations are based on regularities of distribution of forest and desert in climatic space (Lieth, 1975) which have an ecophysiological basis (Woodward, 1987).The vegetation map of Olson et al. (1985) and an updated version (W.Cramer, personal communication, 1996) of the climate dataset of Leemans and Cramer (1991) were used in the validation procedure.
Carbon in the vegetation is aggregated into two compartments: a "fast" pool of green biomass (leaves), C 1 , and a 'slow' pool of structural biomass (stems, roots), C 2 .Dead organic matter is described by two pools: a "fast" compartment (woody residues), C 3 , and a "slow" compartment (humus), C 4 .Variables C i are simulated separately for trees and grass (represented here by φ).The dynamics of the carbon pools are integrated with an annual time step.Net primary productivity (NPP), , is simulated on an annual basis following the semi-empirical parameterization of Lieth (1975), which is often used for first-guess estimates on a global scale (Post et al., 1997).To find out the NPP per model grid cell, the NPP per square meter is multiplied by the vegetation fraction v, which is in most cases equal to one, and the land area.In dry subtropical regions, v is less than 1 (Eq.8), and this helps to correct the bias in the productivity per grid cell in these regions.Dependence of NPP on the atmospheric CO 2 concentration is taken into account by the biotic growth factor in a logarithmic form (den Elzen et al., 1995).
NPP allocation between green and structural biomass is estimated as a function of NPP, with increased allocation to C 2 relative to C 1 as NPP increases.This function was calibrated using an empirical dataset of NPP and carbon storage from about 500 sites in the northern Eurasia collected by Bazilevich (1993).The same data were used for calibrating parameterizations for the turnover time of biomass τ i , i = {1,2}, which is assumed to be a function of NPP.The turnover time of soil carbon τ i , i = {3,4} is a function of the mean annual temperature following the approach by Schimel et al. (1994).The annual maximum of Leaf Area Index (LAI) is assumed to be proportional to the green biomass.
To account for the subgrid-scale processes of vegetation succession, we apply linear ordinary differential equations for simulating the dynamics of the PFT fractions.The model implies that the vegetation cover reacts to any climate change with a relaxation towards a new equilibrium with a timescale determined by the turnover time of the structural biomass.For instance, if the climate becomes more humid and the equilibrium fraction of trees increases, then the trees become more successful in competing with grasses and occupy an additional fraction of land within the large grid cell with the timescale of tree growth.In the vicinity of an equilibrium,   the equation for the time development of vegetation is a linearized version of the evolutionary model for vegetation dynamics (Svirezhev, 1999) which accounts for competition between trees and grasses in the idealized form.With respect to the dynamics of the northern treeline under CO 2 -induced climate change, VECODE shows similar performance to other dynamic global vegetation models (Cramer et al., 2001).matter (DOM and POM), silica (Si), oxygen (O 2 ) as well as organic and inorganic carbon isotopes.The concentration of dissolved CO 2 at the sea surface is controlled by both physical and biological processes (solubility and biological pumps, respectively).Biology exerts a strong control on the surface CO 2 and is responsible for the fast transfer of carbon to the deep ocean.In a somewhat similar approach to that used in HAMOCC 3 (Maier-Reimer, 1993;Heinze et al., 2003), LOCH intends to reproduce export production (i.e.flux of organic carbon out of the surface ocean).The LOCH biological module should therefore not be viewed as a model of ocean ecosystems but rather as a model of biogenically mediated fluxes of constituents in the ocean.The basis for the export-production model is a pool of phytoplankton whose growth is driven by the availability of nutrients (DIP) and light.The evolution of phytoplankton biomass B follows:

LOCH: the ocean carbon cycle component
where µ B the actual growth rate is a function of temperature T , light L and inorganic phosphorus concentration (DIP): where µ Max is the maximum growth rate and K L , K T and K P are half-saturation constants for temperature, light and inorganic phosphorus concentration, respectively (Table 4).
The sink term R B B takes into account grazing and mortality and is defined as: in which m B and G represent the mortality and the maximum grazing rate, respectively.The use of a Michaelis-Menten-like formulation for grazing in Eq. ( 11) allows for a non-linear closure of the system, which is necessary in order to properly reproduce the productivity changes (Fasham, 1993).
Upon death, organisms feed the fast sinking particulate organic matter (POM) pool.The distribution of the POM flux with depth below the productive layers is governed by a power law z −α POM (Martin et al., 1987), with z the depth measured from the bottom of the euphotic zone.In LOCH, the actual vertical profile driving the distribution of the POM flux evolves according to the fraction of the total export production supported by silica shell building organisms; this is achieved by considering different values of α POM for diatoms and other species.
Below the productive layers the POM remineralizes as DIP or transforms into dissolved organic matter (DOM).DOM subsequently decays into DIP.The remineralization rate of organic matter (POM or DOM) depends on the oxygen availability.Anoxic remineralization occurs in O 2 -depleted re-gions but in a less efficient way than oxic processes.The remineralization rate is given by: where x either stands for POM or DOM.In Eq. ( 12), r o x and r a x represent the maximum oxic and anoxic remineralization rates, respectively.
It should be noticed that, although B and POM are prognostic variables, they are not subject to the 3-D transport.The rationale underlying this choice is that the characteristic timescale of these variables is much shorter than the one of interest in the context of climate studies.
The hard tissues (shells) are made up of CaCO 3 or opal, and their precipitation occurs concurrently with the softtissue formation.About half of the export production in the ocean is supported by diatoms (Nelson et al., 1995).Hence we discriminate between these organisms, which rely on silicon for their growth, and other species.A constant Si:P ratio is used to determine the export of opal accompanying the export production.The vertical distribution of biogenic silica below the productive layers upon the death of the organism writes e −βz , where β takes into account the influence of temperature on the dissolution rate with β = β d e κ d T .
Alkalinity and dissolved inorganic carbon are both needed to determine the concentration of dissolved CO 2 in surface waters as well as the CaCO 3 saturation level in deep waters.The total dissolved inorganic carbon (DIC) represents the sum of dissolved CO 2 , bicarbonate and carbonate.Total alkalinity is a measure of the acid neutralizing capacity of seawater, as defined by Dickson (1981).However, in order to reduce the computing time, this definition is simplified by retaining only the essential contributions (bicarbonate, carbonate and borate).The error resulting from the neglect of phosphorus and silica contributions to Alk is far smaller than other uncertainties inherent to climate modelling.The constants required to determine the various chemical equilibria in seawater are expressed on the seawater pH scale.When needed, transformation from the free pH scale to the seawater pH scale is performed with the help of formulations from Millero (1995) and Dickson and Riley (1979).The system is fully determined by using dissociation constants for water from Millero (1995), for borate from Dickson (1990) and for carbonates from Dickson and Millero (1987).
This model, similarly to other simplified models of the ocean carbon cycle (e.g.Maier-Reimer, 1993;Najjar et al., 2007), assumes the stoichiometric constancy of organic material.It means that POM, DOM as well as biomass composition are in Redfield ratio.The sources and sinks terms for DIC and Alk are therefore simply derived from the biological fluxes with the help of the phosphorus to carbon Redfield ratio of Anderson and Sarmiento (1994) and the nitrogen to phosphorus ratio of Redfield et al. (1963).One important factor for the carbon cycle is the rain ratio, which is the amount of organic carbon assimilated during photosynthesis divided by that of inorganic carbon incorporated into shells.The rain ratio R CaCO 3 in LOCH depends on the availability of silica, the latter determining which type of shells will be preferentially built.The influence of temperature and the ubiquity of calcareous organisms are also included in the parameterization of this process.R CaCO 3 is defined as: with R CaCO 3 ≤ R Max CaCO 3 , the maximum rain ratio.Expression (13) includes the following parameters or variables: r CaCO 3 the minimum rain ratio, K CaCO 3 half-saturation constant for CaCO 3 precipitation ( • C), Zoo the rain ratio associated to zooplankton, Phy the rain ratio associated to non siliceous phytoplankton, and f DIA the fraction of siliceous phytoplankton, f DIA ∈ [0,1].A constant fraction f CaCO 3 of calcium carbonate shells is also assumed to be made of aragonite, which is more soluble than calcite.
The dissolution of shells occurs in the deepest ocean layer under the production area at a rate controlled by the CaCO 3 saturation level.Hence LOCH implicitly includes carbonate compensation mechanisms.The expressions for the solubility of calcite and aragonite are from Mucci (1983) and Millero (1995), while the coefficients for the pressure dependence of the chemical equilibrium constants are from Millero (1995).
Some organic matter and shells escape remineralization or dissolution, and are permanently preserved in sediments.On the other hand, river input of alkalinity, silica, organic matter and carbon constitutes a net source for the ocean.In the case of an equilibrium run, this source exactly compensates the permanent preservation in sediments.The main rivers of the world and their respective importance are taken into account in this process.
The magnitude of the air-sea flux of a gas depends on the difference of its partial pressure between the two media, with an exchange rate given by the product of the solubility and the piston velocity.The solubilities are taken from Wanninkhof (1992) for O 2 and from Weiss (1974) for CO 2 .The piston velocity follows the empirical formulation proposed by Wanninkhof (1992), which relates it to the squared wind velocity and the Schmidt number.The latter is gasdependent and is calculated according to Wanninkhof (1992).An additional term accounts for the chemical enhancement of CO 2 exchange at low wind speeds and high temperatures (Wanninkhof and Knox, 1996).
LOCH also includes an atmospheric module which simulates the evolution of the various gases in the atmosphere.It is based on a 1D diffusion equation in the meridional direction, i.e., one implicitly assumes instantaneous mixing in the zonal and vertical directions.Hence the transport in the atmosphere of a constituent with concentration C (ppmv) obeys to: where t is the time and y the position in the meridional direction.The diffusion coefficient K y (m 2 s −1 ) is homogeneous within each hemisphere and allows mixing within a few weeks.A lower value of K y is used at the equator so that inter-hemispheric mixing occurs with a characteristic timescale of 2 years (Bacastow and Maier-Reimer, 1990).P C includes local sink terms where relevant, e.g., radioactive decay for 14 C. F C represents the exchange of gases between the atmosphere on the one hand and the ocean and the continental biosphere on the other hand.If applicable, F C may also include other sources (e.g., anthropogenic emissions).The gases taken into account are carbon dioxide CO 2 , oxygen O 2 as well as the two isotopic forms 13 CO 2 and 14 CO 2 .
Equation ( 14) is discretized with a constant spatial step, at the same resolution as CLIO (3 • ), and with the same time step as LOCH.The atmospheric module offers two options for the study of the carbon cycle: either the concentrations are prescribed in the atmosphere (diagnostic mode) or the concentrations evolve according to the various exchange processes as described above (prognostic mode).

AGISM: the polar ice sheet component
AGISM (Antarctic and Greenland Ice Sheet Model; Fig. 6) consists of two three-dimensional thermomechanical ice dynamic models for each of the polar ice sheets.Both models are based on the same physics and formulations, however with the major distinction that the Antarctic component incorporates a coupled ice shelf and grounding line dynamics.Ice shelf dynamics is missing from the Greenland component as there is hardly any floating ice under present-day conditions, and this can be expected to disappear quickly under warmer conditions.Having a melt margin on land or a calving margin close to its coast for most of its glacial history, ice shelves probably played a minor role for Greenland also during colder conditions.
Both polar ice sheet models consist of three main components which respectively describe the ice flow, the solid Earth response and the mass balance at the ice-atmosphere and ice-ocean interfaces (Huybrechts and de Wolde, 1999;Huybrechts, 2002; to which papers the interested reader is referred to for a full overview of all equations and model formulations as well as for additional references).Figure 6 shows the structure of the model.At the heart of these models is the simultaneous solution of two evolutionary equations for ice thickness and temperature, together with diagnostic representations of the ice velocity components.Conservation of ice volume and heat is expressed as: where H is the ice thickness, v the depth-averaged horizontal velocity field, M the mass balance and t the time.The thermodynamic equation considers heat transfer to result from vertical diffusion, 3-D advection, and internal frictional heating caused by ice deformation (φ).The inclusion of heat conduction in the bedrock gives rise to a variable geothermal heat flux at the ice sheet base depending on the thermal history of the ice and rock.T i and T m are ice and rock temperatures, respectively, and k, c, and ρ are thermal conductivity, specific heat capacity and density for respectively ice and rock (subscripts "i" and "m").The thermal parameters of ice, k i and c i , also take into account the temperature dependence, which effect is not negligible as their values may change by up to 30% for ice temperatures ranging between 0 • C and −50 • C (Huybrechts, 1992).Main parameter values are given in Table 5.
In grounded ice, the flow results from both internal deformation and sliding over the bed in places where the temperature reaches the pressure melting point and a lubricating water layer is present.Ice deformation in the ice sheet domain results from vertical shearing, most of which occurring near to the base.Longitudinal deviatoric stresses (i.e., the difference between the total stress and the hydrostatic stress -actually "cryostatic" in this case) are disregarded according to the widely used "Shallow Ice Approximation" (e.g., Hutter, 1983).This does not treat the rapid component of the otherwise badly understood physics specific to fast-flowing outlet glaciers or ice streams.A flow law of "Glen type" is used with exponent n = 3 (Glen, 1955;Paterson, 1994).For the sliding velocity, a generalized Weertman relation is adopted (Weertman, 1964), taking into account the effect of the subglacial water pressure.Ice shelves are included by iteratively solving a coupled set of elliptic equations for ice shelf spreading in two dimensions, including the effect of lateral shearing induced by sidewalls and ice rises.At the  Isostasy is taken into account for its effect on bed elevation near grounding lines and marginal ablation zones, where it matters most for ice sheet dynamics, and because isostasy enables ice sheets to store 25-30% more ice than evident from their surface elevation alone.The bedrock adjustment model consists of a viscous asthenosphere, described by a single isostatic relaxation time, which underlies a rigid elastic plate (lithosphere).In this way, the isostatic compensation accounts for the effects of loading changes within an area several hundred kilometers wide, giving rise to devia-tions from local isostatic equilibrium.The downward deflection w of the Earth caused by the weight of ice sheets and oceans is determined by the rigidity of the lithosphere and the buoyancy of the mantle, and is a solution of: where g is the acceleration due to gravity, h is the bedrock elevation, ρ w is the density of sea water, and H sl is the eustatic sea-level stand relative to present-day.The standard value for the flexural rigidity D (cf.  which is asymptotically attained using a relaxation formulation schematically representing the Earth's mantle: where the unloaded surface elevation h 0 has been determined by assuming that the Earth is in present-day isostatic equilibrium with both the ice and water loading, and τ is the asthenospheric decay timescale.The isostatic treatment produces results close to those from more sophisticated visco-elastic Earth's models, while at the same time being much more efficient in terms of computational cost.The loading takes into account contributions from both ice and ocean water within the respective grids, but ignores any ice loading changes beyond the Greenland and Antarctic continental areas. For both ice sheets, calculations are made on a 10 km × 10 km horizontal resolution, with 31 vertical layers in the ice and another 9 layers in the bedrock for the calculation of the heat conduction in the crust (Fig. 7).The vertical grid in the ice has a closer spacing near to the bedrock where the shear concentrates.Rock temperatures are calculated down to a depth of 4 km, which is deemed sufficient to capture most of the effect of temperature changes on glacial-interglacial timescales.This gives rise to between 1.85 and 12.6 × 10 6 grid nodes for Greenland and Antarctica, respectively.Geometric datasets for surface elevation, ice thickness and bed elevation incorporate most of the recent observations up to 2001, such as ERS-1 derived satellite heights, BEDMAP and EPICA pre-site survey Antarc-tic ice thicknesses and the University of Kansas collection of airborne radio-echo-sounding flight tracks over Greenland (Huybrechts and Miller, 2005).The grids correspond to those discussed in Huybrechts and Miller (2005).The finitedifference schemes are implicit in time, either alternatively in the x and y directions for the mass conservation equation or only along the vertical for the thermodynamic equations.The 10 km horizontal resolution substantially improves the representation of the fast-flowing outlet glaciers and ice streams which are responsible for the bulk of the ice transport towards the margin.Other physics specific to these features such as higher-order stress components or subglacial sediment characteristics are not included, in common with the current generation of three-dimensional ice-sheet models.
Interaction with the atmosphere and the ocean is effectuated by prescribing the climatic input, consisting of the surface mass balance (accumulation minus ablation), surface temperature and the basal melting rate below the ice shelves surrounding the Antarctic component.The mass balance model distinguishes between snow accumulation, rainfall and meltwater runoff, which components are all parameterized in terms of temperature.The melt-and runoff model is based on the positive degree-day method and is identical to the recalibrated version as described in Janssens and Huybrechts (2000).Following what has become standard practice in large-scale ice sheet modelling, the melting rate is set proportional to the yearly sum of positive degree days at the surface.The expected sum of positive degree days (EPDD) can conveniently be evaluated as: + max 0, where the standard deviation σ is for temperature with respect to the monthly mean surface temperature T sur mon to account for the daily cycle and weather fluctuations.dt represents monthly time steps and the factor 30.44 is the mean amount of days in one month (365.2422/12).The expected number of positive degree days (EPDD) is used to melt snow and ice.Meltwater is at first retained in the snowpack by refreezing and capillary forces until the pores are fully saturated with water, at which time runoff can occur.This method to calculate the melt has been shown to be sufficiently accurate for most practical purposes.It moreover ensures that the calculations can take place on the detailed grids of the ice-sheet models so that one can properly incorporate the feedback of local elevation changes on the melt rate, features which cannot be represented well on the generally much coarser grid of a climate model.The melt model is also implemented for Antarctica, but since current summer temperatures remain generally below freezing, melt amounts are currently negligible there.Because of their very low surface slopes, it is further assumed that meltwater produced on the surface of the Antarctic ice shelves, if any, refreezes in situ at the end of the summer season, and therefore does not run off to the ocean.Below the ice shelves, a uniform melting rate is applied which magnitude is linked to the heat input into the cavity, as explained in Sect.2.7.

The iceberg model
LOVECLIM has an optional iceberg module which has been activated only in a few studies up to now (Jongma et al., 2009;Wiersma and Jongma, 2009).It will not be used in the experiments discussed in Sect.3. In such a case, the transport of the icebergs is not taken into account so that the latent and freshwater fluxes caused by iceberg melting takes place at the oceanic grid point where the icebergs leave the ice sheet.However, the iceberg model is briefly described here for completeness as it is part of the code.
The basic equation for horizontal motion of the icebergs is: for an iceberg with mass M (kg) and velocity V i (m s −1 ), subject to Coriolis force −Mf k ×V i , air drag F a , water drag F w , sea-ice drag F s , horizontal pressure gradient force F p and wave radiation force F r .The general drag relationship is given by (Smith, 1993): where x refers to air (a), water (w) and sea ice (s) respectively, with medium density ρ x (kg m 3 ) and drag coefficient C x (C a = 1.3, C w = 0.9 (following Smith, 1993, who selected drag coefficients to optimize the fit of the modelled to the observed tracks, but good results were also obtained by Smith with a priori values of 1.0) and C s = C w (Bigg, et al., 1997;Gladstone, et al., 2001; see Table 6).A x is the cross-sectional area of the iceberg perpendicular to the stressing medium x, which has velocity V x (m s −1 ).In accordance with Ekman theory (Bigg et al., 1997), the icebergs are assumed to be travelling with their long axis parallel to the surrounding water and sea-ice flow and at an angle of 45 • to the wind flow (A w = A s = 1 and A a = |1.5sin(45 The wave radiation force is given by (Smith, 1993): where L the length of the iceberg perpendicular to incident waves, which have amplitude a and are assumed to have the same direction as wind velocity V a .
The horizontal pressure gradient force exerted on the water volume that the iceberg displaces F p (Bigg et al., 1997) is taken from the free surface ocean model's variable at the iceberg location (Deleersnijder and Campin, 1995).To obtain the strength of the forcing fields at the iceberg's location, linear interpolation from the four surrounding grid corners of the climate model is used.
The icebergs are weakly repelled from the coast using a velocity of 0.003 m s −1 in an orthogonal direction when their keel exceeds water depth.They are assumed to remain tabular, maintaining a constant length to width ratio of 1:1.5 (see Bigg et al., 1997).Keel shape or other turbulence related effects are not accounted for.Added mass due to entrained melt water is neglected.Due to real icebergs inertial rotation and individual shapes, this approach can only be considered as a rough approximation.It describes the general behaviour of icebergs but cannot be expected to work well for individual bergs.The drag coefficients for water stress acting along the lower surface of the iceberg and atmospheric wind stress acting along the top surface are deemed negligibly small (G.R. Bigg, personal communication, 2004).There is no direct interaction between icebergs.
The iceberg thermodynamics must be accounted for in any long-term simulation of its trajectory, since the iceberg mass and shape change due to melting.The iceberg melt is simplified to basal melt, lateral melt and wave erosion (Bigg, et al., 1997).The basal turbulent melting rate (Weeks and Campbell, 1973) is a function of the difference between iceberg (T i = −4 • C) and water temperature (T w ).
The lateral melt due to buoyant convection along the sides of the iceberg is given by an empirical relationship (Eltahan et al., 1983) as a function of water temperature T w ( • C) of the corresponding ocean layer in the local grid cell.Wave erosion (Bigg et al., 1997) M waves = 0.5S s ( 26) is a function of the sea state S s (based on the definition of the Beaufort scale): where V a is the wind speed (km h −1 ).Iceberg deterioration by atmospheric and radiation effects is considered negligible (Loset, 1993).Break-up of icebergs is not modelled.When the ratio between iceberg length L and height H exceeds a criterion of stability, the icebergs are allowed to roll over (Bigg, et al., 1997): To achieve climatic coupling, the freshwater and latent heat fluxes associated with the iceberg melt are added to the corresponding ocean layer of the local grid cell.Direct feedbacks from the icebergs to the atmosphere are relatively small (e.g., Loset, 1993) and are not accounted for.

Coupling the different components
The equations of the atmospheric and ocean models are solved on different grids.An interpolation is thus required during the transfers between the two models.CLIO provides ECBilt with the sea surface temperature, the sea-ice temperature, the fraction of sea ice in each ocean grid cell and the sea-ice and snow thicknesses (in order to compute the snow and sea-ice albedo in ECBilt).ECBilt gives to CLIO the wind stresses over the ocean and sea ice, the shortwave and net heat flux over the sea-ice and ocean fractions of the grid box, and the solid and liquid precipitation (including runoff, evaporation and sublimation).In order to have a conservative interpolation, the surface covered by land, ocean and sea ice is exactly the same in ECBilt and CLIO.This is achieved by decomposing the surface of each atmospheric grid box in three parts.Those fractions are interpolations on ECBilt grid from the one in CLIO.CLIO determines thus the location of the coastlines, and more generally of the land sea mask for the all the components (Fig. 8).The land and ocean fractions and the ocean bathymetry are fixed trough time in LOVECLIM (but can of course be different for simulations devoted to different periods).As a consequence, any change in coastal geometry, for instance implied by a sea level rise, cannot be taken into account explicitly in the model.
No flux correction on stress and heat fluxes is applied between ECBilt and CLIO.However, as precipitation rates in the Atlantic and the Arctic are significantly overestimated in ECBilt, they are reduced by 8.5% and 25%, respectively, before being transmitted to CLIO in order to avoid too large an ocean drift.In order to conserve mass, the excess water is dumped into the North Pacific where ECBilt underestimates precipitation, ensuring an additional constant flux per unit area in the region located between the Bering Strait and the equator.
LOCH and CLIO run on the same grid (Fig. 4).The time step for solute transport in LOCH is the same as the time step for tracer transport in CLIO, thus eliminating the need for any interpolation procedure.However, LOCH uses a numerical scheme for advection which differs from the one of CLIO.The reason for this difference is to be found in the non monotonic behaviour of the CLIO advection scheme.
Transport in LOCH is based on two-dimensional and three-dimensional fields provided by CLIO: downsloping flows and heights, salt and freshwater fluxes at the sea surface, current velocities and vertical and horizontal diffusivities.The chemical constants, the gas exchange coefficients and other parameters of LOCH are computed from the temperature and salinity fields provided by CLIO.The piston velocity is determined from the wind field simulated by ECBilt.The growth rate of the phytoplankton biomass is set according to the same amount of available light at the sea surface (under the ice in ice-covered areas) as in CLIO; we however use a different extinction coefficient with depth.The sea-ice areal coverage modelled by CLIO is also taken into consideration in the calculation of the air sea gas fluxes.
The atmospheric component of LOCH computes the atmospheric CO 2 evolution in zonal bands which are equally spaced in latitude.A spatial interpolation procedure is thus required to transfer to LOCH the annual mean values of the CO 2 fluxes between atmosphere and continents computed by VECODE as well as those from LOCH since over the North Atlantic the grid of CLIO is not regular with latitude.By combining the carbon fluxes from the continents and from the ocean, LOCH computes a globally averaged,  annual mean atmospheric CO 2 concentration which is then transmitted to ECBilt and VECODE, where it impacts on the radiative transfer and fertilization, respectively.
The key atmospheric variables needed as input for AG-ISM are monthly surface temperature and annual precipitation.Because the details of the Greenland and Antarctica surface climate are not well captured on the ECBilt coarse grid, these boundary conditions consist of present-day observations as represented on the much finer AGISM grid onto which climate change anomalies from ECBilt are superimposed.Monthly temperature differences and annual precipitation ratios, computed against a reference climate corresponding to the period 1970-2000 AD (PD), are interpolated from the ECBilt grid onto the AGISM grid and added to and multiplied by the parameterized surface temperatures and observed precipitation rates, respectively.The perturbation ("delta") method for temperature is represented by: T sur mon (φ,λ,t) = T sur ECBilt (φ,λ,t) − T sur ECBilt (φ,λ,PD) +T sur par (φ,λ,PD) − γ H sur ECBilt (φ,λ,t) − H sur ECBilt (φ,λ,PD) (29) where the monthly mean surface temperature is specified as a function of time t and location (φ, λ), the first term on the right-hand side is the mean monthly temperature anomaly from ECBilt, the subscript par denotes the parameterized surface temperature in the ice sheet model, and an additional correction is required to correct for the elevation temperature change in ECBilt (last term) to avoid double counting.γ is a prescribed atmospheric lapse rate.
The treatment of precipitation is similar to that of temperature, except that the ratio is used and not the difference.This is because using the same form of Eq. ( 29) for precipitation might introduce "negative precipitation" into the climate forcing, which has no physical basis.The appropriate relation reads: where the yearly precipitation rate distribution is also given as a function of time and location, and period 1970-2000.The subscript cli refers to the observed precipitation climatology over the ice sheets and is representative for the same reference period.This approach avoids systematic errors in the absolute EC-Bilt fields and ensures that some processes, such as the melting taking place at the ice sheet margin over a spatial extent narrower than the atmospheric model resolution, can be adequately represented.
The ocean heat flux at the base of Antarctic ice shelves is also calculated in perturbation mode based on a parameterization proposed by Beckmann and Goosse (2003): where M is the basal melt rate, Q net an estimate of the total heat flux entering the ice shelves integrated all along the perimeter of Antarctica, and A the total area of Antarctic ice shelves.Here the subscripts t and 0 refer to the actual model time and the reference time taken as 1500 AD, respectively.
In this approach the melt rate below the ice shelves depends on the net heat input from the oceans into the cavity below the ice shelves.The total melt volume is proportional to changes in the net integrated ocean heat input but inversely proportional to the area of the ice shelves.The underlying assumption is that much of the water in the cavity is recycled locally forming a semi-closed circulation cell.Q net is estimated directly from the mean ocean temperature around Antarctica.After performing mass balance and ice dynamic computations, AGISM transmits the calculated changes in orography and land fraction covered by ice to ECBilt and VECODE.This involves accounting for the albedo of the ice but also for the monthly snow cover over ice-free areas of Greenland.Land cover changes over Antarctica are not expected for most periods being studied.In addition, AGISM provides CLIO with the geographical distribution of the annual mean surface freshwater flux resulting from ice sheet runoff, iceberg calving, runoff from ice-free land and basal ice melting.The transfer of data from AGISM to ECBilt is rather straightforward since the grid cells of ECBilt are much larger than the AGISM ones.Each AGISM grid cell is associated with an ECBilt grid cell, and an area average is made to determine the value of a specific variable on the ECBilt grid.For the interpolation of data from the ECBilt grid to the AGISM grid, we opted to first transform the AGISM points on the ECBilt grid and subsequently apply a Lagrangian interpolation.The selected interpolation is a third-order Lagrange polynomial.Four ECBilt grid points are taken into account in latitude and four in longitude to determine the polynomial providing the variable value at each particular AGISM grid point.
Regarding the coupling between AGISM and CLIO, a simple procedure was set up to allocate the total freshwater flux from AGISM to the respective surface ocean grid boxes of CLIO that border Greenland and Antarctica.It must also be mentioned that the latent heat associated with iceberg melt-ing is pumped from these grid boxes.The coupling technique described above leads to heat and water losses/gains in the coupled model.Due to the perturbation method employed and the use of a Lagrangian interpolation, the amount of water received by AGISM in the form of precipitation is not equal to the amount of water leaving ECBilt.Biases are of the order of 10% to 25% of the total runoff from Antarctica and Greenland, respectively.Similarly, the heat available in ECBilt for the ice sheet melting differs from the one in AGISM.Flux adjustments are therefore necessary to ensure strict conservation of heat and water.These are applied uniformly in a given ocean area around each ice sheet.The water correction is treated as an additional freshwater flux and the heat correction as an additional latent heat flux associated with iceberg melting.This ensures the closure of the heat and water balances in the coupled system.

Evaluation of model performance
As LOVECLIM is a model of intermediate complexity, it cannot be expected to reproduce all the observations with the same skill and the same level of detail as a GCM.Indeed, previous studies have underlined some clear and strong model biases in LOVECLIM results.Some of those biases are directly linked to the model formulation and reducing significantly their amplitudes can only be achieved by modifying fundamental model assumptions.This would then be at the expense of some of the main advantages of LOVE-CLIM.As it is not our goal here to modify the philosophy behind the model development, such biases are still present in version 1.2.
Nevertheless, it is instructive to document the regions (and variables) where the discrepancies are the largest and the ones where the agreement between model results and observations is satisfactory because it is an important element when interpreting results of experiments performed with the model.In the following sections, we will thus describe briefly the mean state of the model for present-day conditions and then discuss the model behaviour for 4 key periods: the last decades, the last millennium, the mid-Holocene (6 ky BP) and the Last Glacial Maximum (LGM, 21 ky BP).The last two periods are standard ones in the Paleoclimate Modelling Intercomparison Project (PMIP, see for instance Braconnot et al., 2007).
Idealized experiments have also been performed with the model.They are not described here but it is useful to mention that when the CO 2 concentration is doubled compared to preindustrial conditions, the surface temperature increases by 1.9 • C after 1000 years of integration in LOVECLIM (with fixed ice sheets), giving an estimate of the model climate sensitivity.This is at the lower end of the range of values obtained from GCM results (e.g., Randall et al., 2007).In another experiment, under pre-industrial conditions, a freshwater flux of 0.1 Sv (1 Sv = 10 6 m 3 s −1 ) has been imposed in  the North Atlantic during 1000 years, inducing a 30% decrease of the maximum of the overturning streamfunction in the North Atlantic (see below for a description of this variable).This indicates that LOVECLIM1.2 is slightly more sensitive to freshwater perturbations than an early version of ECBILT-CLIO (Rahmstorf et al., 2005).

Present-day mean climate
In order to compare the model results with recent observations, a transient simulation has first been performed with LOVECLIM over the last 1500 years using all the components of LOVECLIM except the iceberg model.The average over the last decades of this simulation is used first to evaluate the model behaviour for present-day conditions.This simulation will also be analyzed in Sect.3.2 and 3.3 to study simulated changes during the past decades and the past millennium, respectively.
The initial conditions for LOCH, VECODE, ECBilt and CLIO come from a quasi equilibrium run, of several thousand years duration, corresponding to the forcing applied in 500 AD.For AGISM, as the ice sheets cannot be considered in quasi-equilibrium with the climate at that time, the initial conditions are obtained from a run of AGISM in uncoupled mode covering the last glacial-interglacial cycles and the Holocene up to 500 AD.
During the transient experiments, long-term changes in orbital parameters follow Berger (1978) and the long-term evolutions of non-CO 2 greenhouse gas concentrations are imposed.The variations in the emission of CO 2 from fossil fuel burning are derived from Marland et al. (2003).The influence of anthropogenic  sulfate aerosols is represented through a modification of surface albedo (Charlson et al., 1991).The forcing due to anthropogenic land-use changes (including both surface albedo and surface evaporation and water storage) is applied as in Goosse et al. (2005a) following Ramankutty and Foley (1999).Finally, natural external forcing due to changes in solar irradiance and explosive volcanism are prescribed following the reconstructions of Muscheler et al. (2007) and Crowley et al. (2003), respectively.The total solar irradiance changes have been scaled to provide an increase of 1 W m −2 between the Maunder minimum (late 17th century) and the late 20th century.This roughly corresponds to a threefold reduction in amplitude compared to some previous simulations conducted with the model (e.g., Goosse et al., 2005) but is in better agreement with recent reassessments (Lean et al., 2002;Foukal et al., 2006).
When comparing the mean climate over the last decades of this simulation to observations, we see that LOVECLIM1.2reproduces reasonably well the main characteristics of the observed surface temperature distribution (Fig. 9).For instance, the zero degree isotherm is quite close to the observed one in both hemispheres, with a more or less constant latitude in the Southern Hemisphere and a wavy structure in the Northern Hemisphere that displays a more northern position on continents than over the oceans.The strong differences at mid-and high-latitudes between the cold western part of the Atlantic compared to the warmer eastern part is also clearly seen in both model results and observations.In the tropics, the model is too warm, with a 25 • isotherm located too far away from the equator and an overestimation of the temperature over the continents.Furthermore, the temperature is much too high in the Eastern Pacific.As a consequence, the temperature gradient between the Eastern and Western Pacific is underestimated, reaching in the model about 2.5 • C in annual mean at the equator compared to more than 3.5 • C in observations interpolated on the model grid.The simulated zonal mean precipitation has roughly the right magnitude in nearly all the latitude bands (Fig. 10).However, the simulated pattern is much too symmetric between the hemispheres.In particular, the model is not able to reproduce the clear and strong absolute maximum in the Northern Hemisphere associated with the Intertropical Convergence Zone, displaying its maximum near the equator.A similar problem is also seen in many other EMICs (e.g.Petoukov er al., 2005).Furthermore, the precipitation at the observed local minima around 20 • S and 30 • N is clearly overestimated by the model.At some latitudes, the model error can reach 50% of the precipitation in zonal mean.
In both hemispheres, the large-scale structure of the nearsurface circulation (Fig. 11) is well reproduced by the model with, as expected, a general decrease of the geopotential height with latitudes and local minima in the North Atlantic, the North Pacific and in a belt around 70 • S. Except for the Aleutian low, the model underestimates the gradients in both hemispheres, leading to simulated winds weaker than the observed ones.Furthermore, the simulated minimum of the geopotential height in the North Atlantic is located too far eastwards, close to Baffin Bay, while the observations have their minimum near Iceland, inducing a wrong wind direction east of Greenland.LOVECLIM is able to simulate reasonably well the seaice extent in both hemispheres (Fig. 12).In the Northern Hemisphere, the sea-ice edge is very close to the observed one in the Pacific sector, both during summer and winter.In the Atlantic sector, the simulated sea-ice edge is too far northwards in the Baffin Bay and Labrador region in winter, while, in summer, the sea-ice extent is too large.The amplitude of the seasonal cycle of the sea-ice concentration is thus clearly too weak in this region in the model.In the western part of the North Atlantic, the model tends to slightly overestimate the sea-ice concentration, both in summer and in winter.The sea-ice extent is also slightly overestimated in the Southern Ocean in both seasons.Two exceptions are the regions west of the Antarctic Peninsula in summer and off East Antarctica around 45 • E in winter, where the model underestimates the sea-ice extent.
The maximum of the overturning streamfunction in the North Atlantic reaches 22 Sv, with an export towards the Southern Ocean of 13 Sv (Fig. 13).Deep convection in the model occurs both in the Greenland-Norwegian Sea and the Labrador Sea, as observed over the last decades.The maximum of the deep cell close to Antarctica has a value of 12 Sv, while 17 Sv are transported northward close to the bottom in   the global ocean.The shallow wind-driven cells in the tropics are associated with a total upwelling close to the equator of about 50 Sv.All those values are in relatively good agreement with the data-based estimates and in the range of the values given by other models (Ganachaud and Wunsch, 2000;Karsten and Marshall, 2002;Gregory et al., 2005;Rahmstorf et al., 2005).
As the model tends to overestimate precipitation in the tropics, the vegetation cover is also overestimated in those regions (Fig. 14).The vegetation fraction is also too large at high latitudes, mainly because of an overestimation of the temperature over the continent.By contrast, LOVECLIM has too low a vegetation cover in some regions of Australia and Southern America around 30 • S.

The last decades
In response to the forcing applied, the model simulates a clear increase in the global mean surface temperature (Fig. 15) and in the atmospheric CO 2 concentration (Fig. 16) over the 20th century and the beginning of the 21st century.The model is also able to reproduce the observed intensification of the warming trend over the last decades (Table 7).However, the model significantly underestimates the magnitude of this warming.This can be partly explained by the too large increase in ocean heat content in the model, the ocean playing apparently a larger buffering role in the model than in observations.This is a standard model bias that is discussed in detail in Loutre et al. (2010).
For the atmospheric CO 2 concentration, the model is quite close to observations (Fig. 16), with only a slight underestimation of the trend of the last 50 years (Table 7).The observed decrease in the summer sea-ice extent in the Arctic is   Jones et al. (2001), MBH1999 to Mann et al. (1999), MJ2003 to Mann and Jones (2003), MSH2005 to Moberg et al. (2005), PS2004 Pollack and Smerdon (2004), reference level adjusted following Moberg et al. (2005), RMO2005, Rutherford et al. (2005).also reasonably well simulated by the model (Table 7).This underlines that the underestimation of the warming seen at the global scale is mainly related to too weak a response of the model at low latitudes (Driesschaert, 2005). 5Uncertainties on the LOVECLIM results are estimated from the standard deviation of an ensemble of 5 experiments performed with the model using the same forcing but slightly different initial conditions.

The past millennium
The temperatures simulated over the past millennium display decadal to multi-centennial variations as well as a weak cooling trend over the period 1000-1850 before the large warming of the industrial era (Fig. 15).This is broadly consistent with the various reconstructions available as well as with previous model simulations.However, the long-term cooling between the period around 1000-1200 and the one around 1600-1850 is weaker here than in previous simulations performed with the model (e.g., Goosse et al., 2005).This is mainly due to the weaker amplitude of the variations of solar irradiance applied here.The simulated CO 2 concentration is quite stable in the model over the pre-industrial period (Fig. 16).As a consequence, the model is not able to reproduce the small decrease in CO 2 concentration between the periods 1200-1400 and 1700-1800 suggested by the observations (for a recent discussion of this feature, see for instance Frank et al., 2010).
The changes in the volume of ice sheets as simulated by the standard model over this period are relatively weak.Over Antarctica, the ice volume increases by 0.1% in 1000 years, while it decreases by about 1% over Greenland over the same period (Fig. 17).A small acceleration of the retreat is also seen in Greenland over the last decades.It is hard to tell at this stage if the trend in both curves is due to a longterm response of the ice sheets to past climate changes or results from a small drift introduced by the coupling procedure.Anyway, the simulated changes are small and can be neglected when analyzing future changes as they are at least an order of magnitude smaller than the ones simulated by the model for the 21st century and beyond (Driesschaert et al., 2005;Swingedouw et al., 2008).For analyzing past changes over several thousand years, the problem needs to be considered more carefully, but such simulations have not yet been carried out with LOVECLIM including all its components.(Indermühle et al., 1999), Law Dome (Etheridge et al., 1998), Siple (Neftel et al., 1994), South Pole (Siegenthaler et al., 2005), D47 (Barnola et al., 1995), D57 (Barnola et al., 1995), Drauning Maud Land (DML, Siegenthaler et al., 2005).et al., 1999), Law Dome (Etheridge et al., 1998), Siple (Neftel et al., 1994), South Pole (Siegenthaler et al., 2005), D47 (Barnola et al., 1995), D57 (Barnola et al., 1995), Drauning Maud Land (DML, Siegenthaler et al., 2005).

Mid-Holocene conditions
For the mid-Holocene simulation, the orbital parameters have been set at the value corresponding to 6 ka BP and the methane concentration has been reduced to 650 ppbv.The concentrations values for all other greenhouse gases (including CO 2 ) are the same as for pre-industrial conditions.All the other conditions have been chosen identical to preindustrial ones and a quasi-equilibrium multi-millennia run has been carried out.For this simulation experiment, LOCH and AGISM were not activated.
In response to the larger summer insolation, LOVE-CLIM1.2 simulates an increase in JJAS (June-July-August-September) surface air temperature at 6 ka BP over the continents in the Northern Hemisphere and over the Arctic   In this particular example, the Greenland ice volume budget is equivalent to a positive sea level contribution of about 10 cm over the entire period.The Antarctic ice volume budget is slightly positive but cannot be directly related to sea level change because of ice grounded below sea level.Variability in both indices on centennial timescales arises from the climate forcing and dynamical ice-climate interactions.The modelled trend is not a robust feature of AGISM, but contains a significant component from the model coupling procedure at 500 AD and the specific model parameters selected for ECBilt and CLIO.compared to present-day conditions (Fig. 18).The Southern Ocean is also warmer with a local temperature maximum increase of ∼4 • C between 30 • E-40 • E. By contrast, some regions show a small cooling such as seen in Africa just north of the equator, in the Middle East and west of the Japan coast.The JJAS mean precipitation (Fig. 19) produced by the LOVECLIM1.2model, captures well the Mid-Holocene characteristic increase over Northern Africa and in the Middle East, associated with an increase of vegetation there (Fig. 20).In the northeast of South America there is also  an increase of ∼1 mm day −1 .Just southward of the equator, there is less precipitation over the ocean in the mid-Holocene than today.All those results agree reasonably well with those from other models participating in the PMIP2 intercomparison (Braconnot et al., 2007), albeit tropical ocean feedbacks are relatively weak due to the quasi-geostropic approximation in the atmospheric component ECBilt (Zhao et al., 2005).

The last glacial maximum
In order to simulate the last glacial maximum climate, the orbital parameters have been modified to the values corresponding to 21 ka BP and CO 2 , methane and NO 2 concentrations were set to 185 ppmv, 350 ppbv and 200 ppbv, respectively, following the PMIP2 protocol.In addition, the topography of the ice sheets and the geometry of the coastlines have been imposed according to the ICE-5G reconstruction (Peltier, 2004).As for the run devoted to the mid-Holocene, LOCH and AGISM were not activated.The simulation was started from pre-industrial conditions.After 4000 years, the climate reached a quasi-equilibrium state characterized by a huge cooling of more than 25 • C over the Laurentide and Fennoscandian ice sheets (Fig. 21).The model also simulates a large cooling in the Southern Ocean associated with a large increase in the sea-ice extent there.The cooling is larger over the Atlantic than over the Pacific, in particular northward of 45 • N. In the tropics, the signal is weaker.In some regions, such as North Australia, the changes are very close to zero.Those results are similar to the ones of other simulations performed in the framework of the PMIP2 project (Braconnot et al., 2007), except in the Southern Ocean where the signal obtained in LOVECLIM is larger than the one given by most other models.
In the North Atlantic, the simulated cooling is associated with a southward shift of the sea-ice edge, with sea ice covering the majority of the Greenland, Iceland and Norwegian Seas both in summer and winter.Only a small area off the southern coast of Norway remains ice free all year long.In winter, deep convection occurs close to this location as well as south-east of Iceland.The North Atlantic meridional overturning streamfunction is quite similar to the one simulated for present-day conditions (Fig. 22), with a small decrease of the magnitude compared to present-day nearly everywhere except between 40 • and 60 • N in the top 2000 m of the water column.Furthermore, at high latitudes, the maximum is shifted southward, consistently with the change in the location of the convection patterns.Actually, the maximum of the overturning at LGM is lower here than in the previous versions of LOVECLIM that were characterized by a deeper and stronger meridional overturning at the LGM (e.g.Roche et al., 2007), a feature that previous versions of LOVECLIM shared with many of the other models participating in the PMIP2 intercomparison, although it is generally accepted that the circulation associated with North Atlantic Deep Water was shallower at LGM than at present (Weber et al., 2007;Lynch-Stieglitz et al., 2007).In Fig. 21, we also notice a reduction in the inflow of Antarctic Bottom Water in the Atlantic.At the global scale, the simulated deep circulation appears to be particularly weak in the Pacific and Indian ocean at the LGM and the magnitude of the deep cell close to Antarctica is reduced compared to present-day.

Summary and conclusions
In the previous sections, we have summarized the main equations and parameterizations of all the components of LOVE-CLIM.Furthermore, we have documented the model behaviour for present-day conditions and classical model tests.This provides a general overview and a reference for model users as well as for the scientists who want to learn more about the model, for instance after reading a paper using LOVECLIM results.A brief discussion of model performance is provided for several standard cases.As mentioned in the introduction, a deeper analysis was performed using previous versions of the model for all the experiments presented here.Further analysis is planned for the near future, for instance in the framework of PMIP3 (http://pmip3.lsce.ipsl.fr/).
The discussion of model results underlines that the model appears well adapted to study long-term climate changes, in particular at mid-and high-latitudes.However, we recall that it is of course essential to always try to take into account the model limitations and to estimate how they influence the conclusions of a study.Where the biases are strong, in particular in many regions at low latitudes, this requires a particularly careful analysis.In addition to simulations over long periods, the model is suitable and thus more and more used to perform studies that require large ensembles of simulations.This has not been discussed here but recent examples show, for instance, the influence of the choice of parameters in all the components of the model (Loutre et al., 2010;Goetzler et al., 2010) and the way data assimilation in coupled mode could help in reconstructing past climate changes (Crespin et al., 2009;Goosse et al., 2010).

Fig. 1 .
Fig. 1.Sketch of the LOVECLIM model showing the interactions between the five components.

Fig. 1 .
Fig. 1.Sketch of the LOVECLIM model showing the interactions between the five components.

Fig. 2 .
Fig. 2. Vertical discretization of the atmospheric model ECBilt.is the streamfunction, q is the quasi-geostrophic potential vorticity, T is the air temperature, T s is the surface temperature and q a is the total water content.

Fig. 3 .
Fig.3.Location of the various variables on the grid of of the barotropic velocity, η the surface elevation, u, v velocity, S the salinity, θ the potential temperature, q energy, and K s and K u the vertical diffusion and vertical

Fig. 3 .
Fig. 3.Location of the various variables on the grid of CLIO.U , V are the two components of the barotropic velocity, η the surface elevation, u, v, and w the three components of the velocity, S the salinity, θ the potential temperature, q 2 (two times) the turbulent kinetic energy, and Ks and K u the vertical diffusion and vertical viscosity.

Fig. 4 .
Fig. 4. The horizontal grid of CLIO at a resolution of 3° by 3°.The view is centred on the Atlantic.The two spherical subgrids in two different colors are connected in the Atlantic at the "geographical equator".

Fig. 4 .
Fig. 4. The horizontal grid of CLIO at a resolution of 3 • by 3 • .The view is centred on the Atlantic.The two spherical subgrids in two different colors are connected in the Atlantic at the "geographical equator".
www.geosci-model-dev.net/3/603/2010/Geosci.Model Dev., 3, 603-633, 2010 610 H. Goosse et al.: Description of LOVECLIM version 1.2 Fig.5Schematic representation of the main processes described in the LOCH model(Mouchet, 2010).The left panel focuses on purely biological processes, while the right panel shows the processes affecting the ocean carbon cycle.Up and down blue arrows represent transport processes (advection, diffusion, etc).Transported variables include dissolved inorganic carbon (DIC), alkalinity (Alk), dissolved inorganic phosphorus (DIP), dissolved organic matter (DOM), oxygen (O 2 ) and silica (Si).At the air-sea interface CO 2 and O 2 are exchanged with the atmosphere.B1 stands for opal building phytoplankton biomass and B2 represents the biomass of phytoplankton not relying on silica for growth (please note the inversion of B1 and B2 boxes between panels).POM decays at depth either as DOM or DIP.The flux of POM is governed following a power law function of the depth.Opal dissolves while sinking to the bottom.Calcareous shells (CaCO 3 ) reach the deepest layer where chemical conditions drive their dissolution or preservation.Fluxes to sediments, where permanent preservation prevails in this version, are also represented.Rivers (not illustrated) carry Si, DOM, DIC and Alk to the ocean.

Fig. 5 .
Fig. 5. Schematic representation of the main processes described in the LOCH model(Mouchet, 2010).The left panel focuses on purely biological processes, while the right panel shows the processes affecting the ocean carbon cycle.Up and down blue arrows represent transport processes (advection, diffusion, etc).Transported variables include dissolved inorganic carbon (DIC), alkalinity (Alk), dissolved inorganic phosphorus (DIP), dissolved organic matter (DOM), oxygen (O 2 ) and silica (Si).At the air-sea interface CO 2 and O 2 are exchanged with the atmosphere.B1 stands for opal building phytoplankton biomass and B2 represents the biomass of phytoplankton not relying on silica for growth (please note the inversion of B1 and B2 boxes between panels).POM decays at depth either as DOM or DIP.The flux of POM is governed following a power law function of the depth.Opal dissolves while sinking to the bottom.Calcareous shells (CaCO 3 ) reach the deepest layer where chemical conditions drive their dissolution or preservation.Fluxes to sediments, where permanent preservation prevails in this version, are also represented.Rivers (not illustrated) carry Si, DOM, DIC and Alk to the ocean.
H.Goosse et al.: Description of LOVECLIM version 1.2 Fig. 6.Structure of the three-dimensional ice sheet model AGISM.The inputs are given at the left-hand side.Prescribed environmental variables drive the model, which has ice shelves, grounded ice and bed adjustment as major components.For the Antarctic component, the position of the grounding line follows from a flotation criterion and a specific treatment of the force balance.Ice thickness feeds back on surface elevation, an important parameter for the calculation of the mass balance.The main model outputs the time-dependent ice sheet geometry and the coupled temperature and velocity fields.

Fig. 6 .
Fig. 6.Structure of the 3-D ice sheet model AGISM.The inputs are given at the left-hand side.Prescribed environmental variables drive the model, which has ice shelves, grounded ice and bed adjustment as major components.For the Antarctic component, the position of the grounding line follows from a flotation criterion and a specific treatment of the force balance.Ice thickness feeds back on surface elevation, an important parameter for the calculation of the mass balance.The main model outputs the time-dependent ice sheet geometry and the coupled temperature and velocity fields.

Fig. 7 .
Fig. 7.The numerical grid of AGISM has a horizontal resolution of 10 km for both polar ice sheets (left panel: Antarctic ice sheet; right panel: Greenland ice sheet).Major gridlines are for a distance of 100 km, the insets show the detailed meshes employed in the calculations.The numbers along the axes are gridpoint numbers (561 × 561 gridpoints for AISM, 165 × 281 for GISM).For illustrative purposes, we display observed surface elevation of the present-day ice sheets.Ice sheet cover is shaded grey, ice-free range from green to white, and blue colours depict the ocean.Contour lines over the ice sheets are for every 250 m of elevation, major ones for every 1000 m are shown in thick.
Fig. 8. Fraction of ocean surface in each of the grid cell of ECBilt.

Fig. 8 .
Fig. 8. Fraction of ocean surface in each of the grid cell of ECBilt.
Fig. 12 Location of the sea-ice edge averaged over the period 1980-2000, defined by a monthly ice concentration equal to 15% in (a) March in the Northern Hemisphere, (b) September in the Northern Hemisphere, (c) September in the Southern Hemisphere, (d) March in the Southern Hemisphere.The observations are in blue (Rayner et al. 2003) and LOVECLIM1.2results are in red.

Fig. 12 .
Fig. 12. Location of the sea-ice edge averaged over the period 1980-2000, defined by a monthly ice concentration equal to 15% in (a) March in the Northern Hemisphere, (b) September in the Northern Hemisphere, (c) September in the Southern Hemisphere, (d) March in the Southern Hemisphere.The observations are in blue(Rayner et al., 2003) and LOVECLIM1.2results are in red.
Fig. 13 Simulated meridional overturning streamfunction averaged over the period 1980-2000 (in Sv) for (a) the whole World Ocean and (b) the Atlantic.
Fig. 17.Continental ice volume changes during the last millennium simulated by AGI for (a) Antarctic and (b) Greenland ice sheets.In this particular example, the Greenland volume budget is equivalent to a positive sea level contribution of about 10 cm over entire period.The Antarctic ice volume budget is slightly positive but cannot be dire related to sea level change because of ice grounded below sea level.Variability in b indices on centennial timescales arises from the climate forcing and dynamical ice-clim interactions.The modelled trend is not a robust feature of AGISM, but contain significant component from the model coupling procedure at 500 AD and the spec model parameters selected for ECBilt and CLIO.

Fig. 17 .
Fig. 17.Continental ice volume changes during the last millennium simulated by AGISM for (a) Antarctic and (b) Greenland ice sheets.In this particular example, the Greenland ice volume budget is equivalent to a positive sea level contribution of about 10 cm over the entire period.The Antarctic ice volume budget is slightly positive but cannot be directly related to sea level change because of ice grounded below sea level.Variability in both indices on centennial timescales arises from the climate forcing and dynamical ice-climate interactions.The modelled trend is not a robust feature of AGISM, but contains a significant component from the model coupling procedure at 500 AD and the specific model parameters selected for ECBilt and CLIO.
Fig. 18 Simulated difference of summer (JJAS) temperature (in °C) between the mid-Holocene and present-day conditions.
Fig. 19.Simulated difference of summer (JJAS) precipitation (in mm per day) between the mid-Holocene and present-day conditions.

Fig. 19 .
Fig. 19.Simulated difference of summer (JJAS) precipitation (in mm per day) between the mid-Holocene and present-day conditions.

Fig. 20 .
Fig. 20.Simulated difference in total vegetation cover (%) between the mid-Holocene and present-day conditions.

Fig. 20 .
Fig. 20.Simulated difference in total vegetation cover (%) between the mid-Holocene and present-day conditions.

Fig. 21 .
Fig. 21.Simulated difference of annual mean surface temperatures (in °C) between the last glacial maximum and present-day conditions.

Fig. 21 .
Fig. 21.Simulated difference of annual mean surface temperatures (in • C) between the last glacial maximum and present-day conditions.95

Fig. 22 .
Fig. 22. Simulated meridional overturning streamfunction (in Sv) for (a) the whole World Ocean and (b) the Atlantic simulated for the Last Glacial Maximum.

Table 2 .
Major parameters of CLIO.

Table 3 .
Major parameters of VECODE.Bioclimatic parameter related to minimum precipitation required for vegetation P 0

Table 4 .
Major parameters of LOCH.

Table 5
) corresponds to a lithospheric thickness of 115 km.The steady state deflection of the surface of the Earth is used to calculate the degree to which the Earth is in isostatic equilibrium, www.geosci-model-dev.net/3/603/2010/Geosci.Model Dev., 3, 603-633, 2010

Table 6 .
Major parameters of the iceberg model.

Table 7 .
Simulated trends over the last decades of some important variables.Data from the Mauna Loa record (NOAA ESRL; www.esrl.noaa.gov/gmd/ccgg/trends/).