Model Development Development of a system emulating the global carbon cycle in Earth system models

Abstract. Recent studies have indicated that the uncertainty in the global carbon cycle may have a significant impact on the climate. Since state of the art models are too computationally expensive for it to be possible to explore their parametric uncertainty in anything approaching a comprehensive fashion, we have developed a simplified system for investigating this problem. By combining the strong points of general circulation models (GCMs), which contain detailed and complex processes, and Earth system models of intermediate complexity (EMICs), which are quick and capable of large ensembles, we have developed a loosely coupled model (LCM) which can represent the outputs of a GCM-based Earth system model, using much smaller computational resources. We address the problem of relatively poor representation of precipitation within our EMIC, which prevents us from directly coupling it to a vegetation model, by coupling it to a precomputed transient simulation using a full GCM. The LCM consists of three components: an EMIC (MIROC-lite) which consists of a 2-D energy balance atmosphere coupled to a low resolution 3-D GCM ocean (COCO) including an ocean carbon cycle (an NPZD-type marine ecosystem model); a state of the art vegetation model (Sim-CYCLE); and a database of daily temperature, precipitation, and other necessary climatic fields to drive Sim-CYCLE from a precomputed transient simulation from a state of the art AOGCM. The transient warming of the climate system is calculated from MIROC-lite, with the global temperature anomaly used to select the most appropriate annual climatic field from the pre-computed AOGCM simulation which, in this case, is a 1% pa increasing CO2 concentration scenario. By adjusting the effective climate sensitivity (equivalent to the equilibrium climate sensitivity for an energy balance model) of MIROC-lite, the transient warming of the LCM could be adjusted to closely follow the low sensitivity (with an equilibrium climate sensitivity of 4.0 K) version of MIROC3.2. By tuning of the physical and biogeochemical parameters it was possible to reasonably reproduce the bulk physical and biogeochemical properties of previously published CO2 stabilisation scenarios for that model. As an example of an application of the LCM, the behavior of the high sensitivity version of MIROC3.2 (with a 6.3 K equilibrium climate sensitivity) is also demonstrated. Given the highly adjustable nature of the model, we believe that the LCM should be a very useful tool for studying uncertainty in global climate change, and we have named the model, JUMP-LCM, after the name of our research group (Japan Uncertainty Modelling Project).

By adjusting the effective climate sensitivity (equivalent to the equilibrium climate sensitivity for an energy balance model) of MIROC-lite, the transient warming of the LCM could be adjusted to closely follow the low sensitivity (with an equilibrium climate sensitivity of 4.0 K) version of MIROC3.2.By tuning of the physical and biogeochemical parameters it was possible to reasonably reproduce the bulk physical and biogeochemical properties of previously published CO 2 stabilisation scenarios for that model.As an example of an application of the LCM, the behavior of the high sensitivity version of MIROC3.2 (with a 6.3 K equilibrium climate sensitivity) is also demonstrated.Given the highly adjustable nature of the model, we believe that the LCM should be a very useful tool for studying uncertainty in global climate change, and we have named the model, JUMP-LCM, after the name of our research group (Japan Uncertainty Modelling Project).

