MESMO 2: a mechanistic marine silica cycle and coupling to a simple terrestrial scheme

Here we describe the second version of Minnesota Earth System Model for Ocean biogeochemistry (MESMO 2), an earth system model of intermediate complexity, which consists of a dynamical ocean, dynamicthermodynamic sea ice, and energy moisture balanced atmosphere. The new version has more realistic land ice masks and is driven by seasonal winds. A major aim in version 2 is representing the marine silica cycle mechanistically in order to investigate climate-carbon feedbacks involving diatoms, a critically important class of phytoplankton in terms of carbon export production. This is achieved in part by including iron, on which phytoplankton uptake of silicic acid depends. Also, MESMO 2 is coupled to an existing terrestrial model, which allows for the exchange of carbon, water and energy between land and the atmosphere. The coupled model, called MESMO 2E, is appropriate for more complete earth system simulations. The new version was calibrated, with the goal of preserving reasonable interior ocean ventilation and various biological production rates in the ocean and land, while simulating key features of the marine silica cycle.


Introduction
Here we document development of the second version of the Minnesota Earth System Model for Ocean biogeochemistry (MESMO 2).The first version described earlier (Matsumoto et al., 2008) is based on a non-modular version of the Grid ENabled Integrated Earth system model (GENIE; http://www.genie.ac.uk/), which in turn is based on the computationally efficient ocean-climate model of Edwards and Marsh (2005).This work is independent of the efforts of the GENIEfy project to develop different flavors of the modularised GENIE.
MESMO has a 3-D dynamical ocean on a 36 × 36 equalarea grid with 16 vertical levels.It is coupled with a 2-D energy-moisture balanced model (EMBM) of the atmosphere, dynamic and thermodynamic models of sea ice, and a prognostic model of marine biogeochemistry.It is an earth system model of intermediate complexity (EMIC), a group which occupies an important and unique position within the hierarchy of climate models (Claussen et al., 2002).EMICs represent a compromise between high resolution, comprehensive coupled models of atmospheric and oceanic circulation, which require significant computational resources, and conceptual (box) models, which are computationally very efficient, but represent the climate system in a highly idealised manner.As an EMIC, MESMO retains important dynamics, which allow for simulations of transient climate change, while remaining computationally efficient.The efficiency is achieved by reducing spatial resolution as well as the number and detail of processes compared to high resolution coupled models.
Calibration of the earlier MESMO 1 benefited from multiobjective tuning of the physical model parameters, whereby the mismatch between model-simulated fields and equivalent observed fields was minimised.The equilibrium run of MESMO 1 is well calibrated with respect to oceanic uptake of anthropogenic transient tracers, chlorofluorocarbons, anthropogenic carbon, and radiocarbon (Matsumoto et al., 2008).These tracers reflect ocean ventilation over decades to a century by intermediate waters in the upper ocean as well as Published by Copernicus Publications on behalf of the European Geosciences Union.(Josey et al., 2002) used in MESMO 1 and ECMWF (Trenberth et al., 1989) used in MESMO 2/2E.by the relatively rapid North Atlantic Deep Water (NADW).The model is also well calibrated on centennial time scale with respect to the abundance of natural 14 C ( 14 C) in the deep Pacific and Indian.MESMO 1 has been used successfully in a number of carbon and climate process studies (Lee et al., 2011;Matsumoto and McNeil, 2013;Matsumoto et al., 2010;Matsumoto and Yokoyama, 2013;Sun and Matsumoto, 2010;Ushie and Matsumoto, 2012) as well as in model intercomparison projects (Archer et al., 2009;Cao et al., 2009;Eby et al., 2012;Joos et al., 2012;Weaver et al., 2012;Zickfeld et al., 2012).
A strong motivation for developing MESMO 2 was to investigate the climate-carbon feedbacks involving diatoms, which were not represented in MESMO 1. Diatoms are critical in the ocean carbon cycle, because they are by far the most important agent of vertical transport of carbon from the surface to the deep ocean (Armstrong et al., 2002).Because they can very rapidly grow when conditions are favourable and form mats that facilitate fast sinking, diatoms can account for most of the carbon export production that takes place in the Southern Ocean and more than 50 % globally (Sarmiento and Gruber, 2006).Diatom production is often limited by the availability of silicic acid (Si(OH) 4 ).The silicic acid limitation in much of the low latitudes has its origin in the Southern Ocean, where iron (Fe) deficiency causes Antarctic diatoms to utilise proportionally more Si(OH) 4 for a given amount of NO 3 , so that Si : N utilisation ratio is generally elevated (Sarmiento et al., 2004).The high uptake rate of Si(OH) 4 in the Southern Ocean causes its depletion there.With available data, it is possible to trace this depletion signal originating from the Southern Ocean and spreading to the rest of the ocean via Antarctic Intermediate and Mode Waters (Brzezinski et al., 2002;Sarmiento et al., 2004).These are key features of the observed silica cycle that we chose to simulate with MESMO 2, made possible in part by including Fe as a new tracer in the model.
Another motivation for this work is to have a model that includes a terrestrial scheme, so that the global carbon cycle encompassing the atmosphere, ocean, and land can be simulated.A terrestrial scheme with prognostic land surface albedo would also allow a land albedo feedback in global climate change simulations.In this regard, we make use of the existing model ENTS (efficient numerical terrestrial scheme), coupled previously to GENIE (Williamson et al., 2006).MESMO 2 coupled to ENTS is here referred to as MESMO 2E.In this work, the equilibrium simulations of MESMO 2 and MESMO 2E are described.Key diagnostics of their equilibrium states are summarised in Table 1.

New features of MESMO 2
In the following two sections, we describe in steps the main new physical and biogeochemical modifications and additions which were adopted in MESMO 2 (Table 2).We begin in Sect.2.1 with improvements to the physical model.They include a more realistic local surface albedo specification and seasonal reanalysis winds products from EMCWF.In Sect.2.2, we describe the development of the new biogeochemical model, which includes separating export production by two size classes of phytoplankton.The larger phytoplankton class represents diatoms and requires the iron and silica cycles, which are also described in this section.

New physical features
In MESMO 1, planetary albedo is determined as function of latitude, and ranges from 0.2 at low latitudes to 0.5 at high latitudes.Sunlight reflectivity generally increases with latitude, because sun rays strike the surface at oblique angles, must pass a greater column of atmosphere, and the ground surface more frequently has ice or snow.The planetary albedo in MEMSO determines the reflectivity at the top of the atmosphere with respect to incoming shortwave radiation and follows the original and current formulation in GENIE (Marsh et al., 2011).At the surface, local albedo is modified by the presence of sea ice and can reach a maximum value of 0.7.As there is no land model in MESMO 1, a rather unrealistic situation arises, whereby the surface albedo of the ice sheets on both Greenland and Antarctic is too low, even lower than sea ice located at lower latitudes.So in the first step of MESMO 2 development, we correct this by specifying a land mask with higher albedo (Table 2).Since this increases the global average surface albedo, the latitudinal profile of the planetary albedo is reduced uniformly by 3.5 % to compensate.
The second important change made in MESMO 2 is the replacement of the annual mean wind stress field used in  (Antonov, 2010, Locarnini, 2010).The domains of the deep water masses for calculation of 14 C follow the definitions from the Ocean Carbon Cycle Model Intercomparison Project (Matsumoto et al., 2004): NADW: equator to 60 • N, 1000-3500 m; CDW: 90 • S to 45 • S, 1500-5000 m; and NPDW: equator to 60 • N, 1500-5000 m.Note that all model results except for the 1994 tracer inventories are values at equilibrium, while most of the data constraints are from the postindustrial period.In particular, global mean surface air temperature of 14.0 • C from Jones et al. (1999) represents the 1961-1990 period.The equivalent value is 14.2 MESMO 1 with seasonal wind stress fields to drive ocean dynamics (Fig. 1, Table 2).Following GENIE's core climate model (Edwards and Marsh, 2005), MESMO 1 uses the annual mean wind stress data from Southampton Oceanography Centre (SOC) climatology (Josey et al., 1998).This is replaced in MESMO 2 with the monthly mean wind stress fields from the European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis for the period 1980-1989 (Trenberth et al., 1989).As noted by Josey et al. (2002), the reanalysis winds are approximately 40 % stronger over the Southern Ocean than observation-based winds such as the SOC winds (Fig. 1).As a result, the introduction of the seasonal ECMWF wind stress fields in MESMO 2 greatly strengthens its ocean ventilation, so that the deep ocean In the third step and in an effort to correct this excessive ventilation, we removed the artificial wind stress scaling factor everywhere except in the far North Atlantic.A scaling factor of about 2 was introduced originally by Edwards and Marsh (2005) to achieve sufficiently strong wind driven circulation, and that factor is carried forward in MESMO 1.Using the seasonal ECMWF winds and removing the scaling  completely reduce the excessive deep ventilation significantly and improves the deep 14 C distribution.However, the Atlantic meridional overturning circulation (MOC) becomes too weak (6 Sv; Sv = 10 6 m 3 s −1 ).Therefore, in order to maintain the Atlantic MOC to a reasonable strength in MESMO 2, the wind scaling of 2 is kept northward of 40 • N. Also, the Atlantic-to-Pacific freshwater flux adjustment was increased from 0.2 to 0.3 Sv in the Northern Hemisphere, but reduced in the Southern Hemisphere by the same amount, so that the basis-wide adjustment remains unchanged.With these changes, the Atlantic MOC is 12 Sv in MESMO 2 and 17 Sv in MESMO 2E (Table 1).These freshwater flux adjustments can appear ad hoc, but they compensate for the typically weak transport of moisture in EMBM and are necessary to achieve a reasonably strong MOC in the model architecture of GENIE (Marsh et al., 2011).
In MESMO 1, the winds that drove the different components of the model were different wind products.For exam-ple, the wind stress fields that drove ocean dynamics were not consistent with the wind speed fields that drove air-sea gas exchange.And winds that drove evaporation, precipitation, and the transport of heat and moisture in the 2-D atmospheric model were different still.In version 2, all components of the model, including the terrestrial model ENTS are now consistently driven by the same seasonal ECMWF reanalysis winds.

New biogeochemical features
The limiting nutrients in MESMO 1 are phosphate (PO 4 ), nitrate (NO 3 ), and aqueous CO 2 , whose uptake is governed by Michaelis-Menton kinetics with calibrated half saturation rates.Organic carbon production is assumed to be carried out by only one generic class of phytoplankton, and it is modelled as export production in the top two layers above the fixed compensation depth of 100 m.In MESMO 2, Fe and Si(OH) 4 are added as limiting nutrients, and export production is now modeled through two size classes of phytoplankton, large (LP) and small (SP).We clarify that SP and LP represent downward flux of organic matter from the top two layers and not standing stocks of biomass.
The nutrient dependence of growth for the small class SP is limited by total Fe (FeT), PO 4 , NO 3 , and CO 2 : where K X 's are half-saturation concentrations in the usual Michaelis-Menton formulation of nutrient uptake kinetics (Table 2).As noted in our description of MESMO 1, the justification for including CO 2 as a possible limiting nutrient is that CO 2 availability to phytoplankton may be diffusionlimited (Riebesell et al., 1993).Even though our model experiments show that CO 2 does not limit production under current conditions, we keep this mechanism in place for possible use in the future under different conditions.The K X values used in MESMO 2/E were chosen to be consistent with MESMO 1.For any nutrient, the MESMO1 K X value lies between the LP and SP K X values.Also, smaller values of K X are assigned to SP, so they have a competitive advantage in low nutrient environments over LP.This accounts for the larger surface area to volume ratio that facilitates a faster diffusive uptake of nutrients by SP.
The most limiting nutrient is identified and its uptake rate is determined by the minimum of the above equation.The uptake for other nutrients is related to the limiting nutrient by the particulate organic matter (POM) elemental stoichiometry: P : N : C = 1 : 16 : 117 and the ratios involving Si and Fe are variable as noted below.The growth of both SP and LP is also dependent on light, mixed layer depth, temperature, and biomass; these dependencies remain the same as in MESMO 1. CaCO 3 -forming phytoplankton such as coccolithophorids is assumed to be part of SP in MESMO 2. The dependence of CaCO 3 production on carbonate ion saturation concentration remains the same as in MESMO 1.
The large phytoplankton class LP essentially represents diatoms, so its nutrient dependence term is further limited by silicic acid with an additional term for Si(OH) 4 in Eq. ( 1).Since diatoms are commonly assumed to be competitive when Si(OH) 4 is available, LP increases until Si(OH) 4 is nearly completely drawn down.The rate of Si(OH) 4 utilisation follows the Michaelis-Menton kinetics with a half saturation concentration shown in Table 2. Silicic acid utilisation also depends on the bioavailability of iron such that its uptake relative to nitrate (i.e., the Si : N uptake ratio) increases with decreasing iron (Franck et al., 2000;Hutchins and Bruland, 1998;Takeda, 1998).As a result, the stoichiometry of Si to other nutrients P and C also changes.In MESMO 2, Si : N uptake follows an inverse relation (Si : N = 0.2×10 −9 ×[Fe] −1 , where [Fe] is mol-Fe kg −1 ), with minimum value capped at 0.3 following the data-based estimation (see Fig. 2c of Sarmiento et al., 2004)  temperature and the degree to which the local seawater is undersaturated with respect to the solid phase.The formulation for water column dissolution, including parameter values, follows Ridgwell et al. (2002).Because MESMO 2 is not coupled to a sediment model, any particle that reaches the sea floor dissolves completely and Si(OH) 4 is returned to the overlying water.MESMO 2 calculates the seawater δ 30 Si or the relative abundance of 30 Si compared to the more common and lighter isotope 28 Si.Because marine diatoms fractionate against 30 Si during silicic acid fixation, δ 30 Si is believed to reflect the degree to which Si(OH) 4 is utilised by diatoms (De La Rocha et al., 1998).The fractionation factor during Si(OH) 4 fixation in MESMO 2 is set to 0.9989 or −1.1 ‰ (De La Rocha et al., 1997): where R opal and R Si(OH) 4 are the standard ratios of the rare isotope ( 30 Si) to the most abundant isotope ( 28 Si).
The starting point of MESMO 2 iron cycle is the basic Fe code of the GENIE as used by Annan and Hargreaves (2010).The code follows the seminal modeling work of Archer and Johnson (2000), who gave a prominent role for organic ironbinding ligands in controlling the oceanic iron concentration.While the vertical profile of FeT is nutrient-like (depleted at surface), the deep ocean distribution of FeT does not exhibit the classic Atlantic-to-Pacific increase in nutrients.Although recent work indicates slightly lower FeT in the Pacific versus the Atlantic (Boyd and Ellwood, 2010), the deep ocean concentration has long been considered relatively homogeneous.As assumed by Archer and Johnson (2000), the observed FeT distribution is typically explained and modelled in terms of both iron complexation to organic ligands and removal by scavenging (Johnson et al., 1997).
In MESMO 2, iron is introduced into the ocean via atmospheric dust.For the dust deposition, the depositional flux field projected by the atmospheric tracer transport dust model of Mahowald et al. (1999) is used and re-gridded to the 36 × 36 equal grid of MESMO 2. Uniform iron content in dust of 3.5 weight % with fractional solubility of 0.2 % is assumed.However, rather than apply a uniform solubility of iron in dust, solubility is instead allowed to vary inversely as a function of dust loading, consistent with laboratory experiments and observations summarised by Ridgwell (2001) as well as with a solubility that scales inversely to the square root of dust loading (Baker and Jickells, 2006).Most of the, soluble Fe is quickly bound to organic ligands according to the conditional stability or scavenging schemes of Parekh et al. (2005).FeT is the sum of free dissolved iron (Fe ), which occurs in very small concentrations, and the more abundant ligand-bound iron (FeL).It is partitioned according to: where L is the ligand concentration and set at a constant value of 1 nmol kg −1 everywhere.K ligand , is the conditional stability constant (or binding strength, see Table 2).As the binding strength is increased, iron binds more strongly to ligands, increasing FeL and thus FeT in the water column.Any remaining free iron is quickly scavenged by sinking POM and is modeled, also after Parekh et al. (2005), as: where K sc is the scavenging strength.It is the base scavenging strength K 0 (= 0.079 d −1 ), modified by both a particlespecific fraction and an expression for the abundance of POC with an exponential dependence: Here, for simplicity, we allow scavenging only by POC, and set τ = 0.7.These parameters are taken from the empirical determination of Honeyman et al. (1988).
In MESMO 2, phytoplankton can uptake either Fe or FeL, and the uptake occurs with a variable C : Fe ratio, which is dependent on FeT (Sunda and Huntsman, 1995).Phytoplankton uses iron more efficiently in low iron waters, so that the uptake C : Fe ratio can reach as high as 200 000 in the model.As POM sinks, iron is remineralised along with other nutrients and can subsequently become scavenged again by POM or become ligand-bound and exist in the water column.The scavenged iron which escapes water column remineralisation and reaches the sea floor is then assumed to be buried in the sediments and removed from the model domain.This removal over time will balance the aeolian input at the surface.
Finally, as noted previously (Matsumoto et al., 2010), there was a minor error in the code of MESMO 1 with respect to CaCO 3 remineralisation.It has an exponential formulation with a Q 10 = 2 dependence on temperature; the dependence is similar to POC remineralisation but weaker for CaCO 3 .The original formulation incorrectly allowed more dissolution of CaCO 3 than exists under some conditions, so the fix introduces a cap on the maximum dissolution to prevent this unrealistic situation.In MESMO 2, this fix changes the global CaCO 3 production by less than 1 %.Another minor update in MESMO 2 calculates global ocean CO 2 chemistry immediately prior to model output, so that all outputs reflect the model state at the same time.

Terrestrial scheme ENTS and MESMO 2E
MESMO 2E couples MESMO 2 to the terrestrial scheme ENTS (Williamson et al., 2006), which exists as an optional module within the GENIE framework.It is a simple prognostic model of land biosphere that calculates the exchange of energy, moisture, and carbon between land and the atmosphere.Global fluxes of carbon from photosynthesis, plant and soil respiration, and leaf litter drive carbon stocks of land vegetation and soil.Photosynthesis has dependence on atmospheric CO 2 , water stress, air temperature, and biomass or vegetation fraction.Land vegetation is expressed as fractional coverage with corresponding albedo based on vegetation, soil cover, and soil type.Prognostic variables include vegetation and soil carbon as well as land surface albedo and temperature.Also, as noted above, the seasonal NCEP reanalysis winds, which drove ENTS in Williamson et al. (2006), are replaced in MESMO 2E by the same seasonal ECMWF reanalysis winds that now drive the ocean and atmosphere.The two reanalysis products are actually quite similar, so the replacement has little impact on overall model performance, but it elevates the level of consistency in the new model boundary conditions.
The original ENTS did not have carbon isotopes, which were added in MESMO 2E.During photosynthesis, land plants preferentially fix the lighter 12 C over 13 C and 14 C, so that the biosphere becomes isotopically light and ambient air becomes heavy.In the model, a photosynthetic fractionation factor is set to 0.9815 (i.e., −18.5 ‰) for 13 C and twice as large for 14 C (O 'Leary, 1988).No fractionation is assumed for other terrestrial carbon processes.
Following Williamson et al. (2006), MESMO 2E is calibrated by adjusting a set of freely adjustable rate parameters in ENTS that relate to photosynthesis (k18), vegetation respiration (k24), and soil respiration (k29).The sensitivity of the model results to each of these parameters individually was first determined, then the parameter spaces defined by them were swept in an ensemble of simulations.The choice of parameters was guided by our desire to minimise the difference from published estimates of preindustrial global vegetation, stock (Olson et al., 1983), soil carbon stock (Batjes, 1995) and global carbon fluxes from IPCC (Houghton et al., 2001) (Tables 1 and 2).Our calibration of ENTS changed the land surface properties such that the land surface albedo increased and global temperatures became lower by ∼ 1 • C. In order to compensate for this cooling, planetary albedo in MESMO 2E is reduced from MESMO 2, so that the total change relative to MESMO 1 is increased to −5.5 % when combined with the change already introduced to compensate for the new ice sheet albedo values (Table 2).

Equilibrium runs of MESMO 2 and 2E
Since MESMO is a tool developed primarily to investigate ocean biogeochemistry, the new model versions were calibrated with the goal of having reasonable ocean physics with greater attention paid to reproducing key aspects of the marine silica cycle and ocean biogeochemistry in general.
Equilibrium runs were established for both MESMO 2 and MESMO 2E by spinning up the models with present-day seasonal climatological insolation, winds, and preindustrial CO 2 levels.The time step is 0.01 yr for the physical model and 0.05 yr for the biogeochemical model, so that the latter is called once every five time steps taken by the former.Simulations were run for 5000 yr until deep ocean natural radiocarbon concentrations reached a steady state.The distribution of radiocarbon, in turn, is a useful metric to compare and assess model physics on long timescales.In contrast, the uptake of transient tracers gives a good indication of the ventilation characteristics of the model on decadal timescales.
The Atlantic MOC of MESMO 2 has about the same strength as MESMO 1 at 12 Sv (Fig. 2, Table 1).It is stronger and more realistic in MESMO 2E (17 Sv), which makes the deep Atlantic younger in terms of 14 C, slightly closer to observational mean estimates.However, the difference in natural 14 C between NADW and NPDW is not as large in the new models as in observations or MESMO 1 (Fig. 3).A probably important factor controlling the spread in values is the degree to which air-sea gas exchange occurs when Antarctic Bottom Water (AABW) is formed in the models.In reality, various water mass transformation processes occur on the Antarctic continental shelves and under ice (Foster and Carmack, 1976) with very limited air-sea gas exchange.However, such localised processes are not represented in global models, in which open ocean convection, with significantly more gas exchange, typically plays a significant role in AABW formation.As shown in Fig. 1, the new ECMWF reanalysis winds are significantly stronger in the Southern Ocean than those used in MESMO 1. Gas exchange is therefore enhanced, giving higher 14 C in the Southern Ocean both at the surface and thus at depth.This likely contributes to the smaller 14 C difference between NADW and NPDW in MESMO 2 and 2E (Fig. 3).Other causes of a weak interbasin gradient of 14 C may include excessively strong deep horizontal mixing, which may be a result of insufficient resolution of the bottom topography.Rough bottom topography may constrict near-bottom flow.Another factor may be excessively large sea ice formation, which may force open ocean convection too far north.Still, the NADW-NPDW spread in the new models is about 100 ‰, which is comparable to the Ocean Carbon Cycle Intercomparison Project (OCMIP) models (Matsumoto et al., 2004).
In the postindustrial simulations, we continue from the equilibrium runs, but force the models with observed atmospheric pCO 2 and pCFC-11.On decadal timescales, ocean ventilation is quite reasonable, as judged by the uptake of these transient tracers (Table 1).The uptake for the year 1994 is 100 PgC of anthropogenic carbon for MESMO 2 and 101 PgC for MESMO 2E and 0.56 × 10 6 moles of CFC-11 for MESMO 2 and 0.59 × 10 6 moles for MESMO 2E.These are within the observational constraints of 118 ± 19 PgC (Sabine et al., 2004) and 0.55 ± 0.12 × 10 6 moles CFC-11 (Willey et al., 2004) and compare well against OCMIP models (Matsumoto et al., 2004) ocean, is consistently too shallow in OCMIP and other models (Doney et al., 2004).
In MESMO 1, the global mean surface air temperature is 11.4 • C during the preindustrial equilibrium state (Table 1) and 12.0 • C during the 1961-1990 period in a postindustrial transient simulation with anthropogenic CO 2 forcing.These are much colder than the data-derived, estimated temperature of 14.0 • C for the same 1961-1990 period (Jones et al., 1999).The model-observation discrepancy is significantly improved in the new MESMO versions.For the 1961-1990 period, the global mean surface air temperature is 14.2 • C for MESMO 2 and 14.3 • C for MESMO 2E for postindustrial transient simulations.These are 0.7 • C warmer than their respective preindustrial equilibrium temperatures (Table 1).As would be expected, the coarse resolution of MESMO makes it difficult to sustain strong gradients of, for example, surface air temperature (Fig. 4) or sea surface temperature (Fig. 5).Consequently, air temperature tends to be longitudinally smooth, and the western warm pool is not as warm in the models as observed.
In terms of seasonal sea ice extent, the new models show improvement over MESMO 1, but still overestimate it when compared to observations (Table 1, Fig. 6).It is important to keep in mind that as with global mean surface air temperature, sea ice extent presented in Table 1 represents the preindustrial equilibrium state for the models, not some postindustrial period for the observational constraints.The Arctic sea ice extent in February is 20.3 × 10 6 km 2 in MESMO 2 and 18.5×10 6 km 2 in MESMO 2E, compared to satellite derived estimate of 14-16 × 10 6 km 2 (Table 1).The overestimation is smaller, because global sea ice extent is 8-10 % smaller in year 2000 than during the preindustrial state in the two models.Nevertheless sea ice extends much further south in the far North Atlantic in all versions of MESMO compared to observations (Fig. 6).This reflects the fact that open ocean convection that occurs in reality in the Greenland-Ice Land-Norwegian Seas and leads to the formation of NADW occurs too far south in the models.In the Antarctic, the September sea ice extent is 33.0×10 6 km 2 in MESMO 2 and 25.5×10 6 km 2 in MESMO 2E, compared to satellite derived estimate of 17-20 × 10 6 km 2 (Table 1).Again, the discrepancy would be smaller is smaller for the same reference time.

Ocean biogeochemistry
The values for the new biogeochemical parameters (Table 2) were tuned on the basis of experimental experience with the goal of achieving reasonable global production rates, and spatial distributions of POC, CaCO 3 , and opal (Table 1).We note that the starting point of this tuning was MESMO 1, which was objectively tuned (Matsumoto et al., 2008), and that subjective tuning has yielded superior results to objective tuning in an earlier tuning exercise of GENIE (Lenton et al., 2006) and remains a primary method of tuning (Marsh et al., 2011).In our tuning, we also targeted selected features of the iron and silica cycles.One of the iron targets is the deep water FeT concentrations of 0.6-0.7 nmol kg −1 .Another is  surface FeT concentrations that are low enough for the Si : N uptake ratio to be elevated in the Southern Ocean and the far North Pacific as observed (Sarmiento et al., 2004).Finally, iron has to be limiting production in the Southern Ocean.The two important targets with respect to the marine silica cycle are the export of Si(OH) 4 -depleted waters from the Southern Ocean to the rest of the world ocean via Antarctic Intermediate and/or Mode Waters and the consequent Si(OH) 4 limitation for LP outside the Southern Ocean.Regarding the latter silica target, we note that while observations suggest that diatoms in much of the lower latitudes is limited by silicic acid limitation (Sarmiento et al., 2004), the ecosystem modelling work of Moore et al. (2002) indicates that silicic acid limitation is actually quite small and that iron limitation is much more significant.
The global export production of POC is 11.9 PgC yr −1 in MESMO 2 and 12.5 PgC yr −1 in MESMO 2E (Table 1).These are consistent with a recent synthesis of particle export production, in which Dunne et al. (2007) give their best estimate as 9.6 ± 3.6 PgC yr −1 for POC.Previous estimates ranged between 5.8 and 12.9 PgC yr −1 .Of the simulated global production, the majority is due to LP: 8.7 out of 11.9 PgC yr −1 in MESMO 2 and 9.2 out of 12.5 PgC yr −1 in MESMO 2E.Unfortunately, the proportion of total export production which can be attributed to diatoms is not well understood from observations, with estimates ranging between 20 and 90 % (Sarmiento and Gruber, 2006).As expected, total production is lowest in the nutrient poor oligotrophic gyres (Fig. 7a).The contribution of SP to the total production is highest in oligotrophic gyres, given the lower half saturation values in SP that give it competitive advantage over LP (Fig. 7b).
The simulated global CaCO 3 production is 1.0 PgC yr −1 by MESMO 2 and 0.9 PgC yr −1 by MESMO 2E (Table 1).These are comparable to MESMO 1 but higher than the best estimate of Dunne et al. (2007) of 0.52±0.15PgC yr −1 .They note, though, that constraining the global CaCO 3 export has been controversial and that previous estimates ranged from 0.38 to 4.7 PgC yr −1 with most estimates occupying the lower end of the range.
The global export production of opal is also not well constrained by data.According to Dunne et al. (2007), historical estimates have ranged from 70 to 185 Tmol Si yr −1 (Tmol = 10 12 mol), while their best estimate is 101 ± 35 Tmol Si yr −1 .Global export is 130 Tmol Si yr −1 in MESMO 2 and 139 Tmol Si yr −1 in MESMO 2E (Table 1).The majority of the export production occurs in the Southern Ocean and a secondary peak in the North Pacific (Fig. 7c); this spatial pattern is consistent with Dunne et al. (2007).
In both MESMO 2 and 2E, FeT is higher in the Atlantic than in the Pacific.At the surface, FeT is high in the North Atlantic, which reflects the large aeolian input of Fe from the African Sahara, and also around Antarctica, where there is deep upwelling (Fig. 7d).Since the Si : N uptake ratio by diatoms is inversely related to Fe availability, the Si : N uptake ratio by diatoms is very low in the North Atlantic (Fig. 7e).Conversely Fe limitation in the Southern Ocean and the North Pacific, two well know HNLC regions, cause the Si : N uptake ratio in those regions to be high.These spatial features are consistent with data-based uptake ratio (Sarmiento et al., 2004).At depth, the simulated FeT is approximately 0.65 nmol kg −1 in the Atlantic and about 0.6 nmol kg −1 in the Pacific (Fig. 7f).Available data also show lower FeT in the deep Pacific, which has the oldest waters and thus has experienced the most scavenging (Boyd and Ellwood, 2010).
The subtropical gyres evident in terms of various diagnostics of model production (Fig. 7) clearly have depleted PO 4 concentrations in the top 100 m, although the depletion is stronger in observation (Fig. 8).Also, compared to MESMO 1, both MESMO 2 and 2E have a more pronounced and improved expression of eastern equatorial upwelling in terms of PO 4 , driven by seasonally stronger new wind fields.
In the two new models, Fe is the limiting nutrient in the Southern Ocean for both SP and LP (Fig. 9).For SP, Fe is also limiting in the Arctic, but otherwise, nitrate is the limiting nutrient in much of the world ocean.Diatoms (LP), are mostly limited by Si(OH) 4 outside the Southern Ocean, which is consistent with observations (Sarmiento et al., 2004).There is no Si(OH) 4 limitation in the North Atlantic, where the high aeolian dust flux and thus high surface FeT reduce both the Si : N uptake ratio (Fig. 7e) and the demand for Si(OH) 4 .
Partly as a result of this reduced demand, Si(OH) 4 in the surface subpolar North Atlantic is fairly high in the model in contrast to observations (Fig. 10a and b).In observations, Si(OH) 4 is rather depleted there despite the low Si : N uptake (Sarmiento et al., 2004), because large Saharan iron input and seasonal blooms drive strong opal export.In the models however, diatoms in the North Atlantic do not respond readily to the Sahara iron, since they are more NO 3 limited (Fig. 9) and seasonal variability is typically too small.This modeldata disagreement therefore indicates a limitation of our relatively simple ecosystem model.Elsewhere though, the new models capture observations reasonably well.For example, the surface North Pacific does not experience Si(OH) 4 limitation and has elevated concentrations of Si(OH) 4 in both the models and data (Fig. 10a and b).In the Southern Ocean, the high Si : N uptake rate of Si(OH) 4 in the Southern Ocean causes its depletion there.The degree to which Si(OH) 4 is depleted is commonly expressed in relation to nitrate as Si* (Si* = [Si(OH) 4 ]−[NO 3 ]).Negative values of Si* indicate Si(OH) 4 -depleted waters at the surface (Fig. 10c and d).As seen in available data (Brzezinski et al., 2002;Sarmiento et al., 2004), the new model has this depletion signal spreading out from the Southern Ocean to the rest of the ocean via Antarctic Intermediate/Mode Water (Fig. 10e and f).There may also be a surface contribution to this spreading.We believe that some form of southern sourced intermediate water does exist in MESMO 2 despite its coarse model resolution, because there is a tongue of low salinity water originating from the surface Southern Ocean and penetrating northward (Fig. 10g and h).This tongue is prominent in the Atlantic sector, but less so in the Pacific sector.
We note that these Si* features provided the original motivation to develop MESMO 2. There is a hypothesis that the Si(OH) 4 depletion observed in the modern Southern Ocean can be reversed during times of high dust or iron input such as the glacial periods.The consequent reorganisation of the global silica cycle could help explain some of the variability in atmospheric CO 2 levels (Brzezinski et al., 2002;Matsumoto et al., 2002).There have been a number of paleoceanographic studies that directly attempted to test this hypothesis with as yet no conclusive verdict (Matsumoto and Sarmiento, 2008).The ability of MESMO 2 to reasonably reproduce key Si* features makes the model potentially useful in investigating the biogeochemical implications of climatecarbon feedbacks involving diatoms.Another useful tracer in this regard is seawater δ 30 Si, which becomes heavier in a parcel of seawater as Si(OH) 4 utilisation increases.Measurements of δ 30 Si in deep sea sediments offer an exciting possibility to investigate past changes in the marine silica cycle.Already a number of studies have made measurements of different sedimentary fractions such as opal (Horn et al., 2011) and sponge spicules (Ellwood et al., 2010).
To date, there are two modelling studies of seawater δ 30 Si. First, Wischmeyer et al. (2003) used an ocean general circulation model with a prognostic PO 4 -based export production, which is then related to Si production.Second, Reynolds (2009) used a box model, in which diagnostic PO 4based production is again related to Si production.With a more realistic and mechanistic representation of the marine silica cycle, MESMO 2 makes significant improvements over these earlier efforts.As in Wischmeyer et al. (2003), the −1.1 ‰fractionation during silicic acid uptake in our model imparts a heavy δ 30 Si signal to very near surface waters (Fig. 11a) that correspond to a higher degree of Si(OH) 4 utilisation (Fig. 11b).Our simulated vertical profile of δ 30 Si is consistent with recent water column data from the Atlantic, which show δ 30 Si "from approximately +1.2 ‰ in bottom waters to almost +3 ‰ in the surface mixed layer" with strong enrichment occurring near the surface (Souza et al., 2012).The surface distribution of simulated δ 30 Si (Fig. 11c) resembles the simulated patterns of PO 4 and POC production (Figs.7 and 8) in that the subtropical gyres generally have low nutrients and, therefore, lower overall production but higher fraction of SP and higher degree of nutrient utilisation.However, the δ 30 Si map is more complicated, because surface Si(OH) 4 is decoupled from other nutrients though its dependence on Fe and variable Si : N uptake ratio.In the model, there is a deep δ 30 Si gradient that shows enrichment in the Arctic and the North Atlantic compared to relatively low values in the Southern Ocean (Fig. 11d).A qualitatively similar gradient is also observed in the Atlantic albeit larger (Souza et al., 2012).

The terrestrial biosphere
The calibration of the ENTS parameters in MESMO 2E produces 123 PgC yr −1 for global net photosynthesis, 62 PgC yr −1 for vegetation respiration, and 61 PgC yr −1 for leaf litter and soil respiration (Table 1).The leaf litter flux represents the influx for soil carbon reservoir and soil respiration represents its outflux, so they are equal at steady state.and 60 PgC yr −1 for the other fluxes (Houghton et al., 2001).The equilibrium carbon stock in MESMO 2E is 461 PgC in above ground vegetation and 1319 PgC in soils (Table 1).In comparison, Williamson et al. (2006) simulated stocks of 437 and 1317 PgC for vegetation and soil respectively in the original description of ENTS.The data-based targets for these stocks are 451 PgC (Olson et al., 1985) and 1306 PgC (Batjes, 1995), which include postindustrial land use changes and therefore are likely lower than the corresponding preindustrial stocks.
The spatial distribution in carbon vegetation compares reasonably well to observations (Hall et al., 2005): peak carbon storage in tropical rainforests and secondary maximum in boreal forests (Fig. 12a).However, the relatively simple land scheme of ENTS and the fact that precipitation in EMBM is typically underestimated and causes the continental interiors to become too dry (Marsh et al., 2011) may lead to mismatches with observation.For example, the simulated carbon stock in the boreal forest of Northern Asia may be too low.At the same time, the model captures the main desert regions of the world.Soil carbon distribution reasonably shows high values in northern boreal regions, where low temperatures limit soil respiration, and low values in the tropical regions, where temperatures and thus soil respiration rates are high (Fig. 12b).In terms of 14 C, the vegetation is very close to being in equilibrium with the atmosphere (i.e., ∼ 0 ‰), as expected.Soil is more depleted especially in cold regions such as Alaska and Siberia where low respiration rates lead to longer residence times of carbon (Fig. 12c).

Summary
Two new versions of MESMO are presented here: MESMO 2 and MESMO 2E.The only differences between the two are the coupling of a terrestrial biosphere model ENTS in the latter and the planetary albedo adjustment needed to compensate for the change in land surface albedo caused by the coupling.The physical and biogeochemical modifications described here (Table 2) correct unrealistically low surface albedo values on ice sheets, introduce more seasonality, and allow more explicit representation of the marine silica cycle as compared to MESMO 1.
The use of the same seasonal ECMWF winds to drive all aspects of MESMO 2 removes the inconsistency that existed in the previous version, in which different wind products were used to drive ocean dynamics, air-sea gas exchange, atmospheric heat and moisture transport, and ENTS.Compared to the earlier annual winds, the new seasonal winds impart more momentum to the ocean, especially the Southern Ocean, causing the ocean interior to become excessively well ventilated.This necessitated adjustments in the existing wind scaling factor and interbasin freshwater flux in order to realise distributions of natural 14 C and anthropogenic transient tracers that are consistent with observations.
The implementations of the iron cycle, two classes of phytoplankton, and a dependence of the Si(OH) 4 utilisation on Fe availability are sufficient to simulate key features of the marine silica cycle in MESMO 2. They include Si(OH) 4 limitation for diatoms in much of the low latitude oceans, elevated Si : N uptake ratio in the Southern Ocean and the far North Pacific, Si(OH) 4 depletion in the Southern Ocean, and most importantly the export of Si(OH) 4 depletion to the rest of the world ocean via some form of intermediate and possibly surface waters.These features make MESMO 2 wellequipped for future investigation of climate-carbon feedbacks involving diatoms, which have received much recent attention, especially in the paleoclimate context.

Table 1 .
Key model diagnostics and data constraints.Mean SST and SSS based on the 2009 World Ocean Atlas • C for MESMO 2 and 14.3 • C for MESMO 2E for transient simulations in which the models were forced with observed atmospheric pCO 2 with CO 2 radiative feedback.Similarly, the polar sea ice extent for the year 2000 in the same postindustrial transient simulations is smaller by 8-10 % than at the equilibrium, improving the data-model comparison than presented in the table.

Table 2 .
Model parameters.In MESMO 1, there is only one class of phytoplankton that carries out export production and, therefore, only one half saturation constant for each limiting nutrient.In MESMO 2/2E, the half saturation constants for the large (LP) and small (SP) phytoplankton are respectively larger and smaller than the equivalent values for MESMO 1.