Impact of changes in river fluxes of silica on the global marine silicon cycle: a model comparison

The availability of dissolved silica (Si) in the ocean provides a major control on the growth of siliceous phytoplankton. Diatoms in particular account for a large proportion of oceanic primary production. The original source of the silica is rock weathering, followed by transport of dissolved and biogenic silica to the coastal zone. This model study aims at assessing the sensitivity of the global marine silicon cycle to variations in the river input of silica on timescales ranging from several centuries to millennia. We compare the performance of a box model for the marine silicon cycle to that of a global biogeochemical ocean general circulation model (HAMOCC2 and 5). Results indicate that the average global ocean response to changes in river input of silica is comparable in the models on time scales up to 150 kyrs. While the trends in export production and opal burial are the same, the box model shows a delayed response to the imposed perturbations compared to the general circulation model. Results of both models confirm the important role of the continental margins as a sink for silica at the global scale. Our work also demonstrates that the effects of changes in riverine dissolved silica on ocean biogeochemistry depend on the availability of the other nutrients such as nitrogen, phosphorus and iron. The model results suggest that the effects of reduced silica inputs due to river damming are particularly pronounced in the Gulf of Bengal, Gulf of Mexico and the Amazon plume where they negatively affect opal production. While general circulation models are indispensable when assessing the spatial variation in opal export production and biogenic Si burial in the ocean, this study demonstrates that box models provide a good alternative when studying the average global ocean response to perturbations of the oceanic silica cycle (especially on longer time scales). Correspondence to: C. Y. Bernard (cbe024@uib.no)


Introduction
The marine biogeochemical cycle of silica depends on the weathering of rocks on the continents. This temperaturedependant process releases dissolved silica (dSi) into ground waters and rivers which ultimately may discharge into the ocean. All along the transport route over the continents, the dissolved silica may be used by fresh water diatoms and plants to build up biogenic forms of silica (bSiO 2 ) (Conley and Schelske, 2001). Part of the bSiO 2 is deposited in riverine sediments, but a significant fraction reaches the coastal waters thus contributing to the total riverine input of silica (Conley, 1997). In coastal waters, silica input creates favourable conditions for diatom production, which can account for up to 75% of total primary production . The high level of opal production combined with the shallow settings, commonly allow a rapid settling and efficient burial of biogenic silica in coastal sediments (DeMaster, 2002). In the open ocean, an estimated 50% of the opal produced is exported from the euphotic layer (Van Cappellen et al., 2002) and only about 3% is permanently incorporated in seabed sediments (Tréguer et al., 1995). The burial of biogenic silica in coastal and open ocean environments varies spatially. Despite their small surface area (8% of the total area of the ocean), continental margins are suggested to account for a significant percentage of the total accumulation of bSiO 2 in sediments (DeMaster, 2002;Laruelle et al., 2010).
Climate and land surface hydrology play an important role in controlling silicate weathering and the transport of silica to the ocean. As a consequence, river inputs of silica vary on geological time scales (White and Blum, 1995). For example, geological events such as the uplift of the Himalayan plateau may have accelerated silicate weathering during the late Cenozoic by increasing the exposure of crustal rock and by increasing the monsoon regime (Raymo, 1991). Changes in river inputs of silica likely also affected marine Si cycling over glacial-interglacial cycles. Thus, Tréguer and Pondaven (2000) and Ridgwell (2002) suggest that increased river input of silica during the last glacial maximum, possibly combined with increased dust input (Harrison, 2000) may have contributed to the glacial/interglacial changes in the oceanic carbon pump and atmospheric CO 2 concentrations.
More recently, human activities (damming, land use practices, deforestation, and the introduction of invasive species) are affecting the natural terrestrial cycle of silica and its delivery to the ocean (Conley et al., 1993;Humborg et al., 2000;Ragueneau et al., 2005). The effects are particularly evident in changes in nutrient ratios involving Silica (Si), Nitrogen (N) and Phosphorus (P) in surface waters of many coastal areas in temperate regions over the past decades (Conley et al., 1993;Rabalais et al., 1996). The observed decrease of Si:N and Si:P ratios in these areas are the combined effect of increased human inputs of N and P to rivers and decreased river loads of Si due to retention of biogenic Si in reservoirs behind dams. In tropical regions, dissolved Si concentrations appear to be less affected by anthropogenic factors and climatic, geological and geomorphological factors likely are more important (Jennerjahn et al., 2006). Note that deforestation increases the continental input of silica to the ocean by increasing dissolved silicate losses from vegetation (Conley et al., 2008). However, it is uncertain what role deforestation plays in counteracting the worldwide decline in river silica input. The major consequence of changes in nutrient ratios in temperate regions is a shift in planktonic species composition in the near coastal zone, with flagellates, cyanobacteria and other non-siliceous phytoplankton replacing diatoms (Egge and Aksnes, 1992;Humborg et al., 2000). In some cases, toxic blooms of harmful algal species may develop and strongly impact coastal ecosystems and fisheries (Roelke, 2000).
Several modelling tools have been developed to investigate the dynamics of the marine silica cycle at the global scale and its sensitivity to natural and anthropogenic perturbation. These tools consist of mass balance or box models (Laruelle et al., 2010;Yool and Tyrrell, 2003) and Global Biogeochemical Ocean General Circulation Models (Bernard et al., 2009;Heinze et al., 2003). While box models mostly rely on first order kinetic rate laws (Mackenzie et al., 1993), general circulation models include a more robust and mechanistic description of biogeochemical and physical processes in the ocean (Heinze et al., 2003). Both approaches have been used to gain insight in marine silicon cycling and its response to river inputs of silicon. Results of global scale box modelling of the silica cycle for the coming century, for example, indicate that enhanced biogenic silica dissolution due to global warming may enhance silica availability in aquatic systems. Ultimately this may allow coastal siliceous productivity to recover from the downward trend caused by river damming (Laruelle et al., 2010). Inclusion of present-day river inputs of Si in a high resolution general circulation model indicate that the effects of river inputs on coastal marine Si cycling are most pronounced in hot spots, such as the Amazon plume, the Arctic and Southern Asia (Bernard et al., 2009). Given that box models are more easily accessible and less computationally demanding than general circulation models, and to increase their use as a prognostic and predictive tool, it is of interest to compare the response of both types of models to similar perturbations, both on human and geological time scales.
In this study, we make such a comparison using the box model for the global Si cycle of Laruelle et al. (2010) as well as the Hamburg Ocean Carbon Cycle Model in its annually averaged coarse resolution version for long time integrations (HAMOCC2) and the high resolution version that includes the continental margins and a more detailed silica cycle (HAMOCC5). We first compare the marine Si budget inherent to each model and the model assumptions. Then long term (0-15 kyrs) and short term (0-150 yrs) simulations are performed to test the response of the models to changes in river inputs of nutrients. We demonstrate that the trends at the global scale obtained with the general circulation model (GCM) and the box model are surprisingly similar both for short and long term simulations. The critical role of the continental margins in the global cycling of silica as well as their higher reactivity to riverine perturbations is highlighted using results of the box model and HAMOCC5.