Introduction
It is now increasingly common for climate models used for projections of climate change to explicitly include representation of the carbon cycle.While atmosphere-only general circulation models were called AGCMs, and those with coupled oceans termed AOGCMs, models with more coupled components, which may include various different elements such as ice sheets, atmospheric chemistry and the carbon cycle are increasingly called Earth system models (ESMs), and this is the nomenclature we adopt here.
The inclusion of a carbon cycle gives rise to additional sources of uncertainty, on top of those in the physical system, relating to feedbacks in the carbon cycle.The contribution of carbon cycle uncertainty to the uncertainty in the transient Published by Copernicus Publications on behalf of the European Geosciences Union.
climate response has been estimated, by Huntingford et al. (2009) using box models to emulate C4MIP ESMs, to be around 40% of that of the uncertainty in equilibrium climate sensitivity (the equilibrium temperature response to doubled CO 2 ) and heat capacity.Such uncertainties may have substantial implications for mitigation and adaptation policies relating to climate change.Thus, even as the models increase in complexity and therefore computational cost, it is more important than ever before to be able to perform ensemble experiments in order to investigate uncertainties in the physical and biogeochemical processes, and thus in the climate change projections themselves.
Of course, large ensembles of the most costly models (which are generally designed so as to be capable of running only a handful of simulations on current hardware) are not computationally feasible.Therefore, we inevitably have to simplify the model in some way, and a wide range of so-called Earth System Models of Intermediate Complexity (EMICs) have been developed (Claussen et al., 2002;Weber, 2010).
The main distinguishing feature of such models is a reduction in resolution and/or complexity of some model components, resulting in a substantial reduction in computational cost.One common approach is to substitute an energymoisture balance model (EMBM) in place of a fully dynamical atmospheric GCM (e.g.UVic, Weaver et al., 2001;Bern, Plattner et al., 2001;GENIE, Edwards and Marsh, 2005).Such a model can reasonably represent the behaviour of more complex GCMs, at least for large-scale physical variables such as globally-averaged surface air temperature on multidecadal to centennial scales (Raper et al., 2001).
A typical limitation of such an EMBM, however, is the inability to represent spatial details of the current climate, such as patterns of precipitation and cloud cover.This is a particular problem when we include sophisticated representations of the carbon cycle, since precipitation and radiation are essential factors for plant life.We wish to ensure that our reduced-complexity model is as traceable as possible to the full GCM, and therefore we would like to use a detailed terrestrial carbon cycle model such as Sim-CYCLE (Ito and Oikawa, 2002) which forms part of the MIROC3.2ESM (Kawamiya et al., 2005).Thus, we require some way of efficiently reproducing the detailed physical output of GCM in a more efficient EMIC.
Pattern scaling has been proposed as one method for projecting time-varying climate changes of a GCM in a computationally efficient manner (Santer et al., 1990;Mitchell et al., 1999).In this approach, the spatial pattern of climate change anomalies is assumed fixed, and calculated as the difference between a control run (e.g. a pre-industrial climate simulation) and an equilibrium run under different boundary conditions, typically 2×CO 2 with a slab ocean model.For transient simulations, the pattern of climate change is then scaled by the global mean temperature, which can be calculated using a simpler EMIC, or even derived from an energy balance model.Of course, the validity of this approach depends both on the pattern of climate change being constant in time and on it being well represented by the equilibrium integration.While these are reasonable first-order approximations, they introduce a source of error, and therefore additional uncertainty, into the system (Mitchell et al., 1999).A more sophisticated approach is that of Hooss et al. (2001), who use empirical orthogonal functions to make what they call an impulse response model.All these approaches generate a deterministic pattern of change which will not include natural variability.
In this paper, we present an alternative approach, which introduces negligible additional computational cost in comparison to pattern scaling, but which can, in principle, fully represent both the spatial changes and the natural variability of a transient climate change simulation.The innovation is that rather than using a single climate change pattern derived from an equilibrium simulation, we use the transient output from a previous transient simulation such as the 1% pa 4×CO 2 runs of the CMIP3 project (Meehl et al., 2007).For a given global mean surface temperature anomaly (provided by the EMIC), the year in the transient run that best approximates this temperature anomaly is selected, and the year of climate model data are then used to force the state of the art terrestrial carbon cycle model.If the trajectory of CO 2 mixing ratio of the LCM simulation matches that of the transient simulation of the state of the art model, the EMIC-based results should accurately mimic the full ESM at a small fraction of the computational cost.For reasonable deviations in the trajectory of the mixing ratio (as might arise through changes in model parameters or emissions scenario), the pattern of climate change from the transient run should still be more accurate than that provided by the scaled equilibrium pattern.We illustrate the approach by emulating two versions of the MIROC3.2ESM.We introduce the models and coupling methodology in Sect. 2. In Sect. 3 we describe the tuning of the LCM to the lower sensitivity version of MIROC3.2.The discussion and conclusions follow in Sects.4 and 5.

