Reviews and syntheses: The biogeochemical cycle of silicon in the modern ocean

The element silicon (Si) is required for the growth of silicified organisms in marine environments, such as diatoms. These organisms consume vast amounts of Si together with N, P, and C, connecting the biogeochemical cycles of these elements. Thus, understanding the Si cycle in the ocean is critical for understanding wider issues such as carbon sequestration by the ocean's biological pump. In this review, we show that recent advances in process studies indicate that total Si inputs and outputs, to and from the world ocean, are 57 % and 37 % higher, respectively, than previous estimates. We also update the total ocean silicic acid inventory value, which is about 24 % higher than previously estimated. These changes are significant, modifying factors such as the geochemical residence time of Si, which is now about 8000 years, 2 times faster than previously assumed. In addition, we present an updated value of the global annual pelagic biogenic silica production (255 Tmol Si yr⁻¹) based on new data from 49 field studies and 18 model outputs, and we provide a first estimate of the global annual benthic biogenic silica production due to sponges (6 Tmol Si yr⁻¹). Given these important modifications, we hypothesize that the modern ocean Si cycle is at approximately steady state with inputs $=\mathrm{14.8}(\pm \mathrm{2.6})$ Tmol Si yr⁻¹ and outputs $=\mathrm{15.6}(\pm \mathrm{2.4})$ Tmol Si yr⁻¹. Potential impacts of global change on the marine Si cycle are discussed.

1 Introduction

Silicon, the seventh most abundant element in the universe, is the second most abundant element in the Earth's crust. The weathering of the Earth's crust by CO₂-rich rainwater, a key process in the control of atmospheric CO₂ (Berner et al., 1983; Wollast and Mackenzie, 1989), results in the generation of silicic acid (dSi; Si(OH)₄) in aqueous environments. Silicifiers are among the most important aquatic organisms and include micro-organisms (e.g., diatoms, rhizarians, silicoflagellates, several species of choanoflagellates) and macro-organisms (e.g., siliceous sponges). Silicifiers use dSi to precipitate biogenic silica (bSi; SiO₂) as internal (Moriceau et al., 2019) and/or external (Maldonado et al., 2019) structures. Phototrophic silicifiers, such as diatoms, globally consume vast amounts of Si concomitantly with nitrogen (N), phosphorus (P), and inorganic carbon (C), connecting the biogeochemistry of these elements and contributing to the sequestration of atmospheric CO₂ in the ocean (Tréguer and Pondaven, 2000). Heterotrophic organisms like rhizarians, choanoflagellates, and sponges produce bSi independently of the photoautotrophic processing of C and N (Maldonado et al., 2012, 2019; Llopis Monferrer et al., 2020).

Understanding the Si cycle is critical for understanding the functioning of marine food webs, biogeochemical cycles, and the biological carbon pump. Herein, we review recent advances in field observations and modeling that have changed our understanding of the global Si cycle and provide an update of four of the six net annual input fluxes and of all the output fluxes previously estimated by Tréguer and De La Rocha (2013). Taking into account numerous field studies in different marine provinces and model outputs, we re-estimate the Si production (Nelson et al., 1995), review the potential contribution of rhizarians (Llopis Monferrer et al., 2020) and picocyanobacteria (Ohnemus et al., 2016), and give an estimate of the total bSi production by siliceous sponges using recently published data on sponge bSi in marine sediments (Maldonado et al., 2019). We discuss the question of the balance and imbalance of the marine Si biogeochemical cycle at different timescales, and we hypothesize that the modern ocean Si cycle is potentially at steady state with inputs $=\mathrm{14.8}(\pm \mathrm{2.6})$ Tmol Si yr⁻¹ approximately balancing outputs $=\mathrm{15.6}(\pm \mathrm{2.4})$ Tmol Si yr⁻¹ (Fig. 1). Finally, we address the question of the potential impact of anthropogenic activities on the global Si cycle and suggest guidelines for future research endeavors.

2 Advances in input fluxes

Silicic acid is delivered to the ocean through six pathways as illustrated in Fig. 1, which all ultimately derive from the weathering of the Earth's crust (Tréguer and De La Rocha, 2013). All fluxes are given with an error of 1 standard deviation.

Figure 1 Schematic view of the Si cycle in the modern world ocean (input, output, and biological Si fluxes), and possible balance (total Si inputs=total Si outputs=15.6 Tmol Si yr⁻¹) in reasonable agreement with the individual range of each flux (F); see Tables 1 and 2. The white arrows represent fluxes of net sources of silicic acid (dSi) and/or of dissolvable amorphous silica (aSi) and of dSi recycled fluxes. Orange arrows correspond to sink fluxes of Si (either as biogenic silica or as authigenic silica). Green arrows correspond to biological (pelagic) fluxes. All fluxes are in teramoles of silicon per year (Tmol Si yr⁻¹). Details are given in Sect. S1 in the Supplement.

2.1 Riverine (F_R) and aeolian (F_A) contributions

The best estimate for the riverine input (F_R) of dSi, based on data representing 60 % of the world river discharge and a discharge-weighted average dSi riverine concentration of 158 µM−Si (Dürr et al., 2011), remains at $F_{\mathrm{R}\mathrm{dSi}}=\mathrm{6.2}(\pm \mathrm{1.8})$ Tmol Si yr⁻¹ (Tréguer and De La Rocha, 2013). However, not only dSi is transferred from the terrestrial to the riverine system, with particulate Si mobilized in crystallized or amorphous forms (Dürr et al., 2011). According to Saccone et al. (2007), the term “amorphous silica” (aSi) includes biogenic silica (bSi, from phytoliths, freshwater diatoms, sponge spicules), altered bSi, and pedogenic silicates, the three of which can have similar high solubilities and reactivities. Delivery of aSi to the fluvial system has been reviewed by Frings et al. (2016), and they suggested a value of $F_{\mathrm{R}\mathrm{aSi}}=\mathrm{1.9}(\pm \mathrm{1.0})$ Tmol Si yr⁻¹. Therefore, total $F_{\mathrm{R}}=\mathrm{8.1}(\pm \mathrm{2.0})$ Tmol Si yr⁻¹.

No progress has been made regarding aeolian dust deposition into the ocean (Tegen and Kohfeld, 2006) and subsequent release of dSi via dust dissolution in seawater since Tréguer and De La Rocha (2013), who summed the flux of particulate dissolvable silica and wet deposition of dSi through precipitation. Thus, our best estimate for the aeolian flux of dSi, F_A, remains 0.5(±0.5) Tmol Si yr⁻¹.

2.2 Dissolution of minerals (F_W)

As shown in Fig. 2, the low-temperature dissolution of siliceous minerals in seawater and from sediments feeds a dSi flux, F_W, through two processes: (1) the dissolution of river-derived lithogenic particles deposited along the continental margins and shelves and (2) the dissolution of basaltic glass in seawater, processes that work mostly in deep waters. About 15–20 Gt yr⁻¹ of river-derived lithogenic particles are deposited along the margins and shelves (e.g., Syvitskia et al., 2003; also see Fig. 2). Dissolution experiments with river sediments or basaltic glass in seawater showed that 0.08 %–0.17 % of the Si in the solid phase was released within a few days to months (e.g., Oelkers et al., 2011; Jones et al., 2012; Pearce et al., 2013; Morin et al., 2015). However, the high solid-to-solution ratios in these experiments increased the dSi concentration quickly to near-equilibrium conditions inhibiting further dissolution, which prevents direct comparison with natural sediments. Field observations and subsequent modeling of Si release range around 0.5 %–5 % yr⁻¹ of the Si originally present in the solid phase dissolved into the seawater (e.g., Arsouze et al., 2009; Jeandel and Oelkers, 2015). On the global scale, Jeandel et al. (2011) estimated the total flux of dissolution of minerals to range between 0.7–5.4 Tmol Si yr⁻¹, i.e., similar to the dSi river flux. However, this estimate is based on the assumption of 1 %–3 % congruent dissolution of sediments for a large range of lithological composition which, so far, has not been proven.