The box model
The box model used in this study, for both short-term and long-term simulations, is described in Laruelle et al. (2010). The model is based on an updated budget of the global biogeochemical cycle of reactive silicon, including both the terrestrial and oceanic realms. The Earth surface is divided into 4 compartments along the land ocean continuum: continents (box 1), proximal coastal zone (box 2), distal coastal zone (box 3) and open ocean (box 4). Coastal regions receive particular attention given their role as a filter between the continents and the ocean. The proximal coastal zone consists of large bays, the open water parts of estuaries, inner deltas, inland seas and coastal marshes (Smith and Hollibaugh, 1993;Woodwell et al., 1973). The distal zone comprises the rest of the continental margins up to the shelf break.
In each compartment, an estimate of the amount of dissolved (dSi) and biogenic silica (bSiO 2 ) was established for the water column and the upper layer of the sediment. Yearly averaged fluxes between these silica reservoirs were based on previous budgets (Alexandre et al., 1997;Conley, 2002b;DeMaster, 2002;Tréguer et al., 1995) or were estimated independently. The Si cycle was linked to a steady state hydrological cycle and the advection fluxes for dSi and bSiO 2 were calculated as follows; where F ij and Q ij are the fluxes of reactive Si and water from reservoir i to reservoir j , respectively, M i is the mass of dSi or bSi in reservoir i, and V i is the volume of the reservoir. The remaining transport fluxes correspond to sedimentation and deposition of bSiO 2 in the sediment and the efflux of dSi from sediments to the water column. Simple first order kinetic relationships were obtained for all fluxes that were not linked to the water cycle by deriving a rate constant k ij from the steady state Si budget using the values of the flux F ij and the size of the reservoir from which the flux originates; The reaction network was completed by including the input fluxes through terrestrial rock and seafloor weathering and hydrothermal activities (Tréguer et al., 1995) and the ultimate sinks are burial of bSiO 2 within freshwater and marine sediments (Conley, 2002a;DeMaster, 2002;Tréguer et al., 1995) as well as reverse weathering reactions in shelf sediments (Mackenzie and Garrels, 1966).