Basic structure of the loosely coupled model (LCM)
In the LCM system we have developed, three stand-alone components, or sub-models are loosely coupled by unix shell scripting.The three components are: MIROC-lite, a simplified energy balance atmosphere version of MIROC3.2AOGCM with coupled ocean carbon cycle; Sim-CYCLE, the terrestrial carbon cycle component of the MIROC3.2ESM; and output from runs made with MIROC3.2AOGCM.MIROC3.2 was one of the models which contributed to the IPCC AR4 (IPCC, 2007).See Fig. 1 for the basic design.The design makes it relatively straightforward to replace the components of the LCM for new model versions, or even for different models.The three models that comprise the components of the LCM have been described in detail elsewhere.In this section we introduce them briefly, and then describe the coupling methodology in more detail.

MIROC-lite: an EMIC
MIROC-lite (Oka et al., 2001) is a simplified version of MIROC3.2 in which the OGCM component is the same (albeit at low resolution) but the AGCM of MIROC3.2 is replaced by a 2-D energy moisture balance model atmosphere.
The model diagnoses the surface air temperature by solving the vertically integrated energy balance equation below.
Where C a (in J/(m 2 K)) is the heat capacity of the air column, T is the surface air temperature (in Celsius), t is time, Q sw and Q lw (W/m 2 ) are the net incoming shortwave and outgoing longwave radiation at the top of atmosphere, Q t (W/m 2 ) is the convergence of the horizontal heat transport by the atmosphere, and F T (W/m 2 ) is the net downward heat flux at sea/land surface.
The surface wind and the freshwater flux are diagnosed from the distribution of the surface air temperature and the convergence of atmospheric water transport, respectively.The meridional wind is decomposed into the Hadley and Ferrel circulations which are considered separately.Both types of circulation are described as proportional to the North-South temperature gradient, and the latitude dependent empirical coefficients are determined by a mother GCM (MIROC3.2)'sresult.In the current setting, rain and snowmelt on land is returned to the nearest ocean grid, after taking account of evaporation (interception by vegetation is not considered here).
The ocean component of the model is COCO (Hasumi and Suginohara, 1999), an ocean GCM.The version we use in this study includes a sea-ice component.In order to use a long time-step, the acceleration method of Bryan (1984) is used in MIROC-lite.
Unlike the original version, the latest version of MIROClite (Oka et al., 2010) includes the feedback of the marine carbon cycle.There are two options: (1) a simple carbon cycle with no marine ecosystem model (Yamanaka and Tajika, 1996) which considers nutrient, DIC and alkalinity (in addition to the physical ocean's temperature and salinity), or (2) a carbon cycle model considering marine ecosystem (Palmer and Totterdell, 2001) which includes the above 5 tracers plus phytoplankton, zooplankton and detritus making 8 in total.In this manuscript, ( 2) is used as it is closer to the ocean carbon cycle used in the MIROC3.2ESM.
Here we use the latest version (Oka et al., 2010) of MIROC-lite with some adjustments.First, we impose a freshwater flux adjustment to compensate for the poor representation of the freshwater flux from the Atlantic to the Pacific.Following the traditional method for EMICs with 2-D EMBM atmosphere, we set an artificial freshwater flux (FWF) adjustment.Oort (1983) stated 0.32 Sv in total: with 0.18, 0.17, −0.03 Sv for the bands north of 24 N, 20 S to 24 N, and south of 20 S, respectively.The model's own internally-generated flux is negligible, and we use the Oort's values for the FWF adjustment.
A second modification is to modify the outgoing longwave parameterisation in order to both account for forcing through atmospheric CO 2 concentration (following Table 6.2 of IPCC's TAR), and also to allow the climate sensitivity to be varied in the same way as Plattner et al. (2001): where A, and B are the constants (206.8778 and 1.73357); T is the surface air temperature (in Celsius) of the grid concerned; pCO 2 is the atmospheric CO 2 concentration (in ppm); gT is the global average of T ; and gT c is gT for 1×CO 2 (both are in Celsius).RF 2× (W/m 2 ) in the third term is the 2×CO 2 radiative forcing of the GCM concerned.
Figure 2 shows that by varying the coefficient C in Eq. ( 2), we can obtain a range in climate sensitivity large enough to cover the full range found in GCMs.Here, the value 3.09 is RF 2× for the lower sensitivity version (see Sect. 2.4) of MIROC3.2 (MIROC3.2-LS).A curve of the form climate sensitivity = 3.09/(C +X) is fit to the numerical results (black Fig. 2. Climate sensitivity adjustment.For the red curve, 3.09 (W/m 2 ) is the radiative forcing for 2×CO 2 and C+1.61 is the total (equilibrium or effective) climate feedback parameter.symbols).The resulting curve, with X = 1.61, is shown as a red line.This fitted value of 1.61 is within the ranges of the equilibrium and the effective climate feedback parameters for recent AOGCMs mentioned by Yokohata et al. (2008).It should be noted that an EMBM has a same value for the equilibrium and the effective climate sensitivity, while as discussed by Yokohata et al. (2008) generally this is not the case for GCMs (mainly because of difference in feedbacks, although there may be some effect from using a slab ocean to estimate the equilibrium climate sensitivity).For MIROC3.2-LS, for example, the equilibrium climate sensitivity is 4.0 K, while the effective climate sensitivity is 4.7 K.
See the Appendix A for further details of the configuration of MIROC-lite, including the time step, spatial resolution and topography.

Sim-CYCLE
Sim-CYCLE (Ito and Oikawa, 2002) is a process based terrestrial carbon cycle model, which was developed based on an ecosystem scale model by Oikawa (1985).The origin of these models is in the dry-matter production theory proposed by Monsi and Saeki (1953).
In this model, ecosystem carbon storage is divided into plant biomass and soil organic carbon, which are further subdivided into five compartments: foliage, stem, and root for plant biomass; litter and mineral soil for soil organic matter.The model also has a water and radiation process, as carbon dynamics is closely coupled with these processes.The single-leaf photosynthetic rate (PC) is formulated as a Michaelis-type function of the incident photosynthetic photon flux density (PPFD in ): where PC sat is PC for the light saturation condition; QE is light-use efficiency.PC sat and QE are formulated as (maximum value) × (stress function), where as stresses, temperature, CO 2 level, air humidity and soil water (the parameters are different for C3/C4/crop plants) are taken into consideration.
Ecosystem scale gross primary production (GPP) is calculated under an assumption of exponential attenuation of PAR irradiance due to leaves' mutual shading.
Autotrophic plant respiration consists of two components: the maintenance respiration, and the growth respiration.The amount of the maintenance respiration per unit existing carbon is an exponential function of temperature (degree Celsius) with a coefficient of so called Q10, while the growth respiration is proportional to a net biomass gain, where ARM and T are the maintenance respiration per unit biomass, and the temperature.Soil organic carbon is divided into two components: the labile part of litter which circulates once in several months or a year, and the passive part in mineral soil which remains for decades or centuries.Heterotrophic soil respiration is composed of two components for these two.For both, temperature and soil moisture conditions are influential.For temperature dependence, an Arrhenius type function proposed by Lloyd and Taylor (1994) is used.
The fraction of C4, the fraction of crop plants, and the distribution of 19 biomes (based on the classification of Matthews, 1984) are pre-determined.Thus this is not a dynamic vegetation model.The parameters are determined using observational data of 21 sites for a variety of vegetation types.
Sim-CYCLE has daily time steps and thus needs daily input climatic data (air temperature at 2 m height, precipitation, ground surface temperature, soil temperature at 10 cm depth, soil temperature at 200 cm depth, specific humidity, wind speed, and ground surface radiation).The model can be used in both equilibrium and transient mode.

MIROC3.2: description of the model and the dataset used for the LCM
The third component of the LCM is the archive of climate fields from a preperformed run of a GCM.Here we use output from the medium resolution (T42) version of MIROC3.2, a Japanese AOGCM, which includes five physical components: atmosphere, land, river, sea ice, and ocean (K-1 model developers, 2004).The atmospheric model has 20 vertical σ -layers.The model has an interactive aerosol module, simplified SPRINTARS (Takemura et al., 2000(Takemura et al., , 2002)), and a land surface model MATSIRO (Takata et al., 2003).The ocean component is the same as in MIROC-lite, COCO (Hasumi and Suginohara, 1999).However, the resolution here is higher at 256×192 (×44 layers).As for the marine carbon cycle model, a model based on Oschlies and Garcon (1999) and Oschlies (2001) is used.
MIROC3.2 has two sensitivity versions, one with a equilibrium climate sensitivity of 4.0 K (lower sensitivity, LS) and one with a sensitivity of 6.3 K (higher sensitivity, HS).These only differ in the cloud treatment, and both of them provide realistic simulations of the mean present climate (Ogura et al., 2008).The difference between the two versions is in the treatment of cloud microphysics.According to Ogura et al. (2008), there are three differences: (1) mixed phase (i.e., solid and liquid) temperature range, (2) form of melted cloud ice, and (3) values of two parameters included in formulations in precipitation rate and sedimentation of cloud particles.The effect of these differences, is to change the sign of response in cloud condensate to the doubled CO 2 concentration (positive for LS and negative for HS).Yokohata et al. (2005) compared their response to the Pinatubo volcanic forcing and conclude that LS provided more realistic response, while HS's response is too strong.
The output of the standard 1% pa compound CO 2 enrichment experiment from MIROC3.2-LS prepared as part of the CMIP3 experiments (Meehl et al., 2007) for the last IPCC report (AR4) was used as the dataset in this implementation of the LCM.The increment is started from the pre-industrial state.We use one of the three ensemble members.

Coupling method
The coupling process (see Fig. 1) is organised as follows.
(1) MIROC-lite runs one year with a given CO 2 concentration.This may be either directly prescribed (in the case of a concentration scenario), or in the case of an emissions scenario, diagnosed from the previous year's concentration updated with the annual emissions and the feedback from the carbon cycle from previous year; from (3).As the marine carbon cycle is also switched on, the amount of CO 2 which is absorbed (or released) by the ocean is also calculated and stored.
(2) The value of the global mean surface air temperature in MIROC-lite is used to pick a year from the preperformed MIROC3.2GCM 1% per year increasing CO 2 run (described in the previous section).In order to choose the year, a quadratic curve is used to smooth the annual mean temperature time series from the GCM data set.Thus the actual year chosen will have a mean temperature that differs from the smoothed value due to interannual variability.For the appropriate year, daily climatic fields are extracted from the archive.
(3) The year's worth of daily climatic data from (2) and the CO 2 level determined in (1) are used to drive Sim-CYCLE for one year.We calculate the change in the total terrestrial ecosystem carbon storage for evaluating the feedback of the terrestrial carbon cycle (note: the albedo and evapotranspiration feedback is not considered here).In emission scenario experiments, the total feedback of carbon budget is evaluated by calculating the sum of the terrestrial and marine carbon uptake; obtained in (1).Then, return to (1).
On our supercomputer (SGI Altix 4700), the system runs one model year in around 1.3 to 1.4 min on a single processor.Thus century-length ensemble integrations are easily achievable.

Testing and tuning the LCM
Using the results of FWF adjustment, we obtained acceptable climatology in atmospheric temperature, sea surface temperature, sea surface salinity and North Atlantic meridional overturning circulation (see Figs. 3 and 4 for the model after parameter tuning), although precipitation (particularly for land) is still not adequate to be coupled with a terrestrial vegetation model.
The changes in the annual mean surface air temperature for the 1% incremental run of MIROC3.2-LS/HS are presented as blue/red lines in Fig. 5.
In order to check the performance for a transient run, we reproduced the experiments of Miyama and Kawamiya (2009).
These experiments use the state of the art MIROC3.2-LSESM with oceanic and terrestrial carbon cycle, forced by the stabilisation scenarios of Mueller (2004) (represented as Fig. A2).As our main focus is on transient runs, climate sensitivity tuning was carried out to the effective climate sensitivity.After fixing the climate sensitivity by choosing the appropriate value of C, for the MIROC3.2-LS'seffective climate feedback parameter (0.66, Yokohata et al., 2008), in the previous equation (Eq.2), we can reproduce the trend of MIROC3.2-LS'sbehaviour in the transient run (cyan line in Fig. 5).
Although the default parameter set (with adjusted climate sensitivity as described above) provides a good simulation of the physical transient, it also resulted in too much carbon uptake by the ocean.Therefore, we performed some ensemble experiments with perturbed parameters to investigate and tune the response of the model.For the physical parameters, we considered those which have a strong influence on mixing and circulation in the ocean (i.e., vertical diffusivity, horizontal diffusivity, and GM thickness diffusion Gent and McWilliams, 1990), as these should also influence the ocean's carbon uptake.
www.geosci-model-dev.net/3/365/2010/Geosci.Model Dev., 3, 365-376, 2010  In MIROC-lite, however, there is another very important parameter to determine the ocean's carbon uptake.In MIROC-lite, as in MIROC3.2, the air-sea CO 2 exchange (F , in mol/(sm 2 )) is formulated as: where pCO 2 a and pCO 2 o are partial pressure (in ppmv) of CO 2 in the atmosphere and the ocean surface, the CO 2 gas transfer velocity CK (in m/s) = a • u 2 / √ Sc/660 with Sc (a function of SST, dimensionless) being the Schmidt number, and Sl is solubility (in mol/(ppmv m 3 ), depending on T , S).For the detail of Eq. ( 5), including the formulation of Sc and Sl, see the OCMIP protocol (Orr et al., 2000).
Unlike MIROC3.2,however, in MIROC-lite the wind speed u (m/s) in this equation is fixed as a globally and temporally constant value, and this parameter has a large influence in determining the amount of carbon uptake by the ocean.
Thus, in total four parameters were perturbed in the ensemble experiment to select the best-fit parameter set to represent the GCM's oceanic CO 2 uptake.Figure 6a presents the dominant contribution of the wind speed to the maximum carbon uptake by the ocean.Particularly up to 4 m/s, the maximum marine carbon uptake is controlled mostly by the wind speed alone.On the other hand, as depicted in Fig. 6b, the vertical diffusivity has major effect on the Atlantic overturning circulation, but the contribution of this parameter to the marine carbon uptake is to a limited degree in comparison to that of the wind speed.
The output (an equilibrium state after a 3000 year spinup) of some of the variables using the best parameter set is presented as Fig. 3a-d.The model has acceptable performance in latitudinal change (i.e., the zonal mean is well-reproduced), but the longitudinal change including the land-ocean contrast is not so well reproduced.
When looking at the derivation from the reference (observation or reanalysis) data (Fig. 3e-h), the most obvious problem for the basic MIROC-lite model is, as mentioned in Sect. 1, the precipitation, which does not penetrate into the continental interiors.However, this is not used to drive the Sim-CYCLE terrestrial ecosystem model, and thus has little impact on our simulations.As for the Atlantic meridional overturning, the stream function (Fig. 4) is acceptable and the maximum value (20.7 Sv) is close to that of MIROC3.2-LS(20.9 Sv).In order to check how well the final version LCM (JUMP-LCM) can reproduce the behavior of the ESM using MIROC3.2-LS,air temperature, sea surface temperature, terrestrial/marine carbon uptakes are examined for stabilisation scenarios of 450, 550, and 1000 ppm (Fig. 7).As shown in the figure except for the short-term variabilities, the JUMP-LCM reproduced the basic shapes of the curves as well as the magnitude of the peak values.The only noticeable difference is observed in the marine carbon uptake for the 1000 ppm scenario (Figs.7g and 7h), indicating the limitation of a model tuned for 450 ppm scenario.
For comparison we also performed a similar experiment to mimic the results from MIROC3.2-HS.For this experiment we can only compare the physical outputs as there are www.geosci-model-dev.net/3/365/2010/Geosci.Model Dev., 3, 365-376, 2010 no results from a full Earth system model based on this physical model.Thus, we reproduce a 1% pa CO 2 enrichment scenario.As the effective climate feedback parameter (λ eff ) for MIROC3.2-HShas not been reported, we tuned this using the result of 1% pa experiment.From the estimated radiative forcing and the temperature change in the experiment, we obtained λ eff as 0.47.The transient temperature change by the tuned MIROC-lite-HS for the 1% pa CO 2 enrichment experiment is shown as the magenta line in Fig. 5.Although we did not tune ocean physical parameters relating to the transient response (as the ocean component is common for MIROC3.2 and MIROC-lite), tuning ocean heat uptake efficiency may also be needed when the JUMP-LCM is used for mimicking other ESMs.
In order to show the potential of JUMP-LCM to mimic other models, the results from the stabilisation scenario experiments for HS version are shown in Fig. 8.The change we have made to represent the HS version is: (1) to change value of C in Eq. ( 2) so that the resultant transient climate response follows that of MIROC3.2-HS and (2) to replace the data archive ((2) in Fig. 1) with that of MIROC3.2-HS.As the radiative forcing of HS was reported as almost identical to that of LS (the difference is less than 0.1 W/m 2 (Yokohata et al., 2005)), we used the same value for HS too.For this experiment we do not have validation data, such as Miyama and Kawamiya (2009) for the LS version.The results suggest a lower carbon uptake by the land, presumably due to stronger warming.The same tendency is also seen in the ocean too, but to a more moderate degree.
It can be seen from Fig. 7 that the natural variability evident in the GCM output is not reproduced in JUMP-LCM apart for the land carbon uptake, which was driven by the GCM's climatic field.We can attempt to represent natural variability in the physical system, by adding a random number term to the radiative forcing calculation.Natural variability is thought to result from radiative forcing and interactions between different components of the climate system.Pelletier (1997) showed that natural variability has a certain power spectrum and Hoerling et al. (2008) stated that multidecadal variabilities are mainly controlled by external radiative forcing due to GHGs, aerosol, solar and volcanic variations.
On the other hand, Wigley and Raper (1990) demonstrated that because of the ocean's large heat capacity, a random white noise forcing results in a red spectrum in the global mean temperature.With these facts in mind, we concluded that as a start it is reasonable to add a random number term to Eq. ( 2).Here we tested 3.46•(RN[0,1]-0.5),where RN[0,1] is a random number (generated for each year and kept constant in a year) with a uniform distribution between 0 and 1.The coefficient 3.46 was determined so that the resultant standard deviation of the random term is the value (1 W/m 2 ) which Wigley and Raper (1990) mentioned as a suitable standard deviation in interannual radiative forcing.As shown in Fig. 9, by adding this term we could reproduce the natural variability and a good by-product is an improved land carbon uptake curve by 3.46•RN[0,1] (Fig. 9c).
Further investigation will be needed for this issue.

Discussion
The loosely coupled modelling system introduced by this manuscript reproduces the transient carbon cycle tions of a full state of the art Earth system model, at a fraction of the computational cost.Therefore, it should be a powerful tool for investigating uncertainty in climate change, using perturbed parameter ensembles.It is straightforward to change the climate sensitivity of the atmospheric model, and all internal parameters of the ocean GCM and carbon cycle (both terrestrial and ocean).The different spatial patterns of climate change arising from different climate models from around the world could in principle be utilised by simply swapping in the results from the CMIP3 database.Thus we believe that the loosely coupled system we present here can conveniently and efficiently account for all major sources of uncertainty in the climate's response to elevated CO 2 levels.
The limitation of the database may generate some problems.For example, for the coupled run with 1000 ppm scenario, the global temperature went out of range of the database at year 2279.While it would in principle be possible to extrapolate the database, this has not been implemented.For century length integrations, however, this is unlikely to be a problem for the 4×CO 2 database.
For long-term equilibration experiments, the climate field of the full system would approach that of an equilibrium experiment rather than the transient that we use (a 1% per year incremental increasing experiment).However, even in this case, the standard approach of pattern scaling from an equilibrium (slab ocean) run is not an ideal approach either, since this ignores the issue of ocean response.In any case, we expect our approach to be most accurate during a period of steadily increasing temperature, which probably covers most plausible scenarios at least over the next century.
Using transient reference data may be one reason why the LCM overestimates the ocean's carbon uptake after the peak, and thus to get a good performance for the total (accumulated) carbon uptake the peak value should be significantly smaller than the target value (about 2.5 Pg C/y).More discussion on this issue is in the Appendix A.
Currently, the feedback processes from vegetation to the atmosphere apart from change in the total carbon storage (e.g., change in albedo, evapotranspiration and sensible heat flux) are not considered.We should mention that these changes were not included even in the complex and costly MIROC ESM simulation that we emulate here.As well as considering other realistic river maps, this is also future work.However, it is expected that the effect of sum of these feedbacks will not change the results very greatly.

Conclusions
In order to utilize the strengths of both GCMs and EMICs, we developed a loosely coupled model (LCM) system connecting an EMIC, vegetation model and existing GCM output.We expect the result to be a powerful tool for studying the uncertainty in the carbon cycle and its contribution to the future climate change.The LCM reproduced the basic behaviour of the MIROC3.2ESM for transient runs very accurately over the 21st century, with a modest error over longer term equilibration scenarios.Using this system we intend, by varying model parameters, to investigate uncertainty, particularly in the carbon cycle components of MIROC3.2, and also to extend the approach to other versions of MIROC and other ESMs.

Configuration of MIROC-lite
In the present work, the spatial resolution of MIROC-lite is set as 6×6 degree (60×30 gridboxes, to cover the entire globe) and the ocean has 15 layers of unequal thickness (thinner at the surface) down to 5500 m depth.The thickness of Fig. A2.CO 2 stabilization scenarios used here (Mueller, 2004).the first (shallowest) ocean layer is 50 m, functioning as the mixed layer.The time step is 36 h and on a single CPU of our SGI Altix 4700, it takes around 32 (16) h for 3000 year integration with (without) marine ecosystem.The topography and the bathymetry of MIROC-lite are shown in Fig. A1.

A1 The stabilization scenarios
The stabilization scenarios used in the study is presented as Figure A2.

A2 Errors in representing the peak and the total marine carbon uptake
Figure A3 presents the relation between the errors in the peak and the total of the ocean's carbon uptake and shows that the distribution of the ensemble members (of the experiment Geosci.Model Dev., 3, 365-376, 2010 www.geosci-model-dev.net/3/365/2010/described in Sect.3) do not pass though the origin point (0, 0), but instead pass through (−0.2, 0) and (0, 0.2).So, the best performer in one variable is not the best in another.
To check whether this is due to the simplified wind speed treatment (described in Sect.3) we looked at the temporal change in wind speed and its variance for MIROC3.2-ESM'sexperiment (not shown) and found that the wind speed as well as its variance increased before the period of the peak marine carbon uptake but did not change significantly after that (variance is checked since the contribution of the wind speed to the carbon uptake is nonlinear).As the characteristics of the GCM's wind field do not change significantly after the period of the maximum marine carbon uptake, once the peak carbon uptake is tuned to the GCM's value, total carbon uptake should also be close to the GCM's value.Thus, it is not likely that using a constant wind speed is the main reason why good performance for the peak/total amount of the marine carbon uptake is not synchronized.

Fig. 3 .Fig. 4 .
Fig. 3. Output (a-d) and deviance from reference (e-h) for basic variables.After tuning in Sect.3.An equilibrium state after a 3000 year spin-up is presented.

Fig. 5 .
Fig. 5. Change in the annual mean surface air temperature for the 1% incremental run of MIROC3.2-LS/HS and MIROC-lite emulating the two versions of MIROC3.2 (CO 2 (ppm) used the right axis).

Fig. 7 .
Fig. 7. Comparison of LCM's and ESM's outputs for representative variables.Solid/broken lines are coupled/uncoupled runs, and red/green/black are 450/550/1000 ppm scenarios; in (a) and (c) red/green broken lines are hidden under the black broken line.

Fig. 8 .
Fig. 8.Result for HS version (based on the best fit parameters to the LS version).Lines are same as Fig.7.In (c), the black solid line after the year 2233 is not presented, as it became too warm to refer the GCM's archive.

Fig. 9 .
Fig. 9.Result of adding the natural variability term.Each of five types of black lines presented a random number run, while red/green/blue curves are GCM/LCM/ensemble mean of 5 random number runs.
Fig. A1.Land and ocean used in this manuscript.

Fig. A3 .
Fig. A3.Relation between the error in the maximum and the total carbon uptake.