Figure 2 Schematic view of the low-temperature processes that control the dissolution of (either amorphous or crystallized) siliceous minerals in seawater in and to the coastal zone and in the deep ocean, feeding F_GW and F_W. These processes correspond to both low and medium energy flux dissipated per volume of a given siliceous particle in the coastal zone, in the continental margins, and in the abysses and to high-energy flux dissipated in the surf zone. Details are given in S1 in the Supplement.

Another approach to estimate F_W is to consider the benthic efflux from sediments devoid of biogenic silica deposits. Frings (2017) estimates that “non-biogenic-silica” sediments (i.e., clays and calcareous sediments, which cover about 78 % of the ocean area) may contribute up to 44.9 Tmol Si yr⁻¹ via a benthic diffusive Si flux. However, according to lithological descriptions given in GSA Data Repository 2015271 some of the non-biogenic-silica sediment classes described in this study may contain significant bSi, which might explain this high estimate for F_W. Tréguer and De La Rocha (2013) considered benthic efflux from non-siliceous sediments ranging between ≈10–20 $\mathrm{mmol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ in agreement with Tréguer et al. (1995). If extrapolated to the 120 million square kilometer zone of opal-poor sediments in the global ocean, this gives an estimate of $F_{\mathrm{W}}=\mathrm{1.9}(\pm \mathrm{0.7})$ Tmol Si yr⁻¹.

2.3 Submarine groundwater (F_GW)

Since 2013, several papers have sought to quantify the global oceanic input of dissolved Si (dSi) from submarine groundwater discharge (SGD), which includes terrestrial (freshwater) and marine (saltwater) components (Fig. 2). Silicic acid inputs through SGD may be considerable, similar to or in excess of riverine input in some places. For instance, Georg et al. (2009) estimated this input to be 0.093 Tmol Si yr⁻¹ in the Bay of Bengal, which is ≈66 % of the Ganges–Brahmaputra river flux of dSi to the ocean. At a global scale the best estimate of Tréguer and De La Rocha (2013) for F_GW was 0.6(±0.6) Tmol Si yr⁻¹. More recently, Rahman et al. (2019) used a global terrestrial SGD flux model weighted according to aquifer lithology (Beck et al., 2013) in combination with a compilation of dSi in shallow water coastal aquifers to derive a terrestrial groundwater input of dSi to the world ocean of 0.7(±0.1) Tmol Si yr⁻¹. This new estimate, with its relatively low uncertainty, represents the lower limit flux of dSi to the ocean via SGD. The marine component of SGD, driven by a range of physical processes such as density gradients or waves and tides, is fed by seawater that circulates through coastal aquifers or beaches via advective flow paths (Fig. 2; also see Fig. 1 of Li et al., 1999). This circulating seawater may become enriched in dSi through bSi or mineral dissolution, the degree of enrichment being determined by subsurface residence time and mineral type (Techer et al., 2001; Anschutz et al., 2009; Ehlert et al. 2016a).

Several lines of evidence show that the mineral dissolution (strictly corresponding to net dSi input) may be substantial (e.g., Ehlert et al., 2016b). Focusing on processes occurring in tidal sands, Anschultz et al. (2009) showed that they can be a biogeochemical reactor for the Si cycle. Extrapolating laboratory-based dissolution experiments performed with pure quartz, Fabre et al. (2019) calculated the potential flux of dissolution of siliceous sandy beaches that is driven by wave and tidal action. If, according to Luijendijk et al. (2018) one-third of the world's shorelines are sandy beaches, this dissolution flux could be 3.2(±1.0) Tmol Si yr⁻¹. However, this estimate is not well constrained because it has not been validated by field experiments (Sect. S2 in the Supplement). Cho et al. (2018), using a ²²⁸Ra inverse model and groundwater $\mathrm{dSi}{/}^{\mathrm{228}}\mathrm{Ra}$ ratios, estimate the total (terrestrial + marine) SGD dSi flux to the ocean to be 3.8(±1.0) Tmol Si yr⁻¹; this represents an upper limit value for SGD's contribution to the global ocean dSi cycle. Without systematic data that corroborate the net input of dSi through the circulation of the marine component of SGD (e.g., porewater δ³⁰Si, paired dSi and ²²⁸Ra measurements), we estimate the range of net input of dSi through total SGD as 0.7 Tmol Si yr⁻¹ (Rahman et al., 2019) to 3.8 Tmol Si yr⁻¹ (Cho et al., 2018), with an average, i.e., $F_{\text{GW}}=\mathrm{2.3}(\pm \mathrm{1.1})$ Tmol Si yr⁻¹, which is approximately 3 times larger than that of Tréguer and De La Rocha (2013).

2.4 (Sub)polar glaciers (F_ISMW)

This flux was not considered by Tréguer and De La Rocha (2013). Several researchers have now identified polar glaciers as sources of Si to marine environments (Tréguer, 2014; Meire et al., 2016; Hawkings et al., 2017). The current best estimate of discharge-weighted dSi concentration in (sub)Arctic glacial meltwater rivers lies between 20–30 µM, although concentrations ranging between 3 and 425 µM have been reported (Graly et al., 2014; Meire et al., 2016; Hatton et al., 2019). Only two values currently exist for dSi from subglacial meltwater beneath the Antarctic Ice Sheet (Whillans subglacial lake and Mercer subglacial lake, 126–140 µM; Michaud et al., 2016, Hawkings et al., 2020), and there is a limited dataset from periphery glaciers in the McMurdo Dry Valleys and Antarctic Peninsula (≈10–120 µM; Hatton et al., 2020; Hirst et al., 2020). Furthermore, iceberg dSi concentrations remain poorly quantified but are expected to be low (≈5 µM; Meire et al., 2016). Meltwater typically contains high suspended sediment concentrations, due to intense physical erosion by glaciers, with a relatively high dissolvable aSi component (0.3 %–1.5 % dry weight) equating to concentrations of 70–340 µM (Hawkings et al., 2018; Hatton et al., 2019). Iceberg aSi concentrations are lower (28–83 µM) (Hawkings et al., 2017). This particulate phase appears fairly soluble in seawater (Hawkings et al., 2017), and large benthic dSi fluxes in glacially influenced shelf seas have been observed (Hendry et al., 2019; Ng et al., 2020). Direct silicic acid input from (sub)polar glaciers is estimated to be 0.04(±0.04) Tmol Si yr⁻¹. If the aSi flux is considered then this may provide an additional 0.29(±0.22) Tmol Si yr⁻¹, with a total F_ISMW ($=\mathrm{dSi}+\mathrm{aSi}$) input estimate of 0.33(±0.26) Tmol Si yr⁻¹. This does not include any additional flux from benthic processing of glacially derived particles in the coastal regions (see Sect. 2.2 above).

2.5 Hydrothermal activity (F_H)

The estimate of Tréguer and De La Rocha (2013) for F_H was 0.6(±0.4) Tmol Si yr⁻¹. Seafloor hydrothermal activity at mid-ocean ridges (MORs) and ridge flanks is one of the fundamental processes controlling the exchange of heat and chemical species between seawater and ocean crust (Wheat and Mottl, 2000). A major challenge limiting our current models of both heat and mass flux (e.g., Si flux) through the seafloor is estimating the distribution of the various forms of hydrothermal fluxes, including focused (i.e., high temperature) vs. diffuse (i.e., low temperature) and ridge axis vs. ridge flank fluxes. Estimates of the Si flux for each input are detailed below.

Axial and near-axial hydrothermal flux settings. The best estimate of the heat flux at ridge axis (i.e., crust 0–0.1 Ma in age) is 1.8(±0.4) TW, while the heat flux in the near-axial region (i.e., crust 0.1–1 Ma in age) has been inferred at 1.0(±0.5) TW (Mottl, 2003). The conversion of heat flux to hydrothermal water and chemical fluxes requires assumptions regarding the temperature at which this heat is removed. For an exit temperature of 350(±30) ^∘C typical of black smoker vent fluids, and an associated enthalpy of 1500(±190) J g⁻¹ at 450–1000 bar and heat flux of 2.8(±0.4) TW, the required seawater flux is $\mathrm{5.9}(\pm \mathrm{0.8})\times {\mathrm{10}}^{\mathrm{16}}$ g yr⁻¹ (Mottl, 2003). High-temperature hydrothermal dSi flux is calculated using a dSi concentration of 19(±11) mmol kg⁻¹, which is the average concentration in hydrothermal vent fluids that have an exit temperature>300 ^∘C (Mottl, 2012). This estimate is based on a compilation of >100 discrete vent fluid data, corrected for seawater mixing (i.e., end-member values at Mg=0; Edmond et al., 1979) and phase separation. Although the chlorinity of hot springs varies widely, nearly all of the reacted fluid, whether vapor or brine, must eventually exit the crust within the axial region. The integrated hot spring flux must therefore have a chlorinity similar to that of seawater. The relatively large range of dSi concentrations in high-temperature hydrothermal fluids likely reflect the range of geological settings (e.g., fast- and slow-spreading ridges) and host-rock composition (ultramafic, basaltic, or felsic rocks). Because dSi enrichment in hydrothermal fluids results from mineral–fluid interactions at depth and is mainly controlled by the solubility of secondary minerals such as quartz (Mottl 1983; Von Damm et al. 1991), it is also possible to obtain a theoretical estimate of the concentration of dSi in global hydrothermal vent fluids. Under the conditions of temperature and pressure (i.e., depth) corresponding to the base of the upflow zone of high temperature (>350–450 ^∘C) hydrothermal systems, dSi concentrations between 16 and 22 mmol kg⁻¹ are calculated, which is in good agreement with measured values in end-member hydrothermal fluids. Using a dSi concentration of 19(±3.5) mmol kg⁻¹ and water flux of $\mathrm{4.8}(\pm \mathrm{0.8})\times {\mathrm{10}}^{\mathrm{16}}$ g yr⁻¹, we determine an axial hydrothermal Si flux of 0.91(±0.29) Tmol Si yr⁻¹. It should be noted, however, that high-temperature hydrothermal fluids may not be entirely responsible for the transport of all the axial hydrothermal heat flux (Elderfield and Schultz, 1996; Nielsen et al., 2006). Because dSi concentrations in diffuse hydrothermal fluids are not significantly affected by subsurface Si precipitation during cooling of the hydrothermal fluid (Escoube et al., 2015), we propose that the global hydrothermal Si flux is not strongly controlled by the nature (focused vs. diffuse) of axial fluid flow.

Ridge flank hydrothermal fluxes. Chemical fluxes related to seawater–crust exchange at ridge flanks have been previously determined through direct monitoring of fluids from low-temperature hydrothermal circulation (Wheat and Mottl, 2000). Using basaltic formation fluids from the 3.5 Ma crust on the eastern flank of the Juan de Fuca Ridge (Wheat and McManus, 2005), a global flux of 0.011 Tmol Si yr⁻¹ for the warm ridge flank is calculated. This estimate is based on the measured Si anomaly associated with warm spring (0.17 mmol kg⁻¹) and a ridge flank fluid flux determined using oceanic Mg mass balance, therefore assuming that the ocean is at steady state with respect to Mg. More recent results of basement fluid compositions in cold and oxygenated ridge flank settings (e.g., North Pond, Mid-Atlantic Ridge) also confirm that incipient alteration of volcanic rocks may result in significant release of Si to circulating seawater (Meyer et al., 2016). The total heat flux through ridge flanks, from 1 Ma crust to a sealing age of 65 Ma, has been estimated at 7.1(±2) TW. Considering that most of ridge flank hydrothermal power output should occur at cool sites (<20 ^∘C), the flux of slightly altered seawater could range from 0.2 to 2×10¹⁹ g yr⁻¹, rivaling the flux of river water to the ocean of 3.8×10¹⁹ g yr⁻¹ (Mottl, 2003). Using this estimate and Si anomaly of 0.07 $\mathrm{mmol}-\mathrm{Si}{\mathrm{kg}}^{-\mathrm{1}}$ reported in cold ridge flank setting from North Pond (Meyer et al., 2016), a Si flux of 0.14 to 1.4 Tmol Si yr⁻¹ for the cold ridge flank could be calculated. Because of the large volume of seawater interacting with oceanic basalts in ridge flank settings, even a small chemical anomaly resulting from reactions within these cold systems could result in a globally significant elemental flux. Hence, additional studies are required to better determine the importance of ridge flanks to oceanic Si budget.

Combining axial and ridge flank estimates, the best estimate for F_H is now 1.7(±0.8) Tmol Si yr⁻¹, approximately 3 times larger than the estimate from Tréguer and De La Rocha (2013).

2.6 Total net inputs (Table 1A)

Total Si input $=\mathrm{8.1}(\pm \mathrm{2.0})\left(F_{\mathrm{R}(\mathrm{dSi}+\mathrm{aSi})}\right)+\mathrm{0.5}(\pm \mathrm{0.5})\left(F_{\mathrm{A}}\right)+\mathrm{1.9}(\pm \mathrm{0.7})\left(F_{\mathrm{W}}\right)+\mathrm{2.3}(\pm \mathrm{1.1})\left(F_{\text{GW}}\right)+\mathrm{0.3}(\pm \mathrm{0.3})\left(F_{\text{ISMW}}\right)+\mathrm{1.7}(\pm \mathrm{0.8})\left(F_{\mathrm{H}}\right)=\mathrm{14.8}(\pm \mathrm{2.6})$ Tmol Si yr⁻¹.

The uncertainty of the total Si inputs (and total Si outputs, Sect. 3) has been calculated using the error propagation method from Bevington and Robinson (2003). This has been done for both the total fluxes and the individual flux estimates.

3 Advances in output fluxes

3.1 Long-term burial of planktonic biogenic silica in sediments (F_B)

Long-term burial of bSi, which generally occurs below the top 10–20 cm of sediment, was estimated by Tréguer and De La Rocha (2013) to be 6.3(±3.6) Tmol Si yr⁻¹. The burial rates are highest in the Southern Ocean (SO), the North Pacific Ocean, the equatorial Pacific Ocean, and the coastal and continental margin zone (CCMZ; DeMaster, 2002; Hou et al., 2019; Rahman et al., 2017).

Post-depositional redistribution by processes like winnowing or focusing by bottom currents can lead to under- and over-estimation of uncorrected sedimentation and burial rates. To correct for these processes, the burial rates are typically normalized using the particle reactive nuclide ²³⁰Th method (e.g., Geibert et al., 2005). A ²³⁰Th normalization of bSi burial rates has been extensively used for the SO (Tréguer and De La Rocha, 2013), particularly in the “opal belt” zone (Pondaven et al., 2000; DeMaster, 2002; Geibert et al., 2005). Chase et al. (2015) re-estimated the SO burial flux, south of 40^∘ S at 2.3(±1.0) Tmol Si yr⁻¹.

Hayes et al. (2020, 2021) recently calculated total marine bSi burial of 5.46(±1.18) Tmol Si yr⁻¹, using a database that comprises 2948 bSi concentrations of top core sediments and ²³⁰Th-corrected accumulation fluxes of open-ocean locations >1 km in depth. The ²³⁰Th-corrected total burial rate of Hayes et al. (2021) is 2.68(±0.61) Tmol Si yr⁻¹ south of 40^∘ S, close to the estimate of Chase et al. (2015) for the SO. Hayes et al. (2021) do not distinguish between the different analytical methods used for the determination of the bSi concentrations of these 2948 samples to calculate total bSi burial. These methods include alkaline digestion methods (with variable protocols for correcting from lithogenic interferences; e.g., DeMaster, 1981; Mortlock and Froelich, 1989; Müller and Schneider, 1993), X-ray diffraction (e.g., Leinen et al., 1986), X-ray fluorescence (e.g., Finney et al., 1988), Fourier-transform infrared spectroscopy (Lippold et al., 2012), and inductively coupled plasma mass spectrometry (e.g., Prakash Babu et al., 2002). An international exercise calibration on the determination of bSi concentrations of various sediments (Conley, 1998) concluded that the X-ray diffraction (XRD) method generated bSi concentrations that were on average 24 % higher than the alkaline digestion methods. In order to test the influence of the XRD method on their re-estimate of total bSi burial, Hayes et al. (2021) found that their re-estimate (5.46(±1.18) Tmol Si yr⁻¹), which includes XRD data (≈40 % of the total number of data points), did not differ significantly from a re-estimate that does not include XRD data points (5.43(±1.18) Tmol Si yr⁻¹). As a result, this review includes the re-estimate of Hayes et al. (2021) for the open-ocean annual burial rate, i.e., 5.5(±1.2) Tmol Si yr⁻¹.

The best estimate for the open-ocean total burial now becomes 2.8(±0.6) Tmol Si yr⁻¹ without the SO contribution (2.7(±0.6) Tmol Si yr⁻¹). This value is an excess of 1.8 Tmol Si yr⁻¹ over the DeMaster (2002) and Tréguer and De La Rocha (2013) estimates, which were based on 31 sediment cores mainly distributed in the Bering Sea, the North Pacific, the Sea of Okhotsk, and the equatorial Pacific (total area 23 million square kilometers) and where bSi % was determined solely using alkaline digestion methods.

Estimates of the silica burial rates have usually been determined from carbon burial rates using a Si:C ratio of 0.6 in CCMZ (DeMaster 2002). However, we now have independent estimates of marine organic C and total initial bSi burial (e.g., Aller et al., 1996, 2008; Galy et al., 2007; Rahman et al., 2016, 2017). It has been shown that the initial bSi burial in sediment evolved as unaltered bSi or as authigenically formed aluminosilicate phases (Rahman et al., 2017). The Si:C burial ratios of residual marine plankton post-remineralization in tropical and subtropical deltaic systems are much greater (2.4–11) than the 0.6 Si:C burial ratio assumed for continental margin deposits (DeMaster, 2002). The sedimentary Si:C preservation ratios are therefore suggested to depend on differential remineralization pathways of marine bSi and C_org under different diagenetic regimes (Aller, 2014). Partitioning of ³²Si activities between bSi and mineral pools in tropical deltaic sediments indicate rapid and near-complete transformation of initially deposited bSi to authigenic clay phases (Rahman et al., 2017). For example, in subtropical/temperate deltaic and estuarine deposits, ³²Si activities represent approximately ≈50 % of initial bSi_opal delivery to sediments (Rahman et al., 2017). Using the ³²Si technique, Rahman et al. (2017) provided an updated estimate of bSi burial for the CCMZ of 3.7(±2.1) Tmol Si yr⁻¹, higher than the Tréguer and De La Rocha (2013) estimate of 3.3(±2.1) Tmol Si yr⁻¹ based on the Si:C method of DeMaster (2002).

Combining the Hayes et al. (2021) burial rate for the open-ocean zone including the SO and the Rahman et al. (2017) estimate for the CCMZ gives a revised global total burial flux, F_B, of 9.2(±1.6) Tmol Si yr⁻¹, 46 % larger than the Tréguer and De La Rocha (2013) estimate.

3.2 Deposition and long-term burial of sponge silica (F_SP)

The estimate of Tréguer and De La Rocha (2013) for F_SP, the net sink of sponge bSi in sediments of continental margins, was 3.6(±3.7) Tmol Si yr⁻¹. The longevity of sponges, ranging from years to millennia, temporally decouples the process of skeleton production from the process of deposition to the sediments (Jochum et al., 2017). While sponges slowly accumulate bSi over their long and variable lifetimes (depending on the species), the deposition to the sediments of the accumulated bSi is a relatively rapid process after sponge death, lasting days to months (Sect. S3 in the Supplement). The estimate of Tréguer and De La Rocha (2013) was calculated as the difference between the sponge dSi demand on continental shelves (3.7(±3.6) Tmol Si yr⁻¹) – estimated from silicon consumption rates available for a few sublittoral sponge species (Maldonado et al., 2011) – and the flux of dSi from the dissolution of sponge skeletons in continental shelves (0.15(±0.15) Tmol Si yr⁻¹). This flux was tentatively estimated from the rate of dSi dissolution from a rare, unique glass sponge reef in British Columbia (Canada; Chu et al., 2011) and which is unlikely to be representative of the portion of sponge bSi that dissolves back as dSi after sponge death and before their burial in the sediments. To improve the estimate, Maldonado et al. (2019) used microscopy to access the amount of sponge silica that was actually being buried in the marine sediments using 17 sediment cores representing different marine environments. The deposition of sponge bSi was found to be 1 order of magnitude more intense in sediments of continental margins and seamounts than on continental rises and central basin bottoms. The new best estimate for F_SP is 1.7(±1.6) Tmol Si yr⁻¹, assuming that the rate of sponge silica deposition in each core was approximately constant through the Holocene, i.e., 2 times smaller than Tréguer and De La Rocha's preliminary estimate.

3.3 Reverse weathering flux (F_RW)

The previous estimate for this output flux, provided by Tréguer and De La Rocha (2013), $F_{\text{RW}}=\mathrm{1.5}(\pm \mathrm{0.5})$ Tmol Si yr⁻¹, was determined using indirect evidence since the influence of reverse weathering on the global Si cycle prior to 2013 was poorly understood. For example, reverse weathering reactions at the sediment–water interface were previously thought to constitute a relatively minor sink (0.03–0.6 Tmol Si yr⁻¹) of silica in the ocean (DeMaster, 1981). The transformation of bSi to a neoformed aluminosilicate phase, or authigenic clay formation, was assumed to proceed slowly (>10⁴–10⁵ years) owing principally to the difficulty of distinguishing the contribution of background lithogenic or detrital clays using the common leachates employed to quantify bSi (DeMaster, 1981). Recent direct evidence supporting the rapid formation of authigenic clays comes from tropical and subtropical deltas (Michalopoulos and Aller, 1995; Rahman et al., 2016, 2017; Zhao et al., 2017), and several geochemical tools show that authigenic clays may form ubiquitously in the global ocean (Michalopoulos and Aller, 2004; Ehlert et al., 2016a; Baronas et al., 2017; Geilert et al., 2020; Pickering et al., 2020). Activities of cosmogenic ³²Si ($t_{\mathrm{1}/\mathrm{2}}\approx \mathrm{140}$ years), incorporated into bSi in the surface ocean, provide demonstrable proof of rapid reverse weathering reactions by tracking the fate of bSi upon delivery to marine sediments (Rahman et al., 2016). By differentiating sedimentary bSi storage between unaltered bSi (bSi_opal) and diagenetically altered bSi (bSi_altered) in the proximal coastal zone, ³²Si activities in these pools indicate that 3.7 Tmol Si yr⁻¹ is buried as unaltered bSi_opal and 4.7(±2.3) Tmol Si yr⁻¹ as authigenic clays (bSi_clay) on a global scale. Here, we adopt 4.7 Tmol Si yr⁻¹ for F_RW representing about 3 times the value of Tréguer and De La Rocha (2013).

3.4 Total net output (Table 1A)

Total Si output = 9.2 (± 1.6) (F_B(net deposit)) + 4.7 (± 2.3) (F_RW) + 1.7 (± 1.6) (F_SP) = 15.6 (± 2.4) Tmol Si yr⁻¹.

4 Advances in biological fluxes

4.1 Annual bSi pelagic production

4.1.1 From field data

The last evaluation of global marine silica production was by Nelson et al. (1995), who estimated global gross marine bSi pelagic production to be 240(±40) Tmol Si yr⁻¹. Since 1995, the number of field studies of bSi production (using either the ³⁰Si tracer method, Nelson and Goering (1977), or the ³²Si method (Tréguer et al., 1991; Brzezinski and Phillips, 1997)) has grown substantially from 15 (1995) to 49 in 2019, allowing the first estimate based on empirical silica production rate measurements (Fig. 3 and Sect. S4 in the Supplement). It is usually assumed that the silica production, as measured using the above methods, is mostly supported by diatoms, with some unknown (but minor) contribution of other planktonic species.

Figure 3 Biogenic silica production measurements in the world ocean. Distribution of stations in the Longhurst biogeochemical provinces (Longhurst, 2007; Longhurst et al., 1995). All data are shown in Sect. S4 in the Supplement (Annex 1).

The silica production rates measured during 49 field campaigns were assigned to Longhurst provinces (Longhurst, 2007; Longhurst et al., 1995) based on location, with the exception of the SO, where province boundaries were defined according to Tréguer and Jacques (1992). Extrapolating these “time-and-space-limited” measurements of bSi spatially to a biogeographic province, and annually from the bloom phenology for each province (calculated as the number of days where the chlorophyll concentration is greater than the average concentration between the maximum and the minimum values), results in annual silica production estimates for 26 of the 56 world ocean provinces. The annual production of all provinces in a basin were averaged for the “ocean basin” estimate (Table 2) and then extrapolated by basin area. The averages from provinces were subdivided among coastal for the “domain” estimate (Table 2), SO, and open-ocean domains and extrapolated based on the area of each domain. Averaging the ocean basin and the domain estimates (Table 2), our best estimate for the annual global marine bSi production is 267(±18) Tmol Si yr⁻¹ (Table 2).

Table 1 Si inputs, outputs, and biological fluxes at word ocean scale.

References are (1) Nelson et al. (1995) and (2) Tréguer and De La Rocha (2013). (3) Recalculated from our updated dSi inventory value; see Supplement for detailed definition of flux term (in detailed legend of Fig. 1).

Table 2 Biological fluxes (F_Pgross in Tmol Si yr⁻¹).

Global silica production as determined from numerical models and extrapolated from field measurements of silica production (uncertainties are standard errors).

4.1.2 Annual bSi pelagic production from models

Estimates of bSi production were also derived from satellite productivity models and from global ocean biogeochemical models (GOBMs). We used global net primary production (NPP) estimates from the carbon-based productivity model (Westberry et al., 2008) and the vertically generalized productivity model (VGPM) (Behrenfeld and Falkowski, 1997) for the estimates based on satellite productivity models. NPP estimates from these models were divided into oligotrophic (<0.1 µg chl a L⁻¹), mesotrophic (0.1–1.0 µg chl a L⁻¹), and eutrophic (>1.0 µg chl a L⁻¹) areas (Carr et al., 2006). The fraction of productivity by diatoms in each area was determined using the DARWIN model (Dutkiewicz et al., 2015), allowing a global estimate where diatoms account for 29 % of the production. Each category was further subdivided into high-nutrient low-chlorophyll (HNLC) zones (>5 µM surface nitrate; Garcia et al., 2014), coastal zones (<300 km from a coastline), and open-ocean (remainder) zones for application of Si:C ratios to convert to diatom silica production. Si:C ratios were 0.52 for HNLC regions, 0.065 for the open ocean, and 0.13 for the coastal regions, reflecting the effect of Fe limitation in HNLC areas (Franck et al., 2000), of Si limitation for uptake in the open ocean (Brzezinski et al., 1998, 2011; Brzezinski and Nelson, 1996; Krause et al., 2012), and of replete conditions in the coastal zone (Brzezinski, 1985). Silica production estimates were then subdivided between coast (within 300 km of shore), open ocean, and SO (northern boundary 43^∘ S from Australia to South America, 34.8^∘ S from South America to Australia) and summed to produce regional estimates (Table 2). Our best estimate for the global marine bSi production is 207(±23) Tmol Si yr⁻¹ from satellite productivity models (Table 2).

A second model-based estimate of silica production used 18 numerical GOBMs models of the marine silica cycle that all estimated global silica export from the surface ocean (Gnanadesikian and Toggweiler, 1999; Usbeck, 1999; Heinze et al., 2003; Wischmeyer et al., 2003; Jin et al., 2006; Dunne et al., 2007; Sarmiento et al., 2007; Bernard et al., 2011; Ward et al., 2012; Matsumoto et al., 2013; De Souza et al. 2014; Holzer et al., 2014; Aumont et al., 2015; Dutkiewicz et al., 2015; Pasquier and Holzer, 2017; Roshan et al., 2018). These include variants of the MOM, HAMOCC OCIM, DARWIN, cGENIE, and PICES models. Export production was converted to gross silica production by using a silica dissolution-to-production (D:P) ratio for the surface open ocean of 0.58 and 0.51 for the surface of coastal regions (Tréguer and De La Rocha, 2013). Model results were first averaged within variants of the same model and then averaged across models to eliminate biassing the average to any particular model. Our best estimate from GOBMs for the global marine bSi production is 276(±23) Tmol Si yr⁻¹ (Table 2). Averaging the estimates calculated from satellite productivity models and GOBMs gives a value of 242(±49) Tmol Si yr⁻¹ for the global marine bSi production (Table 2).

4.1.3 Best estimate for annual bSi pelagic production

Using a simple average of the “field” and “model” estimates, the revised best estimate of global marine gross bSi production, mostly due to diatoms, is now $F_{\text{Pgross}}=\mathrm{255}(\pm \mathrm{52})$ Tmol Si yr⁻¹, not significantly different from the Nelson et al. (1995) value.

In the SO, a key area for the world ocean Si cycle (DeMaster, 1981), there is some disagreement among the different methods of estimating bSi production. Field studies give an estimate of 67 Tmol Si yr⁻¹ for the annual gross production of silica in the SO, close to the estimate of 60 Tmol Si yr⁻¹ calculated using satellite productivity models (Table 2). However, the bSi production in the SO estimated by ocean biogeochemical models is about twice as high, at 129 Tmol Si yr⁻¹ (Table 2). The existing in situ bSi production estimates are too sparse to be able to definitively settle whether the lower estimate or the higher estimate is correct, but there is reason to believe that there are potential biases in both the satellite NPP models and the ocean biogeochemical models. SO chlorophyll concentrations may be underestimated by as much as a factor of 3–4 (Johnson et al., 2013), which affects the NPP estimates in this region and hence our bSi production estimates with this method. The bSi production estimated by ocean biogeochemical models is highly sensitive to vertical exchange rates in the SO (Gnanadesikan and Toggweiler, 1999) and is also dependent on the representation of phytoplankton classes in models with explicit representation of phytoplankton. Models that have excessive vertical exchange in the SO (Gnanadesikan and Toggweiler, 1999), or that represent all large phytoplankton as diatoms, may overestimate the Si uptake by plankton in the SO. Other sources of uncertainty in our bSi production estimates include poorly constrained estimates of the Si:C ratio and dissolution-to-production ratios (see Sect. S4 in the Supplement). The errors incurred by these choices are more likely to cancel out in the global average but could be significant at regional scales, potentially contributing to the discrepancies in SO productivity across the various methods.

4.1.4 Estimates of the bSi production of other pelagic organisms

Extrapolations from field and laboratory work show that the contribution of picocyanobacteria (like Synechococcus; Baines et al. 2012, Brzezinski et al., 2017; Krause et al., 2017) to the world ocean accumulation of bSi is <20 Tmol Si yr⁻¹. The gross silica production of rhizarians, siliceous protists, in the 0–1000 m layer might range between 2–58 Tmol Si yr⁻¹, about 50 % of it occurring in the 0–200 m layer (Llopis Monferrer et al., 2020).

Note that these preliminary estimates of bSi accumulation or production by picocyanobacteria and rhizarians are within the uncertainty of our best estimate of F_Pgross.

4.2 Estimates of the bSi production of benthic organisms

The above-updated estimate of the pelagic production does not take into account bSi production by benthic organisms like benthic diatoms and sponges. Our knowledge of the production terms for benthic diatoms is poor, and no robust estimate is available for bSi annual production of benthic diatoms at a global scale (Sect. S4 in the Supplement).

Substantial progress has been made for silica deposition by siliceous sponges recently. Laboratory and field studies reveal that sponges are highly inefficient in the molecular transport of dSi compared to diatoms and consequently bSi production, particularly when dSi concentrations are lower than 75 µM, a situation that applies to most ocean areas (Maldonado et al., 2020). On average, sponge communities are known to produce bSi at rates that are about 2 orders of magnitude smaller than those measured for diatom communities (Maldonado et al., 2012). The global standing crop of sponges is very difficult to be constrained. The annual bSi production attained by such standing crop is even more difficult to estimate because sponge populations are not homogeneously distributed on the marine benthic environment, and extensive, poorly mapped, and unquantified aggregations of heavily silicified sponges occur in the deep sea of all oceans. A first tentative estimate of bSi production for sponges on continental shelves, where sponge biomass can be more easily approximated, ranged widely, from 0.87 to 7.39 Tmol Si yr⁻¹, because of persisting uncertainties in estimating sponge standing crop (Maldonado et al., 2012). A way to estimate the global annual bSi production by sponges without knowing their standing crop is to retrace bSi production values from the amount of sponge bSi that is annually being deposited to the ocean bottom, after assuming that, in the long run, the standing crop of sponges in the ocean is in equilibrium (i.e it is neither progressively increasing nor decreasing over time). The deposition rate of sponge bSi has been estimated at 49.95(±74.14) $\mathrm{mmol}-\mathrm{Si}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ on continental margins, at 0.44(±0.37) $\mathrm{mmol}-\mathrm{Si}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ in sediments of ocean basins where sponge aggregations do not occur, and at 127.30(±105.69) $\mathrm{mmol}-\mathrm{Si}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ in deep-water sponge aggregations (Maldonado et al., 2019). A corrected sponge bSi deposition rate for ocean basins is estimated at 2.98(±1.86) $\mathrm{mmol}-\mathrm{Si}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$ assuming that sponge aggregations do not occupy more than 2 % of seafloor of ocean basins (Maldonado et al., 2019). A total value of 6.15(±5.86) Tmol Si yr⁻¹ can be estimated for the global ocean when the average sponge bSi deposition rate for continental margins and seamounts (representing 108.02 million squarekilometers of seafloor) and for ocean basins (253.86 million squarekilometers) is scaled up through the extension of those bottom compartments. If the bSi production being accumulated as standing stock in the living sponge populations annually is assumed to become constant in a long-term equilibrium state, the global annual deposition rate of sponge bSi can be considered a reliable estimate of the minimum value that the annual bSi production by the sponges can reach in the global ocean. The large associated SD value does not derive from the approach being unreliable but from the spatial distribution of the sponges on the marine bottom being extremely heterogeneous, with some ocean areas being very rich in sponges and sponge bSi in sediments at different spatial scales while other areas are completely deprived of these organisms.

5 Discussion

5.1 Overall residence times

The overall geological residence time for Si in the ocean (τ_G) is equal to the total amount of dSi in the ocean divided by the net input (or output) flux. We re-estimate the total ocean dSi inventory value derived from the Pandora model (Peng et al. 1993), which according to Tréguer et al. (1995) was 97 000 Tmol Si. An updated estimate of the global marine dSi inventory was computed by interpolating the objectively analyzed annual mean silicate concentrations from the 2018 World Ocean Atlas (Garcia et al., 2019) to the OCIM model grid (Roshan et al., 2018). Our estimate is now 120 000 Tmol Si, i.e., about 24 % higher than the Tréguer et al. (1995) estimate. Tables 1B and 3 show updated estimates of τ_G from Tréguer et al. (1995) and Tréguer and De La Rocha (2013) using this updated estimate of the total dSi inventory. Our updated budget (Fig. 1, Tables 1B and 3A) reduces past estimates of τ_G (Tréguer et al., 1995; Tréguer and De La Rocha, 2013) by more than half, from ca. 18 kyr to ca. 8 kyr (Table 3C). This brings the ocean residence time of Si closer to that of nitrogen (<3 kyr; Sarmiento and Gruber, 2006) than phosphorus (30–50 kyr; Sarmiento and Gruber, 2006).

Table 3 A total of 25 years of evolution of the estimates for Si inputs, outputs, biological production, and residence times at World Ocean scale.

References are as follows. (1) Tréguer et al. (1995). (2) Tréguer and De La Rocha (2013). (3) This review. (4) Nelson et al. (1995). ⁽⁵⁾ Recalculated from our updated dSi inventory value.

The overall biological residence time, τ_B, is calculated by dividing the total dSi content of the world ocean by gross silica production. It is calculated from the bSi pelagic production only given the large uncertainty on our estimate of the bSi production by sponges. τ_B is ca. 470 years (Tables 1B and 3). Thus, Si delivered to the ocean passes through the biological uptake and dissolution cycle on average 16 times (τ_G/τ_B) before being removed to the sea floor (Tables 1B and 3C).

The new estimate for the global average preservation efficiency of bSi buried in sediments is ($F_{\mathrm{B}}/F_{\text{Pgross}}=(\mathrm{9.2}/\mathrm{255})=)\mathrm{3.6}$ %, which is similar to the Tréguer and De La Rocha (2013) estimate. This makes bSi in sediments an intriguing potential proxy for export production (Tréguer et al., 2018). Note that the reverse weathering flux (F_RW) is also fed by the export flux (F_E) (Fig. 4). So, the preservation ratio of biogenic silica in sediment can be calculated as $(F_{\mathrm{B}}+F_{\text{RW}})/F_{\text{Pgross}}=(\mathrm{13.9}/\mathrm{255})=\mathrm{5.45}$ %, which is ≈30 times larger than the carbon preservation efficiency.

Figure 4 Schematic view of the Si cycle in the coastal and continental margin zone (CCMZ), linked to the rest of the world ocean (open-ocean zone, including upwelling and polar zones). In this steady-state scenario, consistent with Fig. 1, $\text{total inputs}=\text{total outputs}=\mathrm{15.6}$ Tmol Si yr⁻¹. This figure illustrates the links between biological, burial, and reverse weathering fluxes. It also shows that the open-ocean bSi (pelagic) production (F_P(gross)=222 Tmol Si yr⁻¹) is mostly fueled by dSi inputs from below 92.5 Tmol Si yr⁻¹, with the CCMZ only providing 4.7 Tmol Si yr⁻¹ to the open ocean.

5.2 The issue of steady state

Over a given timescale, an elemental cycle is at steady state if the outputs balance the inputs in the ocean and the mean concentration of the dissolved element remains constant.

5.2.1 Long timescales (>τ_G)

Over geologic timescales, the average dSi concentration of the ocean has undergone drastic changes. A seminal work (Siever, 1991) on the biological–geochemical interplay of the Si cycle showed a decline of a factor of 100 in ocean dSi concentration from 550 Myr ago to the present. This decline was marked by the rise of silicifiers like radiolarian and sponges during the Phanerozoic. Then during the mid-Cenozoic diatoms started to dominate a Si cycle previously controlled by inorganic and diagenetic processes. Conley et al. (2017) hypothesized that biological processes might also have influenced the dSi concentration of the ocean at the start of oxygenic photosynthesis, taking into account the impact of the evolution of biosilicifying organisms (including bacteria-related metabolism). There is further evidence that the existing lineages of sponges have their origin in ancient (Mesozoic) oceans with much higher dSi concentrations than the modern ocean. Some recent sponge species can only complete their silica skeletons if dSi concentration much higher than that in their natural habitat is provided experimentally (Maldonado et al., 1999). Also, all recent sponge species investigated to date have kinetics of dSi consumption that reach their maximum speed only at dSi concentrations that are 1 to 2 orders of magnitude higher than the current dSi availability in the sponge habitats, indicating that the sponge physiology evolved in dSi-richer ancestral scenarios. Note that with a geological residence time of Si of ca. 8000 years, the Si cycle can fluctuate over glacial–interglacial timescales.

5.2.2 Short timescales (<τ_G)

In the modern ocean the main control over silica burial and authigenic formation rate is the bSi production rate of pelagic and benthic silicifiers, as shown above. The gross production of bSi due to diatoms depends on the dSi availability in the surface layer (Fig. 1). Silicic acid does not appear to be limiting in several zones of the world ocean, which include the coastal zones and the HNLC zones (Tréguer and De La Rocha, 2013). Note that any short-term change of dSi inputs does not imply modification of bSi production, or export or burial rate. For this reason, climatic changes or anthropogenic impacts that affect dSi inputs to the ocean by rivers and/or other pathways could lead to an imbalance of Si inputs and outputs in the modern ocean.

5.2.3 A possible steady-state scenario

Within the limits of uncertainty, the total net inputs of dSi and aSi are 14.8(±2.6) Tmol Si yr⁻¹ and are approximately balanced by the total net output flux of Si of 15.6(±2.4) Tmol Si yr⁻¹. Figure 1 supports the hypothesis that the modern ocean Si cycle is at steady state, compatible with the geochemical and biological fluxes of Table 1.

Consistent with Fig. 1, Fig. 4 shows a steady-state scenario for the Si cycle in the coastal and continental margin zone (CCMZ), often called the “boundary exchange” zone which, according to Jeandel (2016) and Jeandel and Oelkers (2015), plays a major role in the land-to-ocean transfer of material (also see Fig. 2). Figure 4 illustrates the interconnection between geochemical and biological Si fluxes, particularly in the CCMZ. In agreement with Laruelle et al. (2009), Fig. 4 also shows that the open-ocean bSi production is mostly fueled by dSi inputs from below (92.5 Tmol Si yr⁻¹) and not by the CCMZ (4.7 Tmol Si yr⁻¹) (Sect. S5 in the Supplement).

5.3 The impacts of global change on the Si cycle

As illustrated by Figs. 1 and 4, the pelagic bSi production is mostly fueled from the large recycled deep-ocean pool of dSi. This lengthens the response time of the Si cycle to changes in dSi inputs to the ocean due to global change (including climatic and anthropogenic effects), increasing the possibility for the Si cycle to be out of balance.

5.3.1 Impacts on riverine inputs of dSi and aSi

Climate change at short timescales during the 21st century impacts the ocean delivery of riverine inputs of dSi and aSi (F_R) and of the terrestrial component of the submarine groundwater discharge (F_GW), either directly (e.g., dSi and aSi weathering and transport) or indirectly by affecting forestry and agricultural dSi export. So far the impacts of climate change on the terrestrial Si cycle have been reported for boreal wetlands (Struyf et al., 2010), North American (Opalinka and Cowling, 2015) and western Canadian Arctic rivers (Phillips, 2020), and the tributaries of the Laptev and East Siberian seas (Charette et al., 2020), but not for tropical environments. Tropical watersheds are the key areas for the transfer of terrestrial dSi to the ocean, as approximately 74 % of the riverine Si input is from these regions (Tréguer et al., 1995). Precipitation in tropical regions usually follows the “rich-get-richer” mechanism in a warming climate according to model predictions (Chou et al., 2004, 2008). In other words, in tropical convergence zones rainfall increases with climatological precipitation, but the opposite is true in tropical subsidence regions, creating diverging impacts for the weathering of tropical soils. If predictions of global temperature increase and variations in precipitations of the IPCC are correct (IPCC, 2018), it is uncertain how F_R or F_GW, two major components of dSi and aSi inputs, will change. Consistent with these considerations are the conclusions of Phillips (2020) on the impacts of climate change on the riverine delivery of dSi to the ocean, using a machine-learning-based approach. Phillips (2020) predicts that within the end of this century dSi mean yield could increase regionally (for instance in the Arctic region), but the global mean dSi yield is projected to decrease, using a model based on 30 environmental variables including temperature, precipitation, land cover, lithology, and terrain.

5.3.2 Abundance of marine and pelagic and benthic silicifiers

A change in diatom abundance was not seen on the North Atlantic from continuous plankton recorder (CPR) data over the period 1960–2009 (Hinder et al., 2012). However, studies have cautioned that many fields (e.g., chl) will take several decades before these changes can be measured precisely beyond natural variability (Henson et al 2010; Dutkiewicz et al 2019). The melting of Antarctic ice platforms has already been noticed to trigger impressive population blooms of highly silicified sponges (Fillinger et al., 2013).

5.3.3 Predictions for the ocean phytoplankton production and bSi production

Climate change in the 21st century will affect ocean circulation, stratification, and upwelling and therefore nutrient cycling (Aumont et al., 2003; Bopp et al., 2005, 2013). With increased stratification, dSi supply from upwelling will decrease (Figs. 1 and 4), leading to less siliceous phytoplankton production in surface compartments of lower latitudes and possibly the North Atlantic (Tréguer et al., 2018). The impact of climate change on the phytoplankton production in polar seas is highly debated as melting of sea ice decreases light limitation. In the Arctic Ocean an increase in nutrient supply from river- and shelf-derived waters (at the least for silicic acid) will occur through the Transpolar Drift, potentially impacting rates of primary production, including bSi production (e.g., Charette et al., 2020). In the SO bSi production may increase in the coastal and continental shelf zone as iron availability increases due to ice sheet melt and iceberg delivery (Duprat et al., 2016; Boyd et al., 2016; Herraiz-Borreguero et al., 2016; Hutchins and Boyd, 2016; Tréguer et al., 2018; Hawkings et al., in press). However, Henley et al. (2019) suggest a shift from diatoms to haptophytes and cryptophytes with changes in ice coverage in the western Antarctic Peninsula. How such changes in coastal environments and nutrient supplies will interplay is unknown. Globally, it is very likely that a warmer and more acidic ocean alters the pelagic bSi production rates, thus modifying the export production and outputs of Si at short timescales.

Although uncertainty is substantial, modeling studies (Bopp et al., 2005; Laufkötter et al., 2015; Dutkiewicz et al., 2019) suggest regional shifts in bSi pelagic production with climatic change. These models predict a global decrease in diatom biomass and productivity over the 21st century (Bopp et al., 2005, Laufkötter et al., 2015; Dutkiewicz et al., 2019), which would lead to a reduction in the pelagic biological flux of silica. Regional responses differ, with most models suggesting a decrease in diatom productivity in the lower latitudes and many predicting an increase in diatom productivity in the SO (Laufkötter et al., 2015). Holzer et al. (2014) suggest that changes in supply of dissolved iron (dFe) will alter bSi production mainly by inducing floristic shifts, not by relieving kinetic limitation. Increased primary productivity is predicted to come from a reduction in sea ice area, faster growth rates in warmer waters, and longer growing seasons in the high latitudes. However, many models have very simple ecosystems including only diatoms and small phytoplankton. In these models, increased primary production in the SO is mostly from diatoms. Models with more complex ecosystem representations (i.e., including additional phytoplankton groups) suggest that increased primary productivity in the future SO will be due to other phytoplankton types (e.g., pico-eukaryote) and that diatoms' biomass will decrease (Dutkiewicz et al., 2019; also see model PlankTOM5.3 in Laufkötter et al., 2015), except in regions where sea ice cover has decreased. Differences in the complexity of the ecosystem and parameterizations, in particular in terms of temperature dependences of biological process, between models lead to widely varying predictions (Laufkotter et al., 2015; Dutkiewicz et al., 2019). These uncertainties suggest we should be cautious in our predictions of what will happen with the silica biogeochemical cycle in a future ocean.

5.4 Other anthropogenic impacts

For decades if not centuries, anthropogenic activities directly or indirectly altered the Si cycle in rivers and the CCMZ (Conley et al., 1993; Ittekot et al., 2000, 2006; Derry et al., 2005; Humborg et al., 2006; Laruelle et al., 2009; Bernard et al., 2010; Liu et al., 2012; Yang et al., 2015; Wang et al., 2018a; Zhang et al., 2019). Processes involved include eutrophication and pollution (Conley et al., 1993; Liu et al., 2012), river damming (Ittekot et al., 2000; Ittekot et al., 2006; Yang et al., 2015; Wang et al., 2018), deforestation (Conley et al., 2008), changes in weathering and in river discharge (Bernard et al., 2010; Yang et al., 2015), and deposition load in river deltas (Yang et al., 2015).

Among these processes, river damming is known for having the most spectacular and short-timescale impacts on the Si delivery to the ocean. River damming favors enhanced biologically mediated absorption of dSi in the dam reservoir, thus resulting in significant decreases in dSi concentration downstream. Drastic perturbations in the Si cycle and downstream ecosystem have been shown (Ittekot et al., 2000; Humborg et al., 2006; Ittekot, 2006; Zhang et al., 2019), particularly downstream of the Nile (Mediterranean Sea), the Danube (Black Sea), and the fluvial system of the Baltic Sea. Damming is a critical issue for major rivers of the tropical zone (Amazon, Congo, Changjiang, Huanghe, Ganges, Brahmaputra, etc.), which carry 74 % of the global exorheic dSi flux (Tréguer et al., 1995; Dürr et al., 2011). Among these major rivers, the course of the Amazon and Congo is, so far, not affected by a dam or, as for the Congo River, the consequence of Congo damming for the Si cycle in the equatorial African coastal system has not been studied. The case of Changjiang (Yangtze), one of the major world players in dSi delivery to the ocean, is of particular interest. Interestingly, the Changjiang (Yangtze) River dSi concentrations decreased dramatically from the 1960s to 2000 (before the building of the Three Gorges Dam, TGD). This decrease is attributed to a combination of natural and anthropogenic impacts (Wang et al., 2018a). Paradoxically, since the construction of the TGD (2006–2009) no evidence of additional retention of dSi by the dam has been demonstrated (Wang et al., 2018a).

Over the 21st century, the influence of climate change, and other anthropogenic modifications, will have variable impacts on the regional and global biogeochemical cycling of Si. The input of dSi will likely increase in specific regions (e.g., Arctic Ocean), whilst inputs to the global ocean might decrease. Global warming will increase stratification of the surface ocean, leading to a decrease in dSi inputs from the deep sea, although this is unlikely to influence the Southern Ocean (see Sect. 5.3.3). Model-based predictions suggest a global decrease in diatom production, with a subsequent decrease in export production and Si burial rate. Clearly, new observations are needed to validate model predictions.

6 Conclusions and recommendations

The main question that still needs to be addressed is whether the contemporary marine Si cycle is at steady state, which requires the uncertainty in total inputs and outputs to be minimized.

For the input fluxes, more effort is required to quantify groundwater input fluxes, particularly using geochemical techniques to identify the recycled marine flux from other processes that generate a net input of dSi to the ocean. In light of laboratory experiments by Fabre et al. (2019) demonstrating low-temperature dissolution of quartz in clastic sand beaches, collective multinational effort should examine whether sandy beaches are major global dSi sources to the ocean. Studies addressing uncertainties at the regional scale are critically needed. Further, better constraints on hydrothermal inputs (for the northeast Pacific-specific case), aeolian input, and subsequent dissolution of minerals both in the coastal and in open-ocean zones and inputs from ice melt in polar regions are required.

For the output fluxes, it is clear that the alkaline digestion of biogenic silica (DeMaster, 1981; Mortlock and Froelich, 1989, Müller and Schneider, 1993), one of the commonly used methods for bSi determination in sediments, is not always effective at digesting all the bSi present in sediments. This is especially true for highly silicified diatom frustules, radiolarian tests, or sponge spicules (Maldonado et al., 2019; Pickering et al., 2020). Quantitative determination of bSi is particularly difficult for lithogenic or silicate-rich sediments (e.g., estuarine and coastal zones), for example those of the Chinese seas. An analytical effort for the quantitative determination of bSi from a variety of sediment sources and the organization of an international comparative analytical exercise are of high priority for future research. It is also clear that reverse weathering processes are important not only in estuarine or coastal environments, but also in distal coastal zones, slope, and open-ocean regions of the global ocean (Michalopoulos and Aller, 2014; Chong et al., 2016; Ehlert et al., 2016a; Baronas et al., 2017; Geilert et al., 2020; Pickering et al., 2020). Careful use of geochemical tools (e.g., ³²Si, $\mathrm{Ge}/\mathrm{Si}$, δ³⁰Si: Ehlert et al., 2016; Ng et al., 2020; Geilert et al., 2020; Pickering et al., 2020; Cassarino et al., in press) to trace partitioning of bSi between opal and authigenic clay phases may further elucidate the magnitude of this sink, particularly in understudied areas of the ocean.

This review highlights the significant progress that has been made in the past decade toward improving our quantitative and qualitative understanding of the sources, sinks, and internal fluxes of the marine Si cycle. Filling the knowledge gaps identified in this review is also essential if we are to anticipate changes in the Si cycle, and their ecological and biogeochemical impacts, in the future ocean.