HAMOCC2
In order to test the long-term effect of changes in Si supply to the ocean, we employ the HAMOCC global biogeochemical ocean model (Maier-Reimer, 1993 in its computationally efficient annual average version "HAMOCC2s" (Heinze et al., , 2003. Annual average velocity and thermohaline fields are taken from the Large Scale Geostrophic dynamical ocean general circulation model with climatological atmospheric data (details are given in Winguth et al., 1999, for their "interglacial first guess" circulation). The effect of deep convective mixing at high latitudes is represented in the annual average velocity field, which is used for transporting the dissolved tracer substances within the model water column. The horizontal resolution is 3.5 • × 3.5 • . The water column is structured into 11 layers (centered at 25, 75, 150, 250, 450, 700, 1000, 2000, 3000, 4000, and 5000 m). The bioturbated top sediment zone is structured into 10 layers, with interfaces at 0, 0.3, 0.6, 1.1, 1.6, 2.1, 3.1, 4.1, 5.1, 7.55, and 10 cm down from the water sediment interface.
The biogeochemical model includes the following processes: air-sea gas exchange, biogenic particle export production out of the ocean surface layer, particle flux through the water column and particle degradation by dissolution as well as remineralisation, transport of dissolved substances with ocean currents, deposition of particulate constituents on the ocean floor, pore water chemistry and diffusion, advection of solid sediment weight fractions, bioturbation, and sediment accumulation (export out of the sediment mixed layer). The model predicts the following tracer concentrations in the ocean water column, the sediment pore waters, and the solid sediment: water column -DIC (dissolved inorganic carbon), POC (particulate organic carbon), DOC (dissolved organic carbon), CaCO 3 (calcium carbonate or par-ticulate inorganic carbon) of 12 C and 13 C, dissolved oxygen O 2 , dissolved PO 3− 4 as biolimiting nutrient, silicic acid Si(OH) 4 and opal (biogenic particulate silica bSiO 2 ); sediment pore waters -the same dissolved substances as in the water column; and solid sediment -clay, CaCO 3 , opal, and organic carbon. For inorganic carbon chemistry, the dissociation constants of carbonic and boric acid according to Mehrbach et al. (1973), the solubility product for CaCO 3 after Ingle (1975), and the pressure dependencies of Edmond and Gieskes (1970) were applied.
The ocean surface and water column processes and the pore water chemistry are parameterized as described in previous applications of the model (Heinze et al., , 2003. Opal export production and respective depletion of silicic acid in the ocean surface layer are simulated by Michaelis Menten nutrient uptake kinetics. The higher maximum nutrient uptake velocity is higher for Si than for P, leading to a concentration of opal export in upwelling regions. Through this approach, silicic acid is depleted in surface waters before phosphate is completely used up through primary production. The resulting variable Si:C ratios correspond well with trends as given by Brzezinski (1985). The overall coupling of the P and Si cycles in the model is thus governed by vertical (upwelling) velocity and sedimentation rates. For the opal flux through the water column, an implicit numerical algorithm is used involving an independent choice of the particle sinking velocity and the opal dissolution rate. Pore water chemistry follows Archer et al. (1993) "burial = rain minus re-dissolution", but allows for time dependent exchange with ocean bottom water in the free water column (Heinze et al., 2003) and includes efficient numerics for the vertical sediment advection following Maier-Reimer et al. (2005). The re-dissolution constant of opal within the bioturbated sediment zone was adjusted to be lower than the one for the water column. This is in line with the presumed alteration of biogenic silica particles during their aging process in the water column and surface sediment (Van Cappellen et al., 2002). Bioturbation is implemented through diffusion of solid materials, where non-local mixing and particle size dependent mixing have been neglected. The early diagenetic model is based on the concept of conservation of volume (or geometry) where all sediment layers do not change their geometric shape according to a fixed porosity profile (after Ullmann and Aller, 1982) with time. Thus, gaps in the solid sediment are instantaneously closed, either by shifting of material from the respective sediment layer above (in the case where rain exceeds re-dissolution) or from erosion of material from the layer below.
The model is initialized with clay sediment only and run to quasi-equilibrium after 120 000 years of integration before any sensitivity experiments are started. At model equilibrium, the global burial rate of Si corresponds to the given input rate from riverine material. However, the local opal production, deposition, re-dissolution and burial (sediment accumulation) rates are prognostic.

HAMOCC5
HAMOCC5 was used for short timescale simulations (150 yr). Given the high resolution of its grid and detailed description of the Si cycle, the model is not suitable for long term simulations.
HAMOCC5-MPIOM results from the interactive coupling of the Max Plank Institute Ocean Model (MPI-OM), a full primitive equation dynamical ocean model which computes thermohaline circulation and assures the advection and diffusion of biogeochemical tracers of the HAMburg Oceanic Carbon Cycle Circulation Model. The two fully coupled models receive the same radiative forcing. HAMOCC5 includes a more comprehensive description of biogeochemical processes. The model grid used in this simulation is an orthogonal curvilinear C-grid with an average resolution of 3 • . To optimize calculations, the North Pole is artificially located over Greenland and the South Pole over Antarctica. The resulting resolution is 29 km in the Arctic to about 390 km in the Tropics. The water column is divided into 40 vertical levels whose thickness gradually increases with depth, from 12 m at the surface to a maximum of 600 m in the deep ocean. This resolution resolves the continental margins, although in a coarse manner. In this study, the continental shelf is represented by grid cells shallower than 1500 m, of which the integrated surface covers 8% of the world ocean (27.10 6 km 2 ). The average depth is 530 m, 50% of this surface is shallower than 300 m. The time step is 0.1 day.
Only the features relevant to the silica cycle will be described here. For a full description of HAMOCC5 and MPI-OM, we refer to the technical reports that are available online (Maier-Reimer et al., 2005;Wetzel et al., 2010). Riverine inputs of nutrients were implemented according to Bernard et al. (2009) using the COSCAT segmentation of the coast line .
In HAMOCC5, opal production is computed as a fraction of the dead material produced by zooplankton and phytoplankton (including the pseudo faeces); the model calculates opal and calcium carbonate production as the two fractions of the non living part of the detritus, the shells. A main difference with HAMOCC2 is that the opal export production is dependent on primary production. The competition between opal production and calcium carbonate production is regulated by the silica availability (Leynaert et al., 2001;Paasche, 1980). It is assumed that phytoplankton consists of diatoms, coccolithophorids, and flagellates. As diatoms are known to be the fastest competitors (Egge and Aksnes, 1992), opal production by diatoms will be preferred to CaCO 3 production if sufficient dSi is available in the ocean surface layer.
The flux of particles through the water column redistributes phosphorous, silica and associated tracers during sinking, thus enriching the deep waters in nutrients. In the default version of the code, particles have constant sinking speeds, w DET , w CaCO3 , w Opal and w Dust for organic detritus, CaCO 3 , opal and clay, respectively. The export production is computed as the POC (particulate organic carbon), opal and CaCO 3 leaving the euphotic layer, i.e. the material sinking below 90 m depth in the surface ocean. Remineralisation of opal and CaCO 3 occurs during sinking of particles after they have left the euphotic layer. As the model only computes the last step of the living part of the silica cycle, hardly any of the opal produced in the euphotic layer is remineralized in the surface ocean. This simplification implies that primary production of opal is not computed in HAMOCC5, and we therefore only refer to exported opal production.
Fluxes from the bottom ocean layer in each ocean grid cell provide the boundary condition for the sediment module that includes 4 sediment weight fractions and 12 layers following . The sediment module computes the accumulation of deposited material on the sea floor as well as remineralisation in the sediments and the release of redissolved tracers to the lowest level of the water column.

Model comparison
The box model and the two GCMs compared in this study present many conceptual and structural differences. The major differences between the box model and GCMs are the inclusion of other elements than silica and the spatially-explicit description of biogeochemical processes (Table 1). Note also that the GCMs differ in their description of biogeochemical cycling of silica and in the representation of the continental margins. While coastal zone processes are resolved in the box model and roughly in HAMOCC5, they are not included in HAMOCC2. Thus, given the computational demands of HAMOCC5, only the box model allows the assessment of the effects of coastal zone processes on the long-term silica cycle.
The steady state budgets for Si in all three models show general similarities but also some major differences (Fig. 1). Examples of similarities are the rates of sediment burial and benthic recycling in the coastal zone in HAMOCC5 and the box model and the total ocean burial in the box model and HAMOCC2. Differences are observed in the process rates in the euphotic zone (0-100 m) and intermediate waters (100-1000 m). In addition, the burial in the open ocean in HAMOCC5 is much higher than in the other models. Conceptually, it must be remembered that the steady state Si cycle of the box model is the hypothesis on which the model itself is built. The budgets calculated by the GCMs, in contrast, result from calibration and optimization processes (Heinze et al., 2003). The 3 models are fed with similar inputs of dSi  to allow a better intercomparison. The box model, also includes other sources of Si: bSiO 2 (Conley, 1997), ground water inputs (Slomp and Van Cappellen, 2004) and aeolian dust deposition on the open ocean (Tréguer et al., 1995).
The spatial resolution of HAMOCC2 does not allow the coastal zone to be resolved and its global budget is, therefore, the most simple. 211 Tmol yr −1 of Si are exported as opal from the euphotic layer of the ocean. About two thirds of this amount redissolves in the water column before reaching the sediment. To ensure a steady state, a massive recy-cling takes place, leading to an efflux of 69.7 Tmol yr −1 from the sediments. The net burial, perfectly balances the riverine inputs (6.1 Tmol yr −1 ). HAMOCC5, on the other hand, can not be run to full steady state including the slowest reservoirs because of overwhelming computation costs associated with its higher spatial resolution, a more complex biogeochemical module and a much higher temporal discretisation. In this model, 10% of the export production takes place on the continental margins where 1.2 Tmol yr −1 of Si is buried. Out of the 84.2 Tmol yr −1 of opal exported from the euphotic layer in the open ocean, 35.4 Tmol yr −1 reaches the sediment. This compartment is not at steady state and the burial is 16.7 Tmol yr −1 , which is likely too high. However, given the long residence time of Si in the ocean, this approximation is not expected to affect the trends in simulations on the time scale of a few centuries. The box model not only calculates export production but also the total dSi uptake. This number is not calculated in the GCMs but is based on the study of Tréguer et al. (1995). In the model, it is assumed that 20% of the 240 Tmol yr −1 of dSi uptake takes place on the continental shelf. Due to rapid pelagic recycling, only 7.7 Tmol yr −1 enters the sediment, which is close to the estimate of 6.6 Tmol yr-1 obtained with HAMOCC5. However, the coastal burial flux is more than twice as large in the box model as in HAMOCC5 (2.7 Tmol yr −1 versus 1.2 Tmol yr-1). Note that this former number includes the coastal reverse weathering of 1 Tmol Si yr −1 while the additional shelf burial in the box model is 1.7 Tmol yr-1. In both models burial fluxes are significantly higher than Tréguer's estimate of 0.6 Tmol/yr (1995), thus supporting the suggestion of De-Master (2002) that silica burial in shelf sediments has been underestimated so far.
In the box model, only 50% of the opal primary production is recycled in the euphotic layer of the open ocean (based on Tréguer et al., 1995;Van Cappellen et al., 2002). The resulting export production of opal is similar in both models (103 Tmol yr −1 and 84.7 Tmol yr −1 , respectively). The dissolution in the remainder of the water column is 75%, which compares well with values for the other models (64% in HAMOCC2 and 58% in HAMOCC5). Only 4.1 Tmol yr −1 of a total deposition of 24 Tmol Si yr −1 remains permanently in the sediment. This implies a Si recycling in the sediment of 83% which lies between those of the GCMs (92% and 53% for HAMOCC2 and HAMOCC5, respectively).

Long-term effects of changes in river inputs (HAMOCC2 vs. box model)
The long term response of the box model and HAMOCC2 to changes in river inputs of silica was investigated through a set of simulations run over 150 kyrs. The perturbations selected were taken from  using updated river fluxes as described in Bernard et al. (2009): -Simulation 1: reduction to 75% of Si inputs.
-Simulation 4: step function simulating a 10-fold increase Si inputs during 20 kyrs after 10 kyrs of unmodified inputs and followed by 20 kyrs with no inputs. After 50 kyrs, the Si delivery returns to its original level. For each of the 3 first scenarios, an initial spin-up period of 10 kyrs with unmodified inputs was applied before the perturbation.
Export production of opal and burial of Si in the sediment were used as indicators of pelagic and benthic Si processing in the box model and HAMOCC2 as done by . The limit of 100 m is used to compute export production. Export production provides insight into the biological activity of diatoms. Sediment burial is the ultimate sink for the marine biogeochemical cycle of Si and is a key process in determining the long term mass balance of silica in the ocean. In the box model, all sediment sink terms are included, and the silica lost to reverse weathering on the continental shelf is considered as part of the "burial term". Both models are fed with comparable Si deliveries (6.2 and 7.4 Tmol Si yr −1 for

Fig. 2. Long term scenarios Box Model versus HAMOCC2: opal export production (left) and opal burial (right): -Run 1: The riverine Si input is decreased by 25% (a-b) -Model experiment Run 2: riverine Si input is increased by a factor 4 (c-d).
HAMOCC2 and the box model, respectively), and the inputs are balanced by Si burial at the beginning of the simulations.
Simulation 1 shows that a 25% reduction in Si inputs induces comparable decreases in both models for export production and sediment burial ( Fig. 2a and b). While the response in export production represents −10% for HAMOCC2, it reaches −25% in the box model. The quasiinstantaneous drop following the perturbation in the box model is the result of the dynamics of Si in the coastal zone; the characteristic residence time for Si on the continental margins is significantly shorter (2.3 kyrs) than that of Si in the open ocean (17 kyrs). As a consequence, the box representing the coastal zone in the box model responds to each perturbation much faster. This is reflected in the steeper slope of the lines for the box model in Fig. 2. The same general decrease is observed in opal burial. While HAMOCC2 appears to reach a new steady state after only 20-30 kyrs (despite a residence time for Si of 23 kyrs), it takes about 100 kyrs for the box model (with the shorter residence time of 17 kyrs) to reach a new equilibrium. This suggests that the higher spatial resolution and shorter residence time of the individual grid cells in HAMOCC2 compared to the box model, allows a faster collective response of the whole ocean system to a perturbation. A similar lowering of the response time as the result of collective behaviour, has been demonstrated for box models of the global carbon cycle (Lasaga, 1980). Simulation 2, assuming a 4-fold increase of river Si inputs, (Fig. 2c and d) yields an opposite response, with increases in export production and sediment burial for the box model as well as HAMOCC2. Again, HAMOCC2 reaches a steady state much faster than the box model (after ∼30 kyrs as opposed to ∼120 kyrs), and the magnitude of the increase is less than that of the box model. Export production reaches a maximum at 150% of its original value in HAMOCC2, while it reaches a maximum at 300% in the box model. This is partially explained by the structure of the model itself which only allows dSi availability to limit opal production. Hence, dSi uptake (not shown) and export production in the box model follow the inputs of Si delivery to the ocean almost linearly, whereas in HAMOCC2, phosphorus may also limit diatom growth (see Table 1). Despite the difference in export production, opal burial increases by a factor of 4 in both models. This can be explained by a lower redissolution of the sinking bSiO 2 in HAMOCC2 compared to the box model.
Simulation 3 (without any Si inputs) ( Fig. 2e and f) drives export production and opal burial in both models towards zero as the ocean becomes depleted of Si. Export production decreases faster in the box model than in HAMOCC2, reaching only 25% of its initial value after less than 50 kyrs. The opposite trend is observed with Si burial, which drops significantly faster in HAMOCC2 than in the box model at the beginning of the simulation. However, more than 150 kyrs are required for either model to remove all oceanic Si through burial.
Simulation 4 is a step function, imposing a very strong increase of riverine Si input (10-fold) followed by a shutdown. This results in a strong increase in export production and opal burial followed by a significant drop 50 kyrs after the beginning of the simulation in both models. In agreement with the previous results, export production reaches much higher values; the increases in sediment burial in both models remain comparable in magnitude in the box model. However, the temporal response of the models, strongly differs because of the smoother behaviour of HAMOCC2. In the box model, the short residence time on the continental shelf leads to very steep increases in production and sediment burial during the first centuries following the perturbation and slower rates later on. Between 30 and 50 kyrs, without Si inputs, HAMOCC2 reaches slightly lower values for both fluxes than the initial conditions. Despite sharp initial de-creases, the box model still exhibits higher export production (ca. 200 Tmol Si yr −1 ) and sediment burial (18 Tmol Si yr −1 ) after 50 kyr than at the beginning of the simulation. This difference between the models leads to opposite trends during the remainder of the simulations while they slowly return to their original steady states, reached at 80 kyrs and 150 kyrs for HAMOCC2 and the box model, respectively.
Overall, both models present similar qualitative responses to major long term variations in silica inputs from the rivers. Quantitatively however, discrepancies are observed: HAMOCC2 shows a faster return to steady state than the box model while the latter always exhibits larger changes in export production. In fact, the box model shows a quasilinear response to the imposed perturbations. Hence, its export production varies within a significantly larger range of values than that of HAMOCC2. Nonetheless, the long term responses of both models are comparable and similar with respect to opal burial. The differences in transient behaviour mainly concern export production and essentially affect the first 20-30 kyr of the simulations. The two major differences between both models are due to the higher sensitivity of the export production in the box model to Si inputs and the higher resolution in the GCM which allows a faster return to steady state despite the similarity in the overall residence time of Si (Table 1).

Short-term effects of changes in river inputs (HAMOCC5 vs. box model)
The effects of short term changes in the riverine input of silica were assessed by comparing results of 3 scenarios (Simulations B-D) for river input of nutrients to those of a reference scenario (Simulation A; Table 2). In simulation B, the effects of a perturbation of the terrestrial silica cycle induced by river damming was tested (Laruelle et al., 2010). In the box model, this scenario involves changes in the dSi and bSiO 2 input that rely on projected changes in the number of dams until 2025 (Gleick, 2003), combined with a relation between global water use and the number of new dams (Rosenberg et al., 2000), and an estimate of global population changes throughout the 21st century. The HAMOCC5 equivalent scenario was implemented by imposing a reduction of the riverine dissolved silica flux. The coefficient was computed as the silica flux reduction ratio at the proximal/distal interface in the box model. In addition, a scenario was run without any silica input for the box model and HAMOCC5 (Simulation C), and without any input of other nutrients (nitrogen and phosphorus) for HAMOCC5 alone (Simulation D). HAMOCC5 requires a spin up time of 100 yrs. The damming scenario (Simulation B) assumes a maximum decrease of 17% in riverine silica input to the ocean in the year 2100. At the scale of the global ocean, this does not cause a significant decrease of the opal export production and burial ( Fig. 3a and b). Thus, the box model shows a 1% decrease in export production, while the decrease is insignificant compared to the inter-annual variability for HAMOCC5. The response to a total shutdown of riverine Si (Simulation C) input has a much stronger effect, with a decrease in export production of 7% for the box model and 6% for HAMOCC5. Switching off all riverine nutrients in HAMOCC5 (Simulation D), leads to a slightly stronger response, reflecting the additional dependency of opal production on the availability of other nutrients.
The Si cycle on the continental margins is significantly affected by changes in river inputs of nutrients ( Fig. 3c and d).
Since the box model does not explicitly compute shelf opal export production, only results for HAMOCC5 are shown in Fig. 3c. The decrease in export production is difficult to quantify because of the inter-annual variability but the relative change is larger than for the open ocean. The strongest response to the damming scenario is observed for opal deposition in the box model (Fig. 3d). At the end of the simulation, the box model shows a 0.5 Tmol reduction of the opal deposition (−5.8%), compared to a 0.1 Tmol decrease for HAMOCC5 (−1.5%).
The production and export of bSiO 2 in the box model is only limited by dSi availability. This makes the box model much more sensitive to variations in the riverine input of silica than HAMOCC5. The multi nutrient limitation (N, P, Fe and Si) of HAMOCC5 makes it more stable and reduces the effects of changes in riverine silica alone. The modest response of the GCM shows that opal production on the shelf does not suffer only by silica limitation. The opal production is limited by other nutrients such as N, P or Fe. As a consequence, switching off the riverine input of N and P (Fig. 3d) causes a stronger decrease of the opal export production (−22%) than switching off the riverine silica input alone (−16.6%).
A major difference between the two models is the description of riverine input of silica. The box model explicitly receives both riverine bSiO 2 and dSi, while HAMOCC5 is only fed with dSi. In rivers, lakes and artificial reservoirs, damming leads to a stronger trapping of bSiO 2 than of dSi, despite the fact that loads of both components are affected by the increase in water residence time and the strong link between the cycles of bSiO 2 and dSi. This feature, observed in the field (Humborg et al., 2000), and implemented in the box model, modifies the particulate/dissolved ratio for silica. As a consequence, the lower relative input of bSiO 2 further enhances the drop in bSiO 2 sedimentation in the distal zone of the box model while HAMOCC5 can not capture this process. However, the buffer effect of the proximal zone in the box model, limits the variations of the particulate/dissolved ratio to marginal changes (a few percent at most). Overall, the major factor explaining the differences in the response of the models is the inclusion of multiple nutrient limitations in the GCM.
Shelf opal deposition drops by 16.6% when cutting off the riverine input of silica in HAMOCC5 (Fig. 3d). The drop is 35.7% for the box model, which is higher than when cutting off of all riverine nutrients in HAMOCC5 (21.9%). In comparison, the box model's linear response to changes in silica supply highlights the buffer effect of the multi nutrient limitation of the primary production in HAMOCC5. In the box model, only silica drives bSiO 2 production, and a change in the riverine input has a direct effect on the opal production. In HAMOCC5, changes in the dSi availability will have little effect on the opal production, unless dSi itself is the limiting factor driving the opal production. The explicit spatial component of HAMOCC5 allows the visualization of the complex interaction of lateral advection of riverine silica with regional nutrient limitation. A global map of opal deposition obtained from HAMOCC5 for a 150 yr simulation (Fig. 4a) illustrates the large spatial heterogeneity of opal deposition including the major role of the Austral Ocean, the Equatorial Pacific upwelling as well as coastal upwellings off Mauritania, Peru, Chile, western South Africa (Benguela) and eastern New Zealand. The change in opal deposition relative to the reference run for a model scenario without a riverine Si contribution (Sim C) (Fig. 4c) demonstrates the contribution of terrestrial Si input to the marine silica cycle. Although most of the opal sedimentation supported by the riverine input occurs near to the coast, the response in the Pacific differs from that in the Atlantic Ocean. While opal deposition decreases in the Sargasso Sea and the Western North Atlantic, a minor increase is observed in the central Pacific. Results of the model run without riverine nutrients (Sim D) (Fig. 4d) display a decrease of opal deposition of at least 20 mmol m −2 yr −1 for the whole North Atlantic Ocean. This illustrates the role of regional limitation: in the Atlantic Ocean, primary production is nitrogen limited and thus will be affected by a drop in riverine inputs of N. In contrast, the Pacific is largely limited by iron (except for the equatorial region), thus a change in the Si, N or P will have a smaller impact.
As discussed in Bernard et al. (2009), the effect of changes in silica input from rivers largely depends on the location of its release to the ocean. Thus, a small reduction of the riverine input of silica, as prescribed for the river damming scenario (−17%), will have a minor effect in the Pacific Ocean. The main effect is seen in the Gulf of Bengal where a local drop of more than 50 mmol m −2 yr −1 is observed, and in the Amazon plume and the Gulf of Mexico. The reduced input of the Congo River and other minor rivers of Central Africa lowers opal deposition in the central South Atlantic Ocean by 25 mmol m −2 yr −1 .
As an advantage compared to the box model, the GCM describes the circulation in the North Atlantic Ocean and provides insight in the fate of the Amazon River plume. The strong current moving northwest along the coast of Brazil distributes opal deposition into the Caribbean Sea and mixes with the silica rich waters in the Gulf of Mexico. The Amazon River alone contributes to a large fraction of the offshore North Atlantic pool of silica (Shipe et al., 2006;Froelich et al., 1978); with a flux of Si that is estimated to be on the order of 1-1.5 Tmol Si yr −1 . Such a massive input likely supports most of the North Atlantic siliceous uptake (DeMaster et al., 1996).
While the effect of damming on Si dynamics on the global scale is limited (Fig. 3), the distribution of opal burial for this scenario (Fig. 4d) suggests that locally, effects may be significant, even in the open ocean. The regional details of the riverine silica contribution make the general damming scenario assumed here somewhat unrealistic, however. Indeed, a homogeneous increase in river damming all over the world is unlikely. A realistic approach would take into account regional characteristic (e.g. population, land use, water demand) for a spatially explicit global damming scenario. Unfortunately, such a scenario is unavailable at the present time.

Impact of river inputs of silica on the coastal and open ocean
On time scales of centuries, the fate of riverine dSi inputs in the marine environment strongly depends on the extent of biological utilization of the dSi in the coastal zone (Bernard et al., 2009). When other nutrients, such as nitrogen and phosphorus, are abundant, dSi will be consumed locally and will fuel diatom blooms and biogenic Si production. N and P limitation, in contrast, will limit biological activity and will enable lateral export of dSi. Thus, changes in the riverine and groundwater input of N and P as observed in many coastal systems over the past decades (Rabouille et al., 2001;Ver et al., 1999a, b;Slomp and Van Cappellen, 2004) have a direct impact on Si cycling in coastal seas and continental margins.
Resulting changes in Si/N and Si/P ratios are likely responsible for dramatic shifts in phytoplankton species composition in many coastal systems (e.g. Humborg et al., 2000). Our work with HAMOCC5 highlights this interdependency of the various nutrient cycles and its consequences for opal production. It also confirms results of earlier work on impacts of river nutrients on ocean biogeochemistry emphasizing the role of regional nutrient limitation (Da Cunha et al., 2007). We also show, both with the box model and HAMOCC5, that short term perturbations of river inputs (on time scales of centuries) only significantly affect nutrient cycling on continental margins. The effects of the perturbations of river inputs of nutrients on continental shelves are partly attenuated by nutrient supply from the opean ocean. Only long term perturbations of riverine delivery of nutrients on time scales of kyrs can affect the ocean on a global scale. Future changes of the marine silicon cycle will depend on a multitude of factors including changes in N, P and Si inputs from rivers, global warming (through its effects on solubility of bSiO 2 and oceanic circulation) and river damming.

Conclusions
Despite their different structures, the box model and general circulation model for the marine silica cycle used in this study show surprisingly comparable responses to changes in river input on long and short time scales. Thus, HAMOCC2 and the box model predict a similar export production and opal burial on a time scale of 150 kyrs although the temporal dynamics differ slightly. For simulations of 150 years, the box model and HAMOCC5 forecast comparable decreases in export production and sediment burial on the continental margins and in the open ocean in response to increased Si retention in rivers. Only the amplitude of change in the results of the box model is slightly higher than in the results of the GCM due to the absence of other nutrients in the former model. Results of both models also demonstrate the role of the continental shelf as a major sink of silica at the global scale. Furthermore, coastal waters appear to be more sensitive to perturbations of riverine inputs than the open ocean. Ultimately, the choice of one model over another should be based on the availability of data and technical limitations as much as on model performance. For structural reasons, the box model can not be applied to problems requiring explicit spatial representations. Nevertheless, box models remain suitable tools to evaluate the global scale response of both continental margins and the open ocean to global scale perturbations, especially on longer time scales. Given its low data demand and computational requirements, the box model is the most user friendly of the modeling tools compared here. The performance of our box model could be further improved by coupling it to models for other nutrients, such as N and P (e.g. Wallmann, 2003; Slomp and Van Cappellen, 2007). Such an implementation would potentially show the sensitivity of the marine silica cycle to anthropogenic perturbations of Si:N and Si:P. HAMOCC2 is a well-tested GCM that has been used in many long term studies of biogeochemical cycles of various nutrients (Heinze et al., 2003Heinze and MaierReimer, 1999). The strength of the HAMOCC2 model is that it can be integrated over thousands of years in an acceptable amount time while retaining a relatively high horizontal resolution. Its spatial resolution does however not allow for a full representation of the continental margins, however, and the regional effects of riverine inputs cannot be assessed well. HAMOCC5 is, by far, the most demanding model in terms of computer power and requires high performance parallel clusters to run efficiently. Its spatial resolution is sufficient to spatially resolve current riverine inputs Seitzinger et al., 2005) and paves the way for a detailed assessment of global and regional riverine contributions to marine nutrient cycling. Future work with HAMOCC5 could include spatially explicit scenarios for human-induced alterations in river biogeochemistry. Given the spin-up time required for equilibration with the sediment, its use in global oceanic budgeting is limited.