Carbon dioxide removal via macroalgae open-ocean mariculture and sinking: an Earth system modeling study

In this study, we investigate the maximum physical and biogeochemical potential of macroalgae open-ocean mariculture and sinking (MOS) as an ocean-based carbon dioxide removal (CDR) method. Embedding a macroalgae model into an Earth system model, we simulate macroalgae mariculture in the open-ocean surface layer followed by fast sinking of the carbon-rich macroalgal biomass to the deep seafloor (depth>3000 m), which assumes no remineralization of the harvested biomass during the quick sinking. We also test the combination of MOS with artificial upwelling (AU), which fertilizes the macroalgae by pumping nutrient-rich deeper water to the surface. The simulations are done under RCP 4.5, a moderate-emissions pathway. When deployed globally between years 2020 and 2100, the carbon captured and exported by MOS is 270 PgC, which is further boosted by AU of 447 PgC. Because of feedbacks in the Earth system, the oceanic carbon inventory only increases by 171.8 PgC (283.9 PgC with AU) in the idealized simulations. More than half of this carbon remains in the ocean after cessation at year 2100 until year 3000. The major side effect of MOS on pelagic ecosystems is the reduction of phytoplankton net primary production (PNPP) due to the competition for nutrients with macroalgae and due to canopy shading. MOS shrinks the mid-layer oxygen-minimum zones (OMZs) by reducing the organic matter export to, and remineralization in, subsurface and intermediate waters, while it creates new OMZs on the seafloor by oxygen consumption from remineralization of sunken biomass. MOS also impacts the global carbon cycle by reducing the atmospheric and terrestrial carbon reservoirs when enhancing the ocean carbon reservoir. MOS also enriches dissolved inorganic carbon in the deep ocean. Effects are mostly reversible after cessation of MOS, though recovery is not complete by year 3000. In a sensitivity experiment without remineralization of sunken MOS biomass, the whole of the MOS-captured carbon is permanently stored in the ocean, but the lack of remineralized nutrients causes a long-term nutrient decline in the surface layers and thus reduces PNPP. Our results suggest that MOS has, theoretically, considerable CDR potential as an ocean-based CDR method. However, our simulations also suggest that such large-scale deployment of MOS would have substantial side effects on marine ecosystems and biogeochemistry, up to a reorganization of food webs over large parts of the ocean.

1 Introduction

Anthropogenic emissions are rapidly increasing the global atmospheric CO₂ concentration. In the last decade (2011 to 2020), global fossil CO₂ emissions averaged ∼9.49 PgC yr⁻¹ (equivalent ∼34.8 Pg CO₂ yr⁻¹) with a growth rate of 0.4 % yr⁻¹ (Friedlingstein et al., 2021). In 2019, CO₂ emissions reached a record high of 9.71±0.49 PgC yr⁻¹ (equivalent 35.6±1.8 Pg CO₂ yr⁻¹), and there is no sign of a peak (Edo et al., 2019; Friedlingstein et al., 2021). The slow speed of emission reductions until now makes it difficult to reach the promised climate goals to keep global warming within the guardrail of 2 ^∘C (Peters et al., 2013), much less the recent agreement to seriously consider an even more ambitious 1.5 ^∘C goal (UNFCCC, 2015).

In addition to mitigation efforts to reduce greenhouse gas (GHG) emissions, it is increasingly realized that carbon dioxide removal (CDR), sometimes also called negative emissions technologies (NETs), will likely be a necessary step to achieve the targets of the Paris Agreement (Minx et al., 2017; Rogelj et al., 2018). CDR aims to remove CO₂ from the atmosphere and to store it, ideally permanently, in either the terrestrial, marine or geological carbon reservoirs, thereby mitigating global warming (Glaser, 2010). Due to the limited remaining emission budget (650±130 Pg CO₂ to 1.5 ^∘C and 1300±130 Pg CO₂ to 2 ^∘C), deployment of CDR is required in most pathways studied in the scientific literature to achieve these ambitious targets (Lawrence et al., 2018; IPCC, 2018). As the second-largest inorganic carbon reservoir on the planet, the ocean plays a pivotal role in naturally regulating the atmospheric CO₂ concentration. Since the beginning of the industrial era, the ocean has taken up more than 560 Pg CO₂, about 25 % of the anthropogenic CO₂ emissions (∼2030 Pg CO₂; Gruber et al., 2019; Ciais et al., 2013; Heinze et al., 2015). Its high carbon storage capacity could theoretically match or exceed fossil fuel resources (Scott et al., 2015). Thus, a variety of ocean-based CDR methods have been proposed to take advantage of this potential storage capacity. The proposed ocean-based CDR approaches aim to increase the rate of oceanic CO₂ uptake and storage by either enhancing abiotic processes (i.e., chemical or physical, e.g., ocean alkalinization; Keller et al., 2014; Taylor et al., 2016; Köhler et al., 2013; Albright et al., 2016) or biotic processes (e.g., ocean fertilization; Keller et al., 2014; Smetacek et al., 2012; Oschlies et al., 2010b; Matear and Elliott, 2004; Robinson et al., 2014). Some technologies also seek to remove CO₂ directly from seawater and to store it in some other reservoir, e.g., a geological one (Eisaman et al., 2012).

Macroalgae species (also known as “seaweed” or “kelp”) are highly efficient carbon fixers with a high C:N ratio (Atkinson and Smith, 1983; Fernand et al., 2017) and observed net primary production (NPP) rates of 91–522 gC m² yr⁻¹. In the 1970's, the concept of ocean farming using macroalgae for marine carbon sink and bioenergy production was studied with an actual small test farm established off the coast of southern California. These research activities were abandoned due to the damage of the test farm by winter storms and for several technical and economic reasons (Ritschard, 1992). Utilizing macroalgae for biological ocean-based CDR has recently received renewed interest (Duarte et al., 2017; Chung et al., 2011; Gao et al., 2020; Fernand et al., 2017; Raven, 2018). The macroalgae aquaculture industry is well established globally, with an annual harvest of over 30 million t wet weight (WW, FAO, 2018). Thus, some proposals have focused on using harvested macroalgae for producing biochar (Roberts et al., 2015; Bird et al., 2011) or bio-energy combined with carbon capture and storage (BECCS, Chung et al., 2011; Buschmann et al., 2017; Gao and McKinley, 1994; Chen et al., 2015; Fernand et al., 2017). However, as current macroalgae aquaculture facilities are mainly located in coastal regions, the scope to expand macroalgae aquaculture is limited by the shortage of suitable coastal areas due to nutrient availability and shifting temperature regimes (Duarte et al., 2017; Oyinlola et al., 2020). To address these issues, several offshore macroalgae aquaculture facilities have been designed and evaluated (e.g., the SeaweedPaddock by Sherman et al. (2019), the offshore ring by Buck and Buchholz (2004) and the depth-cycling strategy by Navarrete et al. (2021), in which macroalgae are physically towed into the deep nutrient-rich water at night). Moreover, the Advanced Research Projects Agency–Energy (ARPA-E) of the US Department of Energy (DOE) has committed more than 60 million dollars on the Macroalgae Research Inspiring Novel Energy Resources (MARINER) program to develop the technologies for macroalgal biomass production, including integrated ocean cultivation and harvesting systems (APAR-e, 2021). Thus, the ideas of expanding macroalgae cultivation to the open oceans (mariculture) are ambitious but no longer fictional, and they provide a theoretical possibility to expand macroalgae aquaculture to the open ocean for CDR.

In this study, we evaluate “macroalgae open-ocean mariculture and sinking (MOS)” as an ocean-based CDR method that is designed to artificially enhance the macroalgae-based carbon dioxide removal. The aim of this study is to investigate (1) the maximum physical and biogeochemical CDR potential of MOS, (2) the side effects of such large-scale deployment, and (3) to understand where offshore macroalgae farming would be viable if done at a large scale. This information is needed to help prioritize further research into CDR, to understand if there are potential MOS side effects that become evident only at a large scale and to provide information on the viability of large-scale offshore macroalgae farming in different regions over time by accounting for the implications of nutrient utilization and climate change.

To do this, simulated macroalgae are seeded and cultivated on offshore floating platforms that are moored to the seabed (e.g., see platform designs in Buck and Buchholz, 2004). The platforms are also assumed to float below the open-ocean surface (at 5 m depth) to avoid storm damages. At the end of an annual cycle, platforms with matured macroalgae are rapidly sunk to the seafloor and unload the biomass there. This can be thought of as a short circuiting of the biological pump by bringing marine biomass directly to the seafloor without having it remineralized along the way. Afterwards, the sunken biomass is assumed to continue remineralization at the seafloor, consuming oxygen and releasing dissolved inorganic carbon (C), nitrogen (N) and phosphate (P) into the deep ocean where it ideally remains for centuries to millennia (Fig. 1). The macroalgae used here is an idealized genus. The assumed constant $\mathrm{C}:\mathrm{N}:\mathrm{P}$ ratio is $\mathrm{400}:\mathrm{20}:\mathrm{1}$, which is higher than the stoichiometric ratio of the general phytoplankton in the University of Victoria Earth system climate model (UVic ESCM; $\mathrm{C}:\mathrm{N}:\mathrm{P}$ $=\mathrm{106}:\mathrm{16}:\mathrm{1}$, the Redfield ratio). In practice, some of the biomass may also be permanently buried in sediments (Luo et al., 2019; Sichert et al., 2020), and we will explore the extreme case of zero remineralization in the water column in a sensitivity experiment. In another sensitivity experiment, we investigate combining MOS with artificial upwelling (AU) to alleviate nutrient limitation in the open-ocean surface (Duarte et al., 2017; Kim et al., 2019; Laurens et al., 2020).

Figure 1 Schematic illustrating the biogeochemical fluxes and physical impacts of MOS on nutrients (NO₃ and PO₄), oxygen, dissolved inorganic carbon (DIC), ordinary phytoplankton (P_O in green), diazotrophs (P_D in pale brown) and zooplankton (Z).

To investigate the biogeochemical and climatic implications of MOS, we use an Earth system model of intermediate complexity. Though the idea of massive macroalgae cultivation and biomass offsetting for CDR has been assessed in some earlier publications (Orr and Sarmiento, 1992; Gao and McKinley, 1994; Froehlich et al., 2019; Lehahn et al., 2016), as far as we are aware it has not been evaluated using an Earth system model (ESM). ESM-based assessments are required for studying the response of the global carbon cycle to such perturbations and for estimating the efficacy of such methods in a global carbon cycle context (with regards to atmospheric CO₂ removal). Furthermore, such models can dynamically simulate macroalgae growth and the permanence of carbon storage (i.e., the fate of sunken biomass on the seafloor), as well as their interactions with global marine biogeochemistry. It is essential to clarify these issues before any decisions about eventual implementation of MOS can be made.

2 Methods

2.1 Model description

In this study, we employ the University of Victoria Earth system climate model (UVic ESCM) version 2.9 (Weaver et al., 2001; Eby et al., 2009; Keller et al., 2012), which consists of three dynamically coupled components: a three-dimensional ocean circulation model (Pacanowski, 1996) including a dynamic–thermodynamic sea-ice model (Bitz and Lipscomb, 1999), a terrestrial model (Meissner et al., 2003; Weaver et al., 2001) and a simple one-layer atmospheric energy–moisture balance model (Fanning and Weaver, 1996). The model has a fully coupled carbon cycle including dynamic terrestrial, atmospheric and oceanic carbon inventories. The horizontal resolution of all components is 3.6^∘ longitude $\times \mathrm{1.8}{}^{\circ }$ latitude, and the ocean component has 19 vertical layers. The descriptions of air–sea gas exchange and seawater carbonate chemistry are based on the Ocean Carbon Cycle Model Intercomparison Project (OCMIP) abiotic protocol (Orr et al., 1999). The ocean biogeochemistry is presented with a nutrients–phytoplankton–zooplankton–detritus (NPZD) model that includes one general phytoplankton, diazotrophs, and one zooplankton type (Keller et al., 2012; Eby et al., 2013). The UVic ESCM has been evaluated in several recent studies (e.g., Keller et al., 2014; Mengis et al., 2016; Reith et al., 2016; Kvale et al., 2021).

2.2 Modeling MOS in the UVic ESCM

In this study, the modeling of macroalgae is done with a macroalgae growth model coupled into the UVic ESCM. The aim of the macroalgae model is to investigate the carbon sequestration capacity of MOS and the potential impacts on marine biogeochemistry. In the macroalgae model, the net growth rate is affected by several limiting factors, including nutrients, temperature and solar radiation intensity. The cellular $\mathrm{C}:\mathrm{N}:\mathrm{P}$ ratio of macroalgae is fixed. The loss of macroalgal biomass includes erosion and grazing by zooplankton. The deployment of MOS is done with an algorithm that considers spatial and temporal conditions.

The macroalgae model is also connected to global marine biogeochemical processes, including the inorganic carbon and nutrient pools. In the surface layers, it impacts phytoplankton via nutrients competition and canopy shading. The single, aggregated zooplankton compartment of the biogeochemical model, which represents higher trophic levels, is also designed to graze on macroalgae. In the bottom layers, the remineralization of sunken macroalgal biomass will consume the dissolved oxygen, which in turn limits the rate of remineralization.

2.2.1 Macroalgae model

The macroalgae model is an idealized generic model of genus Laminaria and Saccharina, mainly based on Martins and Marques (2002) and Zhang et al. (2016). The rate of biomass change is governed by Eq. (1) as the imbalance of NGR (net growth rate, d⁻¹) and LR (loss rate, fraction of daily biomass loss due to mortality, erosion and grazing by zooplankton, d⁻¹).

Modeled macroalgae is seeded 5 m below sea surface, considering the light requirement and reduction of damaging risks by surface wave activity (Eq. 11). The deployment of macroalgae considers ambient nutrients availability and avoidance of winter periods (Sect. 3.1).

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{dBiomass}}{\mathrm{d}t}}=(\mathrm{NGR}-\mathrm{LR})\times \mathit{Biomass}\end{array}$$

The NGR is regulated by

$$\begin{array}{}\text{(2)}& \mathrm{NGR}=R_{\mathrm{growth}}-R_{\mathrm{resp}},\end{array}$$

where R_growth is the gross growth rate (d⁻¹), and R_resp is the respiration rate (d⁻¹). The growth rate of macroalgae (R_growth) is given by Eq. (3), regulated by water temperature (T), solar irradiance (I) and dissolved nutrient concentrations (NO₃ and PO₄, NP).

$$\begin{array}{}\text{(3)}& R_{\mathrm{growth}}={\mathit{\mu }}_{\max }\times f\left(T_{\mathrm{w}}\right)\times f\left(\mathrm{NP}\right)\times f\left(I_{\mathrm{ma}}\right)\end{array}$$

In the current model, the macroalgal growth rates are controlled by external concentrations of available nutrients via assumed Michaelis–Menten kinetics with half-saturation constants K_N and K_P for NO₃ and PO₄, respectively:

$$\begin{array}{}\text{(4)}& {\displaystyle }f\left(\mathrm{N}\right)={\displaystyle \frac{{\mathrm{NO}}_{\mathrm{3}}}{K_{\mathrm{N}}+{\mathrm{NO}}_{\mathrm{3}}}},\text{(5)}& {\displaystyle }f\left(\mathrm{P}\right)={\displaystyle \frac{{\mathrm{PO}}_{\mathrm{4}}}{K_{\mathrm{P}}+{\mathrm{PO}}_{\mathrm{4}}}},\text{(6)}& {\displaystyle }f\left(\mathrm{NP}\right)=\mathrm{Min}\mathit{\left\{}f\right(\mathrm{N}),f(\mathrm{P}\left)\mathit{\right\}}.\end{array}$$

The need for iron is not considered in our macroalgae growth model. Although iron is utilized during macroalgae growth (e.g., Suzuki et al., 1995; Kuffner and Paul, 2001), iron limitation on macroalgae is not widely discussed, especially for the genus Saccharina. Besides, as iron is a micronutrient needed in low quantities, the MOS platform could be designed with an iron supply for the macroalgae, in which case MOS could be considered to include a targeted variant of the ocean iron fertilization concept.

The temperature limiting factor used here is an optimum curve following Bowie et al. (1985). T_opt is the species-specific optimum temperature at which the growth rate is maximized. T_max and T_min define the upper and lower temperature limits above and below which macroalgae growth ceases. The temperature optimum curve of the macroalgae is shown in Fig. A3.

$$\begin{array}{}\text{(7)}& {\displaystyle }f\left(T_{\mathrm{w}}\right)=e^{-\mathrm{2.3}\times X_T^{\mathrm{2}}}\text{(8)}& {\displaystyle }{X}_T={\displaystyle \frac{T_{\mathrm{w}}-T_{\mathrm{opt}}}{T_x-T_{\mathrm{opt}}}}\text{(9)}& {\displaystyle }{T}_x=\left\{\begin{array}{lr}{T}_{\min }& \text{ if }{T}_{\mathrm{w}}\le T_{\mathrm{opt}}\\ T_{\max }& \text{ if }{T}_{\mathrm{w}}>T_{\mathrm{opt}}\end{array}\right.\end{array}$$

Respiration is described by an Arrhenius function considering water temperature T_w in degrees Celsius (Duarte and Ferreira, 1997; Martins and Marques, 2002):

$$\begin{array}{}\text{(10)}& R_{\mathrm{resp}}=R_{\max \mathrm{20}}\times r^{(T_{\mathrm{w}}-\mathrm{20})},\end{array}$$

where R_max20 is the maximum respiration rate of the simulated macroalgae species at 20 ^∘C, and r stands for the empirical coefficient for macroalgae respiration (Table 1).

The limiting factor of solar irradiance density for macroalgae photosynthesis (f(I_ma)) is given in Eq. (11) (Steele's photo-inhibition relationship – Kirk, 1994):

$$\begin{array}{}\text{(11)}& f\left(I_{\mathrm{ma}}\right)={\displaystyle \frac{I_{\mathrm{ma}}}{I_{\mathrm{opt}}}}\times e^{(\mathrm{1}-\frac{I_{\mathrm{ma}}}{I_{\mathrm{opt}}})},\end{array}$$

where I_ma stands for the shortwave radiation intensity reaching the depth Z (given by Eq. 12), and I_opt stands for the optimum light intensity for macroalgae growth (constant, Table 1).

$$\begin{array}{}\text{(12)}& I_{\mathrm{ma}}=I_{\mathrm{s}}\times e^{\left(-k_{\mathrm{w}}{Z}_{\mathrm{m}}-{\int }_{\mathrm{0}}^{Z_{\mathrm{m}}}\left(P_o+P_{\mathrm{D}}\right)k_{\mathrm{c}}\times \mathrm{d}{Z}_{\mathrm{m}}\right)}\times \mathrm{d}f\end{array}$$

Equation (12) calculates the shortwave radiation (I_ma) reaching the depth Z_m. This is modified from Keller et al. (2012) and Schmittner et al. (2005), with Z_m denoting the depth of MOS macroalgae platforms beneath the water surface. Z_m is assumed to be 5 m, compromising the empirical depth with sufficient light for macroalgae photosynthesis (1 to 2 m for cultivation (Buck and Buchholz, 2004), 0 to 10 m for wild communities (Eriksson and Bergström, 2005)) and the depth to reduce the risks of damage by stressful turbulence or severe weather events (e.g., hurricanes). df denotes the day length as a fraction of 24 h. I_s stands for the shortwave radiation density at the top of the layer. P_O and P_D are the biomass of ordinary phytoplankton and diazotrophs, respectively, in the layers above the macroalgae. k_w is the light attenuation coefficient for water. I_opt is the optimum light intensity for macroalgae growth. k_c is the light attenuation coefficient of phytoplankton and also accounts for co-varying particulate and dissolved inorganic and organic materials (Kvale and Meissner, 2017). As described in Sect. 2.2.1, the morphology of the frond will not be considered, and the self-shading effects by fronds are not considered here (Duarte and Ferreira, 1997; Brush and Nixon, 2010).

The loss rate LR is regulated by the following:

$$\begin{array}{}\text{(13)}& {\displaystyle }\mathrm{LR}=\mathrm{ER}+{\mathrm{Graze}}_{\mathrm{ma}},\text{(14)}& {\displaystyle }\mathrm{ER}=\mathrm{Biomass}\times R_{\mathrm{erosion}},\text{(15)}& {\displaystyle }{\mathrm{Graze}}_{\mathrm{ma}}={\mathit{\mu }}_{\mathrm{Z}}^{\max }Z\times {\mathit{\psi }}_{\mathrm{ma}}\times \mathrm{Biomass},\end{array}$$

where the erosion of biomass (ER) is controlled by the individual erosion rate R_erosion. As the frond morphology of macroalgae is not modeled here, we set the R_erosion as a constant independent of physical impacts (Trancoso et al., 2005; Zhang et al., 2016). The eroded macroalgal biomass will be directly converted back to nutrients and DIC (dissolved inorganic carbon) according to the macroalgae stoichiometry ratios without remineralization or further degradation by zooplankton. We set this parameterization of erosion, a small biomass loss of 0.01 % per day, pragmatically as instantaneous remineralization rather than introducing another finite remineralization and finite sinking parameterization with difficult-to-constrain parameters. It was used to minimize the computational expense of the model and to avoid having to add another state variable that is subjected to physical transport.

Graze_ma is the biomass loss due to grazing by zooplankton. Z is the zooplankton biomass which is calculated by the NPZD model. ${\mathit{\mu }}_{\mathrm{Z}}^{\max }$ stands for the maximum potential growth rate of zooplankton, defined in Keller et al. (2012, Eq. 28). The zooplankton grazing preference on macroalgae (ψ_ma) is defined in Sect. 2.2.3.

For simplicity and to limit the number of state variables, we made the following modifications to the macroalgae model:

1. We did not include a dynamic $\mathrm{C}:\mathrm{N}:\mathrm{P}$ ratio or a representation of luxury nutrient uptake and storage (Broch and Slagstad, 2012; Hadley et al., 2015). Instead, the $\mathrm{C}:\mathrm{N}:\mathrm{P}$ ratio of the macroalgae biomass was set as a constant (Table 1), which is based on seasonally averaged measurements of the biomass composition of these genera (Zhang et al., 2016; Martins and Marques, 2002).

2. The macroalgae life cycle processes (e.g., alternations of generations) are also not considered in our model (Brush and Nixon, 2010; Trancoso et al., 2005; Duarte and Ferreira, 1997). We thus assumed that the plantlet (e.g., sporophytes for Saccharina) will be reseeded annually on the MOS infrastructure. The assumed deployment strategy, i.e., timing of seeding and sinking, of MOS is latitude dependent according to the seasonality of solar irradiance (see Sect. 3.1). Whenever conditions are unfavorable for macroalgae and no growth occurs during an annual cycle, no re-seeding of macroalgae will occur in these regions.

The parameterization of DOC release by macroalgae could not be included because of the lack of enough information. Few studies exist, and the uncertainties about the release of DOC are large. For example, Barrón et al. (2014) reported a release of DOC by macroalgae from a few species of 23.2±12.6 mmol C m⁻² d⁻¹ with no information on bioavailability. Meanwhile, refractory DOC dynamics are difficult to include in a global Earth system model and are beyond the scope of this study (Anderson et al., 2015; Mentges et al., 2019; Zakem et al., 2021). Thus, the DOC release from macroalgae is not included in this study.

2.2.2 Remineralization of sunken macroalgal biomass

Biomass sinking is simulated by instantly transferring the macroalgal biomass from the surface grid cell to the deepest grid cell at the respective location at the end of each cultivating period. This assumes that the harvested biomass could be engineered to sink to the seafloor in a rapid and efficient manner with no remineralization along the way. Afterwards, the next macroalgae generation will start to grow in the surface layer. Equation (16) calculates the temperature-dependent remineralization rate of sunken macroalgal biomass (μ_ma) following the function of remineralization of detritus in the UVic ESCM (described in Schmittner et al., 2008, Eq. A16). Remineralization consumes oxygen and returns DIC, PO₄ and NO₃ from the sunken macroalgal biomass to the sea water and is described as

$$\begin{array}{}\text{(16)}& {\mathit{\mu }}_{\mathrm{ma}}={\mathit{\mu }}_{{\mathrm{ma}}_{\mathrm{0}}}\exp \left(T_{\mathrm{w}}/T_{\mathrm{b}}\right)\left[\mathrm{0.65}+\mathrm{0.35}\tanh \left({\mathrm{O}}_{\mathrm{2}}-\mathrm{6}\right)\right],\end{array}$$

where μ_ma0 is the remineralization rate of sunken macroalgal biomass at 0 ^∘C. T_w and T_b represent the sea water temperature and e-folding temperature of biological rates, and O₂ is the dissolved-oxygen concentration in mmol m⁻³. When the dissolved oxygen is insufficient (<5 mmol m⁻³), aerobic remineralization will be replaced by oxygen-equivalent, but slower, denitrification via reduction of NO₃ (Keller et al., 2012). Note that remineralization will cease when NO₃ is also completely consumed.

There are considerable uncertainties concerning the fate of sunken macroalgae (Sichert et al., 2020; Krause-Jensen and Duarte, 2016; Luo et al., 2019). A sensitivity simulation explores the situation where ${\mathit{\mu }}_{{\mathrm{ma}}_{\mathrm{0}}}$ is set to zero, which would assume permanent deposition of the sunken biomass on (or in) the seafloor without decaying (Sect. 3.2).

2.2.3 Interactions with pelagic microbial ecosystems

Besides the competition for nutrient resources, the macroalgae canopies may also reduce solar irradiance downward (“canopy shading”) and thus limit phytoplankton photosynthesis beneath the macroalgae (Jiang et al., 2020). Equation (17) (modified from Eq. 14, Keller et al., 2012) describes the shortwave radiation attenuation (I_phyt) through the macroalgae layer, as well as through phytoplankton and water (MOS is not deployed in areas covered by sea ice):

$$\begin{array}{}\text{(17)}& I_{\mathrm{phyt}}=I_{\mathrm{s}}\times e^{-k_{\mathrm{w}}Z-{\int }_{\mathrm{0}}^{\mathrm{Z}}\left(P_{\mathrm{O}}+P_{\mathrm{D}}\right)k_{\mathrm{c}}\times \mathrm{d}Z-k_{\mathrm{ma}}\times h_{\mathrm{ma}}\times \mathrm{Biomass}},\end{array}$$

where k_ma, the macroalgae light extinction coefficient (m⁻¹), is calculated based on the biomass of macroalgae in carbon as follows:

$$\begin{array}{}\text{(18)}& k_{\mathrm{ma}}=a_{\mathrm{ma}}\times \mathrm{Biomass}\times {\mathrm{MR}}_{\mathrm{C}:\mathrm{N}}.\end{array}$$

Here, a_ma is the macroalgae carbon-specific shading area (m² kgC⁻¹; Trancoso et al., 2005), h_ma is the thickness of macroalgae layer, and MR_C:N stands for the molar C:N ratio of macroalgal biomass.

The original NPZD model in Keller et al. (2012) is extended by allowing zooplankton to graze on macroalgae. Our assumption that zooplankton can graze on macroalgae is based on the notion that the marine biogeochemical component of “zooplankton” in UVic ESCM represents all higher trophic levels, including known macroalgae grazers such as amphipods (Jacobucci et al., 2008), gastropods (Chikaraishi et al., 2007; Krumhansl and Scheibling, 2011), sea urchins (e.g., Yatsuya et al., 2020) and fishes (e.g., Peteiro et al., 2014). Thus, we included this food web pathway to assess the sensitivity of macroalgae to potential grazers in the ocean, assuming that with large macroalgae farms the pelagic larvae of some grazing organisms like fish or urchins would settle within the farms.

The grazing preference for macroalgae (ψ_ma) is set to 1$\times {\mathrm{10}}^{-\mathrm{4}}$ according to observational studies (Trancoso et al., 2005). Macroalgae thus provide a grazing option for zooplankton in addition to the traditional NPZD-type model food sources (phytoplankton, diazotrophs, detritus and zooplankton via self-grazing). Therefore, the four original grazing preferences (0.3 on phytoplankton, 0.1 on diazotrophs, 0.3 on detritus and 0.3 on zooplankton; Keller et al., 2012, Table 1) are reduced by $\frac{\mathrm{1}}{\mathrm{4}}$ ψ_ma each. In the areas where MOS is absent (i.e, in the ice-covered ocean surface), the zooplankton grazing will follow the original description in Keller et al. (2012, Table 1) without the preference for macroalgae. No CaCO₃ formation by macroalgae is simulated here (Bach et al., 2021; Macreadie et al., 2017, 2019), as calcareous macroalgae species and epibiont calcifiers are not considered. Therefore, the only alkalinity impact of growing and remineralizing macroalgae comes via changes in nitrate and phosphate.

2.2.4 Mass conversions

In order to parameterize and validate the model, it is necessary to convert from commonly measured macroalgae variables (often in wet- and dry-weight units) to the model units. These conversions include the calculation of carbon and CO₂ sequestered in macroalgal biomass (C_ma, gram carbon), as well as the conversions of dry weight (DW, gram) and wet weight (WW, gram):

$$\begin{array}{}\text{(19)}& {\displaystyle }{C}_{\mathrm{ma}}=\mathrm{Biomass}\times {\mathrm{MR}}_{\mathrm{C}:\mathrm{N}}\times \mathrm{12.011},\text{(20)}& {\displaystyle }{\mathrm{CO}}_{\mathrm{2}\mathrm{ma}}={\mathrm{C}}_{\mathrm{ma}}\times \mathrm{3.67},\text{(21)}& {\displaystyle }\mathrm{DW}={\mathrm{C}}_{\mathrm{ma}}÷{\mathrm{MR}}_{\mathrm{C}:\mathrm{DW}},\text{(22)}& {\displaystyle }\mathrm{WW}=\mathrm{DW}\times {\mathrm{MR}}_{\mathrm{DW}:\mathrm{WW}},\end{array}$$

where 3.67 is the ratio between the atomic mass of CO₂ (44 g mol⁻¹) and carbon (12 g mol⁻¹), Biomass is in moles of nitrogen, and 12.011 is the relative molecular weight of carbon (g mol⁻¹).

2.2.5 MOS carbon retained in the ocean and outgassing

The DIC from remineralization of sunken biomass will eventually be conveyed back to the ocean surface and may leak back to the atmosphere. Equation (23) calculates the ocean-retained fraction (FR, %) of MOS-captured carbon (MOS-C), where the C_captured is carbon in the cumulative sunken biomass, C_{Sunk Biomass} is the carbon in sunken macroalgal biomass that still remains on the seafloor.

$$\begin{array}{}\text{(23)}& \mathrm{FR}={\displaystyle \frac{{\mathrm{C}}_{\mathrm{retained}}}{{\mathrm{C}}_{\mathrm{captured}}}}={\displaystyle \frac{\left({\mathrm{DIC}}_{\mathrm{remineralized}}+{\mathrm{C}}_{\mathrm{Sunk}\mathrm{Biomass}}\right)}{{\mathrm{C}}_{\mathrm{captured}}}}\end{array}$$

In order to track the leakage of MOS-C after remineralization, a tracer of remineralized MOS-C (MOS_DIC) is added to the UVic ESCM in addition to the original DIC tracer. MOS_DIC participates in the inorganic ocean carbon cycle (Weaver et al., 2001, Sect. 3e). When reaching the surface, the outgassing of MOS_DIC will follow the air–sea gas exchange process in UVic ESCM, which is given in Weaver et al. (2001, Sect. 3e). The air–sea exchange flux of MOS-C is also calculated for analyzing the location and quantity of outgassing. The results of MOS-C outgassing are shown in Sect. 4.6.

Trancoso et al. (2005)Van Der Molen et al. (2017)Zhang et al. (2016)Zhang et al. (2016)Keller et al. (2012)Keller et al. (2012)Zhang et al. (2016)Zhang et al. (2016)Keller et al. (2012)Keller et al. (2012)Trevathan-Tackett et al. (2015)Atkinson and Smith (1983)Atkinson and Smith (1983)Aldridge and Trimmer (2009); Conover et al. (2016); Van Der Molen et al. (2017)Chung et al. (2011)Duarte et al. (2003)Martins and Marques (2002); Zhang et al. (2016)Martins and Marques (2002)Van Der Molen et al. (2017)Schmittner et al. (2008)Zhang et al. (2016); Martins and Marques (2002)Breeman (1988)Martins and Marques (2002)Zhang et al. (2016)Oschlies et al. (2010b)Trancoso et al. (2005)Partanen et al. (2016)

Table 1 Model parameters.

3 Experiment design

The UVic ESCM is spun up for >10 000 years to an equilibrium state under pre-industrial (year 1765) atmospheric and astronomical boundary conditions and is then integrated for another 250 years without prescribing atmospheric CO₂ concentrations to allow the carbon cycle to equilibrate. Afterwards, the model is run from 1765 until 2005 and forced with historical fossil fuel emissions and land-use changes (crop and pastureland; Keller et al., 2014). From year 2005 to 2100, simulations are forced with CO₂ emissions represented as a direct adjustment to radiative forcing, land-use change by agriculture, volcanic radiative forcing and sulfate aerosols which are prescribed according to the Representative Concentration Pathway 4.5 (RCP 4.5; Meinshausen et al., 2011; Thomas, 2014; Keller et al., 2014; Partanen et al., 2016). Solar insolation at the top of the atmosphere, wind stress and wind fields are varied seasonally. After the year 2300, CO₂ emissions are assumed to decrease linearly until the end of year 3000 with other forcing held constant to zero.

The full list of simulations is given in Table 3. To test the maximum potential, as well as the global carbon cycle and biogeochemical responses, we simulate MOS for 1000 years, beginning in year 2020 (MOS_Conti). Additionally, termination experiments (MOS_Stop) are performed to analyze the response of the ocean and climate to an abrupt termination of MOS at year 2100.

3.1 Deployment strategies of MOS

The current study focuses on estimating the maximum carbon sequestration potential of MOS and assumes instantaneous seeding on floating infrastructure in the open ocean. The macroalgae is represented as a biogeochemical tracer (Eq. 1) that is not subject to physical transports and remains fixed in the top (first) ocean layer of the UVic ESCM, which is assumed to be well mixed. In our idealized experiments, MOS deployment must fulfill the following requirements:

• The water depth must be ≥3000 m; according to the assessment by Reith et al. (2016), leakage of dissolved inorganic carbon added to deep waters (in this case from remineralization of sunken macroalgal biomass) is small at such depths compared to shallower ones.

• The ambient surface NO₃ concentration is greater than Seed plus K_N (Table 1); this ensures sufficient nutrients for initial growth, as Seed is directly transferred from dissolved NO₃, and K_N is the half-saturation constant for NO₃ uptake. Note that in this calculation, Seed has been converted from the unit of kgC km⁻¹ in Table 1, the unit of concentration of nitrogen (mmol N m⁻³).

• It must be spatially located between 57^∘ N and 72^∘ S to remain in sea-ice-free waters.

Note that the DIC, N and P components of the initial Seed are directly removed from the inorganic-matter pool of the respective grid box in order to maintain model mass balance and to avoid adding extra nutrients and/or carbon to the ocean at the time of seeding.

During the MOS simulations, the seasonality of temperature and of solar radiation is an essential limiting factor of the primary productivity of MOS in various latitudinal regions. In order to avoid the unnecessary loss of macroalgal biomass during winter periods when solar radiation is insufficient and when the ambient water temperature is low, we partitioned the global ocean surface into three belts (N, M and S) and pragmatically applied farming strategies according to Table 2. The period between the seeding and sinking of macroalgae is set as 6 months, from May to October in belt N and from November to the next April in belt S. In belt M, the macroalgae is seeded at the beginning of the year, and sinking occurs after 12 months. The geographical locations of the three belts are shown in Fig. 2.

Table 2 Latitudinal division of MOS deployment regions.

^* In the following year.

The maximum yield of Biomass is set to a constant value of Y_max (Table 1). When the biomass reaches Y_max in a grid cell, macroalgae will stop growing and wait for sinking. After an annual cultivation cycle, the macroalgal biomass is instantaneously delivered to the seafloor, except for a small fraction (equivalent to Seed) that remains at the surface for re-seeding. In some regions where conditions are unfavorable and no net macroalgae growth had occurred during the last cultivation period, the total Biomass will be sunk once without any further re-seeding. In order to prevent MOS from removing too much atmospheric CO₂ in long-term simulations where emissions eventually reach zero, MOS deployment will be terminated once atmospheric CO₂ concentration hits 280 ppm, assuming that there is no need for more CDR once pre-industrial CO₂ values have been reached.

3.2 Sensitivity studies

As test simulations indicated that the CDR potential of MOS is, in many ocean regions, limited by the availability of nutrients in the surface layer, sensitivity simulations were performed with MOS combined with artificial upwelling (AU) that pumps up nutrient‐rich deeper waters to the surface and thereby relaxes nutrient stress and enhances the macroalgae growth.

The simulated MOS-AU system is based on Oschlies et al. (2010b) and Keller et al. (2014). We placed modeled “pipes” that pump deeper water to the ocean surface in areas where MOS is deployed. The simulated upwelling works by transferring water adiabatically from the grid box at the lower end of the pipe to the surface grid box at a rate of 1 cm d⁻¹. These pipes will function continuously until the termination of MOS (in year 2100 or 3000). However, because these earlier studies have revealed a dominant effect associated with the low temperatures of the upwelled colder waters, we here concentrate on the nutrient aspect and simulate a hypothetical MOS-AU system that keeps temperatures at ambient levels (e.g., via heat exchangers).

In the simulated AU system, water, together with dissolved tracers, is transferred from the grid box at the lower end of the pipes to the surface grid box, resulting in a model grid box-average upwelling rate (w, set to 1 cm d⁻¹; Table 1). The lower end of the pipes is fixed at a depth of 1000 m. Similarly to the normal MOS simulations without AU, the MOS-AU simulations are deployed from year 2020 and then terminated at either year 2100 in a discontinuous run or at year 3000 in a continuous one (Table 3).

The MOS-AU joint system is deployed using the following strategies: AU pipes will be deployed everywhere with depth ≥3000 m and start upwelling immediately. If surface nutrient concentrations are raised to the initial seeding condition (Sect. 3.1) in any grid box, MOS will be deployed, thereby expanding the range where MOS can grow.

Another model parameter selected for sensitivity studies is the remineralization rate of sunken macroalgal biomass (μ_ma, Eq. 16). μ_ma is a critical factor impacting the residence time of MOS-captured carbon in the ocean and associated benthic oxygen consumption by remineralization. Macroalgal biomass has been reported to be recalcitrant to microbial degradation; however, the fate of macroalgal biomass in the deep sea is uncertain (Krause-Jensen and Duarte, 2016; Luo et al., 2019; Sichert et al., 2020). Thus, an extreme and idealized situation with μ_ma set to zero is tested in sensitivity simulations (MOS-NoRe, MOS with zero mineralization of sunk macroalgal biomass). This can be thought of as a case where all biomass is permanently buried upon reaching the seafloor. This sensitivity study can also simulate an extreme case of infinitely slow remineralization, which can help in estimating the range of possible fates of remineralized organic matter. Meanwhile, this sensitivity study also represents a different macroalgae farming approach – that of harvesting the biomass to create bioenergy with carbon capture or storage (BECCS) or biochar (e.g., Kerrison et al., 2015; Laurens et al., 2020; Roberts et al., 2015), with the assumption that all harvested biomass was permanently removed from the ocean. While this is a very idealized case, it serves the useful purpose of providing information on how marine biogeochemistry is impacted by the permanent removal of fixed C, N and P.

The stoichiometric C:N ratio of macroalgal biomass (MR_C:N) may also influence the CDR capacity of MOS. In the current study, the MR_C:N (400:20, Table 1) is set as 20, nearly 2 times higher than the phytoplankton stoichiometric biomass C:N ratio in the UVic ESCM ($\mathrm{C}:\mathrm{N}=\mathrm{106}:\mathrm{16}$, the Redfield ratio). However, the difference between the macroalgae and phytoplankton stoichiometric C:N ratio may have strong influences on the CDR potential of MOS. For instance, Bach et al. (2021) have indicated that the CDR potential of a floating-macroalgae (Sargassum) belt may be reduced by 7 % to 50 % due to the nutrient reallocation caused by the variation of gaps of C:N between macroalgae and phytoplankton. Thus, sensitivity experiments of macroalgal C:N ratios (MOS_Conti_CN_High and MOS_Conti_CN_Low) have been performed to investigate the impacts of MR_C:N on the MOS CDR capacity.

Table 3 Description of the model experiments. “Stop” represents the termination of the simulation in year 2100; “Conti” represents the continuous MOS deployment till year 3000; “AU” represents artificial upwelling; NoRe represents zero remineralization of sunken macroalgal biomass; CN represents the molar C:N ratio of macroalgal biomass (MR_C:N in Table 1).

4 Results

4.1 Evaluation of MOS

To evaluate if the simulated MOS systems have plausible macroalgae growth characteristics, we evaluate the seasonal dynamics of the simulated MOS system for a 30 d averaged time slice from 2020 to 2024 under the RCP 4.5 emission scenario and without artificial upwelling.

4.1.1 Distribution of MOS

The red contours in Fig. 2a delineate the occupied area that basically follows the pattern of the simulated NO₃-rich ocean surface (Keller et al., 2012, Fig. 9; Garcia et al., 2010, WOA2009 Dataset) in the northern and the equatorial Eastern Pacific, as well as in the Southern Ocean. Except for the coastal regions and Arctic areas, which are not considered for MOS here, the distribution pattern of MOS agrees with the other estimation of potential open-ocean macroalgae farming locations (e.g., Lehahn et al., 2016, Fig. 2a; Froehlich et al., 2019, Fig. 1). Another powerful limiting factor is the ocean surface temperature (Garcia et al., 2010, WOA2009 Dataset) which is too warm for our idealized species in many places, i.e., temperatures are above the T_opt (20 ^∘C) and nearly reach T_max (35 ^∘C).

Figure 2 Annual vertically integrated macroalgae biomass of MOS (a) and MOS-AU (b) in year 2024. Solid red lines outline the MOS-occupied area at year 2024 in both, while dashed red lines outline the initial MOS seeding area at year 2020 in (a). The simulated MOS area generally covers the NO₃-rich ocean surface (a) and can be expanded with nutrients supplemented by AU (b), making it larger than the estimated adequate area for macroalgae in previous studies Lehahn et al., 2016; Froehlich et al., 2019). Results for Areas I (blue circle), II (yellowish pentagon) and III (cyan rectangle) are discussed in the text and displayed in Fig. 3. Braces indicate the belts of N, M and S with various seeding strategies of macroalgae (Table 2), which are designed to avoid winter periods.

Figure 3 Vertically integrated macroalgae NPP simulated by experiment MOS (solid lines) for year 2024 and with representative Areas I (dark blue circle), II (yellow pentagon) and III (cyan rectangle) highlighted by rectangles with corresponding colors in Fig. 2. The NPP of macroalgae of experiment MOS_AU (dashed) shows an enhancement of NPP as expected.

At the beginning of year 2020, a surface area S_MOS of 72×10⁶ km² was selected by the MOS algorithm according to the requirements described in Sect. 3. This is equivalent to a total cultivated rope length (L_MOS) of 7.2×10⁶ km (Eq. A1, Sect. A1). When the macroalgae start to grow and to consume nutrients, regions with nutrient levels insufficient for further growth are gradually abandoned. By the end of year 2024, the MOS coverage has declined by about 3 % to 69.6×10⁶ km² (Table 4).

Despite the similar distribution patterns, MOS-occupied area (69.6×10⁶ km², ∼19.7 % of the world ocean) is larger than the assessments of ∼10 % of the world ocean by Lehahn et al. (2016) and $\sim \mathrm{48}\times {\mathrm{10}}^{\mathrm{6}}$ km² by Froehlich et al. (2019). Compared to the static estimation based on historical nutrient levels and temperature suitability (Lehahn et al., 2016, Fig. 2a), the dynamic processes redistributing nutrients, as well as the explicit macroalgae growth module in our simulations, contribute to the simulated larger potential area for MOS, especially in the equatorial Eastern Pacific and the Southern Ocean (Table 3). Besides, the adequate area for macroalgae cultivation was limited to the Economic Zones (EEZs) delineated by Froehlich et al. (2019) due to limitations of cost and political feasibility. This constraint has been ignored in the current study. Thus, our simulated MOS-adequate area is 45 % larger than the estimation by Froehlich et al. (2019).

4.1.2 Macroalgae model validation

Validation of the macroalgae model is crucial, as the productivity and macroalgal biomass yield are vital for CO₂ sequestration. Here, we examine the simulated seasonality, NPP rate and biomass yield of MOS in comparison with available observations and assessments. Simulated NPP is high in the first year of deployment in many regions because nutrients are abundant, and then it sharply declines in the following years as a new local biogeochemical state is reached. Thereafter, NPP gradually reaches a relatively steady state by 2024 (Fig. A2). To provide some validation of the macroalgae model, we select three areas named Area I, II and III from Belt N, M and S and analyze their performance in year 2024 (Fig. 3). Each area covers four grid boxes in the uppermost ocean layer of the UVic ESCM.

According to Sect. 3.1 and Table 2, the seeding date for MOS is 1 May in Area I, 1 January in Area II and 1 November in Area III, while the sinking dates are 31 October, 31 December and 30 April the next year, correspondingly. As a result, shown in Fig. 3, the macroalgae NPP in Area I peaks around September with the accumulation of macroalgae biomass in that area. In Area II, due to the nutrient limitation and nutrient competition with ambient phytoplankton, macroalgae biomass grows slower than in the other two areas, leading to a later peak of NPP around October. In Area III, which is located in the Southern Ocean where the nutrients are rich, the macroalgae NPP peaks in February. These results indicates a plausible seasonality of our macroalgae model.

In our simulations, simulated macroalgae NPP rates are comparable to the observed ranges in the productive areas that we selected here. Observed wild macroalgae NPP varies widely, ranging from 91 to 750 gC m⁻² yr⁻¹ (Krause-Jensen and Duarte, 2016). Our model reproduces the macroalgae NPP of 159.2–199.3 gC m⁻² yr⁻¹ in the selected areas (Table 4). Simulated biomass yields in these areas are in the previously reported range as well. Reports of the biomass yield of aqua-cultured Laminaria saccharina (now regarded as a synonym of Saccharina latissima) range from 40 t DW km⁻² yr⁻¹ in an offshore cultivation experiment by Buck and Buchholz (2004) to 456 t DW km⁻² in a coastal cultivation experiment by Peteiro et al. (2014). In our simulations, the yield of selected areas ranges from 492.4 to 648.2 t DW km⁻² yr⁻¹. The selected Area I yields 648.2 DW km⁻² yr⁻¹. In regions with similar latitudes to that of Area I, the biomass yield of aqua-cultured Saccharina japonica (formerly classified as Laminaria japonica) was ∼300 t DW km⁻² yr⁻¹ in China (Zhang et al., 2016) and reached 7280 t DW km⁻² yr⁻¹ in Japan (Yokoyama et al., 2007). Nevertheless, some simulated low values from the globally averaged and latitudinal belt-averaged results are not surprising considering that the open ocean tends to be more nutrient limited than coastal or near-shore regions where the aforesaid observed macroalgae NPP was measured. Our results provide some confidence that our idealized model can simulate macroalgae well enough with respect to typical biomass yield, seasonality and geographical distribution.

Table 4 Properties of globally implemented MOS. Selected areas are from data of year 2024, whereas belt areas are values averaged from 2020 to 2024. Areal NPP rates and “Biomass yield” refer to the respective MOS area. The observational data come from the references.

^a Krause-Jensen and Duarte (2016); ^b Buck and Buchholz (2004); Peteiro et al. (2014); ^c Zhang et al. (2016); Yokoyama et al. (2007)

4.2 Evaluation of MOS with artificial upwelling (AU)

As expected, AU increases the area occupied by MOS from 69.6×10⁶ km² in the run without AU to 129.6×10⁶ km² in the run with AU (Fig. 2). Obvious expansions of areas with suitable growing conditions are found in the eastern tropical Pacific and the North Atlantic. AU also expands S_MOS to the Indian Ocean, which was almost abandoned in regular MOS simulations. In Areas I, II and III, both NPP rate and biomass yield are enhanced due to the upwelled nutrients (belt column, Table 4 and Fig. 3). A closer look into the belt N, M and S areas shows that both the NPP rate and biomass yield per square meter of the deployment area decrease in simulation MOS_AU when compared to the standard MOS simulation (belt column, Table 4). This is related to a “dilution effect”: in the new adequate areas made accessible for MOS by AU, the available nutrients are limited, thus the MOS NPP is relatively low compared to the original nutrient-rich NPP areas. Despite of this, the expanded MOS area in MOS_AU increases the total CO₂ captured by about ∼30 % (Table 4).

4.3 MOS deployment until year 2100

This section showcases the CDR and climate change mitigation capacities of MOS within the 21st century. Impacts of MOS on marine biogeochemistry (nutrients, dissolved oxygen and pelagic ecosystem) and global carbon cycles will also be examined.

4.3.1 CDR & climate change mitigation capacities

Over the 80 years between year 2020 and year 2100, MOS will mainly be deployed in nutrient-rich regions such as the Southern Ocean and the northern and eastern equatorial Pacific (Fig. B1a), although some contraction of initially occupied areas occurred due to the removal of nutrients.

By the year 2100, MOS (MOS_Stop and MOS_Conti) has sequestered 270 PgC (990 Pg CO₂; Table 5), representing ∼37 % of the cumulative CO₂ emissions in the RCP 4.5 pathway. Essentially all of MOS-captured carbon is retained in the ocean over this period as either remineralized dissolved inorganic carbon or organic carbon in the sunken biomass.

The CDR capacity of MOS is sensitive to the MR_C:N (molar C:N ratio of macroalgal biomass; Table 1). Compared to MOS_Conti, the carbon captured by MOS (MOS-C) in MOS_Conti_CN_High is raised by 22 % (by the year 2100) and by 19 % (by the year 3000) when the MR_C:N is increased by 20 %, to 24. When the MR_C:N is 20 % lower than the original value, MOS-C decreases by 13 % (by the year 2100) and by 18 % (by the year 3000). Our results agree with the range of CDR potential reduction by nutrient reallocation (7 %–50 %) reported in Bach et al. (2021).

Table 5 Model simulations under the RCP 4.5 emission scenario. MOS-C represents the carbon sequestered via MOS. C_atm, C_oc and C_ter stand for atmospheric, oceanic and terrestrial carbon reservoirs respectively. ΔSAT stands for surface-averaged temperature relative to 13.18^∘C, the pre-industrial.

In the model, MOS thus gradually reduces atmospheric CO₂ and thereby also limits global warming with respect to the pre-industrial period (ΔSAT; Fig. 4); i.e., the temperature increase of 2.14 ^∘C by the year 2100 is 0.38 ^∘C lower than ΔSAT of Control_RCP 4.5 but is still missing the 2 ^∘C target.

When AU is deployed in conjunction with MOS, the CDR capacity and mitigation effects of MOS are enhanced (Figs. 4a, c, 5). By the end of year 2100, 446.8 Pg carbon is sequestered by MOS_AU, an increase of 39.5 % relative to normal MOS. Correspondingly, MOS_AU successfully achieves the 2 ^∘C target of the Paris Agreement by maintaining a ΔSAT at 1.89 ^∘C relative to pre-industrial (Fig. 4c and Table 5). As in the run without AU, essentially all of the carbon captured via MOS is stored in the ocean until the end of the 21st century (FR, Table 5).

Figure 4 Simulations of (a, b) annual global mean atmospheric CO₂ concentration and (c, d) surface-averaged temperature relative to the pre-industrial (average of year 1850 to year 1900) level of 13.18 ^∘C (ΔSAT). Under RCP 4.5 scenario, MOS reaches the 2 ^∘C target in conjunction with AU, while the 1.5 ^∘C cannot be met in all MOS simulations. Note that MOS is terminated whenever pre-industrial concentrations of atmospheric CO₂ are reached, as seen for MOS_AU_Conti (solid orange) and MOS_Conti_NoRe (dotted blue) in (b, d). Both atmospheric CO₂ and ΔSAT remain lower than control after MOS termination.

Figure 5 Temporal evolution of globally integrated sunken macroalgal biomass on the seafloor. Biomass generally increases with fertilization by AU. In the idealized zero-remineralization simulations, all sunken macroalgal biomass remains on the seafloor, and globally integrated sunken macroalgal biomass shows a monotonous increase.

4.3.2 Global carbon cycle impacts

The net effect of the MOS-induced climate–carbon cycle perturbation is an increase of the oceanic carbon reservoir (C_oc) and a decrease of the atmospheric and terrestrial carbon reservoirs (C_atm, C_ter). MOS enhances oceanic carbon uptake by increasing the atmosphere to ocean carbon flux (Fig. B11), which is driven by the DIC removal by MOS in the oceans' surface layers. However, the terrestrial carbon reservoir declines (relative to Control_RCP 4.5) in all MOS simulations (Table 5). The atmosphere-to-land carbon flux is reduced in MOS simulations (Fig. B10). One cause is the photosynthesis that is reduced by lower CO₂ fertilization of land biota (Keller et al., 2018). This result is in line with other studies showing that CDR can lead to a weakening and even reversal of natural carbon sinks (Keller et al., 2018). For instance, by the year 2100 (Table 5), due to the declined terrestrial carbon pool (C_ter, −29.7 PgC), the reduction in atmospheric carbon pool (C_oc, −142.6 PgC) is less than the gain in the ocean (C_oc, 171.8 PgC). Besides, it is also worth noting that the increment of C_oc in MOS and MOS_AU is 171.8 and 283.9 PgC, which is less than the cumulative amount of carbon sunk out of the surface layer via MOS by year 2100 (Table 5). One reason is that the reduced oceanic carbon uptake by declined PNPP (phytoplankton net primary production; Sect. 4.3.4) offsets the MOS-induced carbon sequestration. As shown in Fig. 9, by the year 2100, global PNPP has been reduced by 20 %, while POC export has been reduced by 30 % in MOS_Stop.

MOS also impacts the distribution of DIC in the ocean. The DIC profiles in Fig. 6 illustrate that the general effect of MOS is to move more DIC to a greater depth (z≥3000 m). By the end of year 2100, MOS simulations show an increased total DIC concentration in the deeper oceans when compared to Control_RCP 4.5 (except for the zero-remineralization sensitivity runs discussed below). For instance, in the deep Southern Ocean, the DIC concentration is, on average, nearly 80 µmol kg⁻¹ higher than the Control_RCP 4.5 in year 2100 (Fig. 6). The conjunction of MOS with AU increases average deep-ocean DIC even more. An example is the simulated increase of DIC in the deep Pacific Ocean and Atlantic Ocean basins by MOS_AU in year 2100 (orange line in Fig. 6, DIC panel). In contrast, DIC concentrations are reduced in shallower waters (depth <1000 m), as the air–sea carbon flux is unable to fully compensate for the carbon removal by MOS. The current model results provide additional evidence that the CDR potential of MOS is partly offset by its negative impacts on the pelagic biological production and the biological carbon pump. In an additional model run (not shown) without MOS but with CO₂ emissions reduced by the annual equivalents of the MOS-induced carbon exports, yielding a total amount of 270 PgC by year 2100, the 270 PgC emissions removal yields a reduction of atmospheric CO₂ by 171.5 GtC by year 2100. This reduction is 20 % higher than the atmospheric CO₂ reduction of 142.6 PgC realized in the original MOS simulation where the MOS-induced shading and removal of nutrients from the surface layers reduce the biological carbon pump and the associated carbon storage in the ocean. When CO₂ emissions are instead cut by an amount corresponding to 80 % of the MOS-induced carbon export, atmospheric CO₂ concentrations simulated by the MOS-free emission-cut runs agree closely with those of the respective MOS experiments. That is, each tonne of CO₂ sequestered in the ocean by MOS is, in our model and on a 100-year timescale, equivalent to an emission cut of about 0.8 t of CO₂.

Figure 6 Global and basin-wide averaged vertical profiles of various model tracers in year 2100 and year 3000 under the RCP 4.5 emission scenario. In general, MOS (except for the zero-remineralization one) transports DIC and nutrients in the surface layer to the deep ocean. The oxygen levels are increased in the mid layers due to the declined downward organic-particle flux (Sect. 4.3.4), but they are decreased in the deep ocean as a result of the remineralization of sunken biomass. These impacts are strengthened when MOS is deployed continuously and/or in conjunction with AU.

4.3.3 Impacts on global nutrients distributions

By the year 2100, the deployment of MOS has changed the global patterns of NO₃ and PO₄. At the surface, NO₃ and PO₄ concentrations decrease due to MOS nutrient consumption. In the deep ocean (depth≥3000 m), PO₄ and NO₃ increase due to the remineralization of sunken macroalgal biomass (except for in the MOS_NoRe simulations). The largest increase in deep-ocean PO₄ appears in the Southern Ocean, while the smallest increase is found in the Indian Ocean (PO₄ yr2100 groups in Fig. 6). This is caused by the distribution of MOS in the surface layer, which, in our simulations, occupies large areas in the Southern Ocean but only a relatively small region in the Indian Ocean (Fig. 2).

The remineralization of sunken biomass consumes dissolved oxygen and releases NO₃ and PO₄. However, low-oxygen environments and the associated switch from aerobic remineralization to denitrification occupy relatively small areas so that this impact is not easily detectable in global nutrient profiles.

In addition to the localized depletion of nutrients by MOS, the MOS-induced Southern Ocean uptake and transport of N and P to the deep ocean act as a type of “nutrient trapping” (Fig. B2). These dynamics thereby reduce nutrients and productivity in middle to low latitudes because less N and P are available to be transported out of the Southern Ocean. A similar dynamic has been seen in modeling studies of ocean iron fertilization (Oschlies et al., 2010a; Keller et al., 2014).

4.3.4 Impacts on simulated pelagic ecosystems and the organic particle export

In our simulations, large-scale deployment of MOS has an impact on pelagic ecosystems, mainly on phytoplankton NPP (PNPP) and biomass.

In the MOS simulations (MOS_Conti/Stop), globally integrated annual PNPP decreases by 20 % (9.5 PgC yr⁻¹) by year 2100 (Table 5). One reason is the canopy shading effect of the floating macroalgae farms, which reduces downward solar radiation available for the phytoplankton community below. In addition, there is nutrient competition between macroalgae and phytoplankton. As shown in Fig. 7a, by the end of the 21st century, PNPP will declined in MOS areas, e.g., the northern and eastern equatorial Pacific and the Southern Ocean. Intriguingly, in a few regions outside the MOS deployment region, PNPP has increased instead. For instance, a “halo” of enhanced PNPP can be observed surrounding the eastern equatorial Pacific MOS region (Fig. 7a). Similar circumstances are simulated in the North Pacific, in the Southern Ocean (60^∘ E:120^∘ E, 30^∘ S) and off the equatorial west coast of Africa. This PNPP enhancement is sustained by the outflow of residual nutrients from MOS deployment regions (see Fig. B9). One reason is that the macroalgae growth is constrained by the maximum biomass yield, as described in Sect. 2.2.1. Macroalgae nutrient uptake thus cannot compensate for the loss of nutrient consumption by light-limited PNPP within the MOS region, which results in enhanced surface nutrients compared to the simulation without MOS, especially when AU supplies additional nutrients to the surface.

Figure 7 Vertically integrated annual PNPP in year 2100. (a, b) MOS minus Control_RCP 4.5, with red boundaries contouring the MOS-occupied area; (c) Control_RCP 4.5; (a) illustrates a decline in PNPP in MOS-occupied areas accompanied by a “halo” of enhanced PNPP surrounding MOS areas, particularly in the ETP caused by the leakage of residual nutrients (Sect. 4.3.4). These impacts on PNPP are amplified in MOS_AU (b).

In the MOS_AU simulations, the PNPP “halo” can be seen in almost the entire MOS-free ocean surface (Fig. 7b). The AU fertilization effect enhances the nutrient leakage from the MOS area. This leads to a higher PNPP in MOS_AU than in the normal MOS simulations (Fig. 7c), but one that is still lower than in the Control_RCP 4.5 run (Fig. 9a1).

Changes in the global particulate organic carbon (POC) export flux generally follow the pattern of PNPP changes (Fig. 9b1). Thus, when MOS is present, the PNPP reduction results in a weakened POC flux.

Figure 8 Proportion of MOS NPP in the global oceanic NPP by year 2100 ($\mathrm{MOS}\mathrm{\_}\mathrm{NPP}/(\mathrm{MOS}\mathrm{\_}\mathrm{NPP}+\mathrm{PNPP})\times \mathrm{100}$). Note that the NPP values are converted to carbon using the respective C:N ratio. The MOS_NPP generally amounts to more than 70 % of total oceanic NPP where MOS is deployed, indicating an obvious NPP shift from phytoplankton PNPP to MOS_NPP.

Figure 8 illustrates the shift of oceanic NPP from PNPP to MOS_NPP. In regions where MOS is deployed, 70 % of the total NPP is macroalgae NPP. The macroalgal NPP is thus nearly twice as high as PNPP. This may lead to additional ecological and biogeochemical issues. One of them is the decline of zooplankton led by the reduced PNPP in this study. We performed an additional simulation, in which the zooplankton grazing on MOS is turned off, and the grazing preferences follow the original settings in Keller et al. (2012). As shown in Fig. B12, the grazing by zooplankton on MOS has no significant effect on either the zooplankton biomass or the MOS_NPP. As the zooplankton grazing preference for macroalgae is lower than for phytoplankton, the zooplankton community is still mainly fed by phytoplankton. Therefore, the decline in zooplankton biomass (Fig. B4) follows the declining phytoplankton biomass trend (Fig. B6).

Figure 9 Temporal evolution of globally integrated PNPP (a1, a2); downward POC flux at 2 km depth (b1, b2). The termination runs branch off of the continuous ones in year 2100 and are identical up to that point. Through the 21st century, MOS reduces PNPP and POC export due to canopy shading and competition for nutrients. Obvious rebounds followed by quick decline can be observed right after terminations of MOS.

4.3.5 Impacts on dissolved oxygen

The two major impacts of MOS on oceanic dissolved oxygen are as follows: (1) increased deoxygenation at the seafloor by the remineralization of sunken macroalgal biomass (except for in MOS_NoRe) and (2) increased dissolved oxygen at middle depths (e.g., 300 m depth) caused by the reduction of the downward POC flux and the associated decline in oxygen consumption by POC remineralization.

In Control_RCP 4.5, the global oceanic dissolved oxygen inventory decreases throughout the simulation. The two main driving mechanisms are the reduced solubility in the warming ocean and the decelerating overturning circulation. The long-term decline of oxygen is especially obvious at depth (1200 m), which is induced by increasing deep-water residence times and the accumulation of respiratory oxygen deficits under global warming (Oschlies et al., 2019; Oschlies, 2021).

As a result of reduced respiration in the upper water column, the size of the oxygen-minimum zone (OMZ) in the eastern tropical Pacific (ETP) shrinks substantially, and the volume of water with O₂<80 mmol m⁻³ in the North Pacific even disappears. In the Southern Ocean, dissolved oxygen increases as well (Fig. 10c). This is more pronounced when AU is applied (Fig. 10e). The increase in dissolved oxygen is caused by decreased microbial remineralization of POC, a consequence of the reduced downward POC flux resulting from the inhibition of PNPP in the surface layer (Sect. 4.3.4). Some decrease in oxygen concentrations occurs in the western Pacific and the Indian Ocean (Fig. 10g), where the surface PNPP is enhanced by the surplus nutrients that leaked out from the MOS-occupied area (see Sect. 4.3.4 and Fig. 7c).

Figure 10d, f and the O₂ yr2100 panel of Fig. 6 illustrate how MOS changes dissolved oxygen in the deep ocean. Within the normal MOS simulations, the decline of benthic dissolved oxygen has mainly happened in the Southern Ocean by year 2100, with the appearance of a few new areas with oxygen concentrations of less than 80 mmol m⁻³ (Fig. 10d). However, when AU is also deployed, the increased macroalgal biomass sinking and remineralization creates even more benthic low-oxygen zones (Fig. 10f) in the ETP and North Pacific Ocean. These new locations correspond to MOS-occupied surface areas.

Figure 10 Dissolved O₂ concentration distribution in year 2100 at 300 m depth (a, c, e) and at the seafloor (b, d, f); lines delineate boundaries of OMZs at any location within the water column with less than 80 mmol m⁻³ oxygen (dashed white) and less than 20 mmol m⁻³ oxygen (solid red). At 300 m depth, elevated dissolved-oxygen levels in MOS simulations are caused by the decline in POC export (c); exceptions are the reduced oxygen concentrations in regions outside the MOS_AU deployment (e). On the seafloor, remineralization of sunken biomass creates several new low-oxygen areas (d, f).

4.4 Long-term effects of MOS

Here, we will analyze the long-term effects of hypothetical massive MOS deployment beyond the Paris Agreement timeframe on a millennial timescale.

Even after a simulated continuous millennial-scale deployment, the distribution of MOS in year 3000 is nearly identical to the one in year 2100 with only a minimal decrease in biomass (Fig. B1b). When deployed beyond the year 2100 (MOS_Conti), MOS will continue to sequester carbon and reduce atmospheric CO₂ on millennial timescales or, in our setup, until atmospheric CO₂ falls back to the pre-industrial level of 280 ppm. The MOS_Conti simulation ultimately sequesters 2533 PgC and decreases atmospheric CO₂ to 318.5 ppm CO₂ by the year 3000 but never achieves the pre-industrial CO₂ level. Notably, atmospheric CO₂ stops decreasing by year 2780 and rebounds afterwards even though MOS continues to sequester carbon. This can be explained by a recurrent deep convection in the Southern Ocean around year 2800 that accelerates oceanic carbon leakage back to the atmosphere (Martin et al., 2013; Reith et al., 2016; Oschlies, 2021). Meanwhile, the leakage of MOS-captured carbon eventually offsets the MOS carbon sequestration (Sect. 4.6).

In the sensitivity simulations MOS_Conti_NoRe and MOS_AU_Conti, atmospheric CO₂ reaches 280 ppm by the years 2820 and 2475, respectively. After reaching 280 ppm, MOS is stopped, and atmospheric CO₂ increases again as remineralized carbon leaks out of the ocean and as the surface ocean adjusts to the no-MOS situation. The largest increase in CO₂ is found in MOS_AU_Conti. Meanwhile, when MOS is deployed (uninterruptedly or till the CO₂ 280 ppm trigger), the land carbon uptake is constantly lower than the control level owing to the reduced CO₂ fertilization effect. Due to the permanent storage of MOS-C in sunken biomass, rebounds of atmospheric CO₂ are relatively gentle in MOS_Conti_NoRe (Fig. B11). Nevertheless, the atmospheric CO₂ levels in continuous MOS simulations are significantly lower (35 % to 50 % of Control_RCP 4.5) by the end of year 3000.

The side effects of MOS also persist and often grow in magnitude with continuous deployment. Though PNPP is enhanced around MOS areas by means of nutrient leakage (PNPP “halo”; see Sect. 4.3.4), the global reduction of surface nutrients and local canopy shading by MOS leads to continuous but gentle lowering of global PNPP after the sharp decreases in the initial 20 years (Fig. 9a2). For instance, in MOS_Conti, PNPP drops by ∼60 % by the end of year 3000 (Table 5). Correspondingly, in MOS_Conti POC export eventually declines by 50 % relative to Control_RCP 4.5 (Fig. 9b2). In sensitivity run MOS_AU_Conti, the nutrient supply by AU, which initially maintains a higher phytoplankton biomass and NPP than in MOS without AU (Sect. 4.3), declines with time as source waters of the upwelling become reduced in nutrients. Therefore, PNPP and the POC export levels drop after year 2200 (Fig. 9a2, b2).

The redistribution of DIC and nutrients is intensified in the continuous simulations. As shown in Fig. 6, when remineralization of MOS sunken biomass is turned on, the Pacific deep ocean and the Southern deep ocean show the highest DIC and PO₄ enrichment by year 3000. The accumulations depend on ocean circulation (e.g., thermohaline circulation) and the distribution of MOS at the surface. In MOS_Conti_NoRe, the ocean DIC decreases globally relative to Control_RCP 4.5. This results from the continuous DIC removal into biomass via MOS with no remineralization. Another cause of the declined global DIC is the declining downward POC flux owing to the PNPP reduction caused by the declining nutrient levels in the surface layer (see Sect. 4.3.4). NO₃ enrichment in the deep ocean is considerably smaller than that of DIC and PO₄ because of enhanced denitrification in the developing benthic low-oxygen regions (Sect. 4.3.5, Fig. B3). In the zero-remineralization situations, deep-ocean PO₄ and NO₃ concentrations decrease compared to the control levels due to the reduced remineralization of POM resulting from the weakened downward flux of POM.

As shown in Fig. 11 and the O₂ yr3000 panel of Fig. 6, dissolved oxygen concentrations at middle depth (e.g., 300 m) increased during millennial MOS deployment due to reduced PNPP and associated downward flux and remineralization of POM in the water column. In benthic waters, regions with very low dissolved oxygen are shown in Fig. 11f, j in the Pacific and Southern oceans. In contrast, increased oxygen concentrations are found in MOS_Conti_NoRe (Fig. 11h), especially in the Atlantic, the Indian and the Southern oceans. Besides the absence of oxygen consumption by macroalgal biomass remineralization, another reason for these oxygen increases lies in the reduction of POC downward flux described in Sect. 4.3.4.

Figure 11 Dissolved O₂ concentrations at depth ∼300 m (left panel) and the ocean bottom (right panel) in year 3000: contour lines indicate boundaries of OMZs with less than 80 mmol m⁻³ oxygen (dashed white) and less than 20 mmol m⁻³ oxygen (solid red). Continuous MOS deployment further shrinks the OMZ at 300 m depth (e, g, i) but expands them at the bottom (f, j), except for the zero-remineralization one in (h).

4.5 Termination effects

After termination of MOS_Stop in year 2100, the atmospheric CO₂ concentrations and SAT (surface-averaged temperature) both rise but generally remain lower than those of the Control_RCP 4.5 simulation (Fig. 4b, d). More than half (FR ranges from 58.6 % to 64.4 %, Table 5, calculated by Eq. 23) of MOS-captured carbon is still stored in the ocean by the end of year 3000. As shown in Fig. 4b, d, in the termination simulations MOS_Stop and MOS_AU_Stop, CO₂ concentration and SAT gradually converge against the Control_RCP 4.5 run as a result of DIC from remineralized macroalgal biomass being transported to the ocean surface and into the atmosphere (see subsequent Sect. 4.6). By year 3000, the atmospheric CO₂ in MOS_Stop is only 28.5 PgC less than in Control_RCP 4.5, while ΔSAT rebounds slightly from −0.38 to −0.23 ^∘C. In MOS_AU_Stop, the differences between pCO₂ and ΔSAT are smaller than the normal MOS, as AU has augmented MOS carbon sequestration, and 64.4 % of it is retained in the ocean. As expected, the idealized non-remineralization condition (MOS_Stop_NoRe) is able to permanently store the sequestered carbon; thus, the rebounds of CO₂ and ΔSAT are less than the normal MOS_Stop.

In all MOS termination simulations, PNPP and POC export rebound abruptly following the cessation of MOS, but they sharply drop over the subsequent decades (Fig. 7a1). The sharp increase in PNPP and POC export results from the sudden absence of macroalgae as a main competitor for nutrients and light. The subsequent decline in PNPP and POC export results from the consumption of the surface nutrients and the lack of subsurface nutrients that had previously been exported directly to the seafloor with the sinking of macroalgae biomass (Fig. B7). Afterwards, PNPP recovers gradually due to the slow returning of remineralized nutrients to the upper ocean. By the year 3000, PNPP in MOS_Stop recovers to 97 % of the control level (Table 5), with the only differences being attributable to the slightly different climate state. In the MOS_NoRe simulation, the PNPP recovery is slower due to the permanent nutrient removal from the upper water column. In the MOS_AU_Stop simulation, PNPP rebounds to higher levels than the normal MOS but drops to lower levels afterwards. The amplified oscillation of PNPP results from the simultaneous termination of AU and MOS: when MOS_AU_Stop is suddenly terminated, the canopy shading and nutrient competition by MOS are removed. Meanwhile, the surplus of nutrients from AU still remains. This boosts PNPP rapidly. However, once these nutrients are consumed, the natural nutrient supply to surface waters is insufficient to maintain the high PNPP.

After the termination of MOS, the rate of oceanic carbon uptake falls abruptly (Fig. B11). After a short peak caused by the abrupt rebound of PNPP, it remains slightly lower than the control level due to the declined PNPP rates and lower atmospheric CO₂ levels. Oppositely, the MOS-induced reduction in terrestrial carbon uptake starts to rebound after MOS cessation (Fig. B10) due to the rise of atmospheric CO₂, which tends to enhance the terrestrial photosynthesis.

When MOS deployment is stopped, the elevated (compared to Control_RCP 4.5) dissolved-oxygen concentrations at middle depth generally decline as the downward POC flux recovers (O₂ yr3000 panel in Fig. 6). The lowered oxygen concentrations in the deep ocean are also reversible after cessation of MOS. For instance, by year 3000, the benthic dissolved oxygen of MOS_Stop (Fig. 10d) is similar to that of Control_RCP 4.5 (Fig. 10b).

4.6 Leakage of MOS-sequestered carbon

The leakage of MOS-sequestered carbon (MOS-C) occurs mostly in the Southern Ocean (Fig. 12, B5). The explanation lies in the dynamics of Southern Ocean upwelling, where Pacific deep water (PDW), Indian deep water (IDW) and Antarctic deep water (AADW), laden with DIC of remineralized MOS biomass, reach the subantarctic ocean surface (Talley, 2013; Weber and Bianchi, 2020; Anderson and Peters, 2016). Moreover, the recurrent deep convection (see Sect. 4.4) in the Southern Ocean around year 2600 accelerates the carbon leakage, which can be observed in Fig. 12 as an enhanced outgassing around year 2600.

The outgassing of MOS-C in discontinuous simulations (e.g., MOS_Stop) starts from year 2100, while the continuous ones starts from year 2300 (Fig. 12). Thus, by the end of the 21st century, nearly all the MOS-C in all MOS simulations is retained in the ocean. However, by the year 3000, even in the continuous MOS simulation, only about 75 % of MOS-C remains in the ocean, while the accelerated vertical water transport by AU slightly reduces this portion to 73 %. In run MOS_Stop, 59 % of MOS-C remains in the ocean by year 3000, whereas the additional sunken macroalgal biomass in MOS_AU_Stop results in more MOS-C (64 %) being retained (Table 5). When sunken biomass is free from remineralization (NoRe runs), the contained carbon is permanently isolated from the atmosphere and stored in the ocean.

Figure 12 Zonally averaged MOS-captured carbon outgassing in the simulation (a) MOS_Conti and (b) MOS_Stop. When conveyed back to the surface, the DIC from MOS remineralization participates in the air–sea exchange (Sect. 2.2.5). MOS-C outgassing starts in year 2100 when MOS is terminated (a) or after year 2300 when MOS is continuously deployed. The outgassing mainly happens in the Southern Ocean. The outgassing is strengthened around year 2800 when a Southern Ocean deep-convection event accelerates the upwelling of deep waters with high concentrations of remineralized DIC (Sect. 4.4).

5 Concluding discussions

In this study, we have tested the potential of the “macroalgae open-ocean mariculture and sinking (MOS)” as a carbon dioxide removal method. Although environmental conditions (e.g., nutrients, temperature) in the open oceans differ considerably from the coastal and/or near-shore regions where macroalgae aquaculture is currently applied in reality, our simulations suggest that, in certain open-ocean regions, macroalgae may successfully grow and sequester carbon (if engineering constraints can be overcome).

Even for continuous deployment at a maximum scale currently deemed possible, MOS alone is not able to reduce the warming to the 2 ^∘C target by the end of the 21st century under the RCP 4.5 moderate-mitigation scenario. This finding is consistent with conclusions from previous studies that no single carbon dioxide removal (CDR) method alone can ensure that the current climate goals are reached (Keller et al., 2014; Anderson and Peters, 2016; Fuss et al., 2020; IPCC, 2018). Clearly, emissions reductions must be the primary means of mitigation; a portfolio of CDR options can only complement these other efforts.

In order to estimate the maximum CO₂ removal potential, the possible synergy of MOS with artificial upwelling (AU) has been investigated. The employed AU concept aims at piping nutrient-rich deep water to the surface to enhance the growth of macroalgae in MOS. As expected, AU is found to have the potential to successfully enlarge the growing area of MOS and to enhance the CDR capacity of MOS.

In the first 80 years of deployment, the maximum MOS carbon sequestration potential is 3.38 PgC yr⁻¹ for regular MOS, but it can be boosted up to 5.56 PgC yr⁻¹ with assistance from AU. If deployment is discontinued from year 2100, about 58.6 % to 70.2 % (normal MOS and MOS_AU, respectively) of MOS-sequestered carbon would be retained in the ocean by year 3000.

Several potential side effects have also been revealed and analyzed. One side effect is the reduction in phytoplankton NPP (PNPP) due to canopy shading and nutrient removal from the sea surface to the bottom by MOS. The declined PNPP in turn offsets ∼37 % of the MOS CDR. Intriguingly, some areas with enhanced PNPP (PNPP “halo”, Sect. 4.3.4) are found surrounding the major areas occupied by MOS, fueled by the residual NO₃ that leaks from MOS areas.

Another strong side effect of large-scaled MOS is the impact on oxygen distributions. Dissolved-oxygen concentrations increase in near-surface and intermediate waters, as MOS reduces the downward flux of plankton-derived organic matter by restraining the surface PNPP. On the other hand, the massive amount of sunken biomass from MOS at the ocean bottom and its subsequent remineralization consumes oxygen and can create large benthic OMZs.

An uncertain factor is the fate of the sunken biomass. It will affect the benthic fauna by depositing large amounts of organic matter, as well as by expanding low-oxygen regions on the seafloor upon oxygen consumption by remineralization. Therefore, we performed additional sensitivity simulations focusing on the macroalgal biomass remineralization rate. When macroalgal biomass does not undergo microbial remineralization, the captured CO₂ can be permanently stored without leakage. This increases the CDR potential of MOS. The benthic OMZs created by remineralization of sunken biomass would also be avoided, while the shrinking of intermediate water OMZs persists. However, other side effects cannot be neglected: in zero-remineralization simulations, the constant removal of nutrients in the surface will impede the recovery of PNPP. This may eventually affect the marine surface ecology and ocean services such as food provision. These potential side effects are also noteworthy for another macroalgae farming approach, i.e., harvesting the macroalgae biomass for bioenergy with carbon capture or storage (BECCS) or biochar (Sect. 3.2).

The impacts of MOS on oxygen distributions may also influence the oceanic sources of N₂O, an atmospheric GHG gas and a major ozone-depleting compound (Ravishankara et al., 2009). The increased and decreased oxygen levels in the middle and bottom layers impact denitrification and may weaken the N₂O sources in the subsurface but increase those in the deep waters (Bange et al., 2019).

Moreover, attention should also be paid to the calcification by calcareous macroalgae (if cultivated) and/or associated epibionts that grow on macroalgae. Bach et al. (2021) have suggested that epibionts living on Sargassum offset 16.5 % of the POC fixed by Sargassum and therefore decrease its natural carbon sequestration potential if the biomass was intentionally sunk for CDR purposes. These calcification rates and the response to ocean acidification of macroalgae are also species specific (Koch et al., 2013). These factors need to be investigated with further research if macroalgae are to be considered for ocean-based CDR methods.

The production and export of DOC is also an area of further study for large-scale farmed macroalgae for carbon sequestration. Macroalgae have been reported to release considerable amounts of DOC and to contribute to the global DOC export from coastal to open-ocean waters. The estimated global DOC release of macroalgae habitats is 355 Tg C yr⁻¹ (Krause-Jensen and Duarte, 2016). The averaged DOC release rate by macroalgae is of 23.2±12.6 mmol C m⁻² d⁻¹ (equivalent to 8.5±4.6 mol C m⁻² yr⁻¹) but with a wide range of 8.4 to 71.9 mmol C m⁻² yr⁻¹ (Barrón et al., 2014). If we simply multiply this annual-averaged DOC release rate with the MOS-occupied area (S_MOS, Table 5), the estimated annual DOC release by MOS would be 7.1±3.8 PgC (MOS) or 12.9±7.0 PgC (MOS_AU). Although the refractory DOC released by macroalgae could potentially be an additional contribution of carbon sinking by MOS, the available information on the generation and composition of the macroalgae DOC is not enough to parameterize a model of this process, and more research on the topic is needed (Krause-Jensen and Duarte, 2016; Barrón et al., 2014; Barrón and Duarte, 2015).

Another side effect not investigated here is the production and emission of halocarbons from macroalgae farms. Macroalgae species have been reported to generate halocarbons in polar, temperate and tropical coastal regions, with a highest producing rate of 6000 pmol CHBr₃ gFW⁻¹ h⁻¹ (Leedham et al., 2013; Baker et al., 2001; Carpenter and Liss, 2000; Latumus, 1995). These volatile low-molecular-weight halocarbon compounds (e.g., CH₃I, CHBr₃ and CHCl₃) are potent greenhouse gases (Meinshausen et al., 2011). They also influence stratospheric ozone destruction when transported by deep atmospheric convection into the stratosphere (Ziska et al., 2013; Tegtmeier et al., 2012, 2013), therefore enhancing radiative forcing (Ramaswamy et al., 1992; Daniel et al., 1995). Large-scale MOS cultivation might release a significant quantity of halocarbons. However, as MOS also reduces global phytoplankton NPP, it is likely that the production of halocarbons by phytoplankton decreases. Further studies are needed to investigate the possible effects of halocarbon emissions from large-scale macroalgae cultivation and how this is offset by a potential decrease in phytoplankton halocarbon production.

Meanwhile, here, we only discuss the CDR potential of the MOS-AU combination under an idealized situation, where the artificial upwelling (AU) upwells nutrients without changing the ambient water temperature on the surface. The aim of deploying akin AU to MOS is to assess the maximum potential of CDR, as AU can upwell nutrient-rich deeper water. However, it can be expected that the NPP of MOS will be slower if the ambient water temperature is reduced by the upwelled cold water, which will further limit the CDR potential of MOS_AU. Moreover, several dominant side effects due to the upwelled cold water in ordinary AU have been revealed, such as the quick rebound or even surpassing of CO₂ concentrations and surface temperatures after the termination of AU (Oschlies et al., 2010b). If the cold-water upwelling was included, there would be extra analysis on the associated impacts on planetary radiation budget imbalance, marine biogeochemistry and global carbon pool (e.g., the enhanced terrestrial carbon sequestration due to the cold effects by AU), which are beyond the scope of this study. Thus, we used the hypothetical AU systems excluding thermal exchange to address the nutrients supplementation of AU to MOS and to avoid the unnecessary complexities and side effects in the meantime. However, further studies are required if ordinary AU should be considered in association with macroalgae farming.

Besides the CDR effect of MOS, in the case of large-scale deployment, the macroalgae farms are likely to increase the albedo of the oceans’ surfaces, especially when they occur near the sea surface (Fogarty et al., 2018; Bach et al., 2021). Meanwhile, we have not considered possible hydrodynamic impacts on ocean circulation, as the thick macroalgae layers and the farming infrastructures may influence the momentum and mixing of the ambient flow fields (Liu and Huguenard, 2020; Nepf, 2012; Thomas and McLelland, 2015). The MOS model analysis presented here clearly has some limitations that future studies might improve on. One of the most critical issues is to improve the realism of the model design by including more representative macroalgae species for various regions. Another aspect is the consideration of dynamic cellular stoichiometry of macroalgae. With a better simulation of the cellular quota, we could improve our understanding of the relation of nutrient and carbon fluxes to MOS and the environment. Explicit consideration of the variable morphology of the macroalgae, as well as of the impacts of currents on frond erosion (Broch and Slagstad, 2012), would also improve the representation of macroalgae loss rates in the model. Further optimization of deployment timing and location for MOS are achievable by evaluating data from field tests or implementations of macroalgae mariculture in the open oceans. Another aspect that needs improvement is the modeling of benthic macroalgal biomass remineralization. Here, we treated macroalgal biomass homogeneously as particulate organic matter. Though the degradation of macroalgal fragments under deep-sea conditions (e.g., low temperature and unique microbial colonies) remains unclear, it might be different from POC in terms of remineralization rate and oxygen consumption. Tracking of macroalgae-sequestered carbon will be required to record its fate and possible carbon leakage after sinking (e.g., by an eDNA method to trace macroalgae carbon in marine sediments, such as that by D'Auriac et al., 2021). The economic perspectives of developing and deploying MOS also need to be investigated. Furthermore, more research on how large-scale macroalgae mariculture will impact human activities (e.g., ocean shipping, fisheries) needs to be undertaken. Associated legal and political issues regarding the usage of international waters for MOS deployment should be considered as well.

Overall, this study adds to the rapidly expanding field of considering macroalgae cultivation for CO₂ removal. The evidence from this study suggests that macroalgae mariculture and sinking has a considerable CDR potential under these idealized conditions, but it brings about substantial side effects on marine ecosystems and marine biogeochemistry. Given this, the concept requires further research with less idealized experimental settings to determine if its CDR benefits outweigh the side effects (Dean et al., 2021).

Appendix A: MOS validations

A1 MOS yield calculation

For the convenience of calculation, we assume that, when MOS occupies a surface grid cell, the area is covered by parallel cultivation ropes (lines) with an interval distance of d (see Table 1 and Fig. A1). The total length of the cultivation lines (L_MOS, in meters) of MOS in the grid cell is then

$$\begin{array}{}\text{(A1)}& L_{\mathrm{MOS}}=S_{\mathrm{MOS}}÷d,\end{array}$$

where S_MOS (m⁻²) is the area of ocean surface occupied by MOS. Accordingly, the conversion between macroalgal biomass yield on ropes (Y_rope in kg DW m⁻¹) or in fields (Y_field in kg DW m⁻²) can be calculated as follows:

$$\begin{array}{}\text{(A2)}& Y_{\mathrm{rope}}=Y_{\mathrm{field}}\times d.\end{array}$$

Figure A1 Sketch of key features of MOS.

A2 MOS NPP

Figure A2 Hovemoeller plot of latitudinally and vertically integrated MOS NPP. High NPP are found in the Southern Ocean. The change of MOS NPP follows the seasonal solar radiation in UVic ESCM.

Figure A3 Temperature optimum curve of the macroalgae in MOS.

Appendix B: Impacts of MOS on global biogeochemistry

Figure B1 MOS biomass distributions. Red lines contour the maximum MOS-occupied area during the previous years. The annual macroalgal biomass of MOS in this figure is an average over a 10 year period, which includes times of low and high biomass due to the sinking of biomass. Thus, the biomass shown here is less than the biomass shown in Fig. 2.

Figure B2 Nutrients horizontal-distribution changes in the Pacific & Atlantic basins (MOS_Conti) relative to RCP 4.5 at year 2100. The nutrients trapping by MOS can be observed on the upper layers at, e.g., the Southern Ocean and mid–high latitudes in the Northern Hemisphere. The nutrients are enriched at the ocean bottom by the remineralization of MOS biomass, especially in the Southern Ocean deep waters.

Figure B3 Denitrification rate at depth 3000–6000 m in year 3000, where the oxygen level is lower than 5 µmol m⁻³, caused by the remineralization of continuously sunken MOS biomass.

Figure B4 Plot of global averaged biomass of zooplankton.

Figure B5 Cumulative (year 2020 to 3000) leakage of MOS-captured carbon in the simulation MOS_Stop.

Figure B6 Plot of global averaged phytoplankton biomass.

Figure B7 Plot of upper-ocean 0–2000 m averaged detritus remineralization rate.

Figure B8 Upper-ocean 3000–6000 m averaged dissolved-oxygen concentrations.

Figure B9 Redistributions of NO₃ avg. 0–200 m depth relative to RCP 4.5.

Figure B10 Global averaged carbon flux from atmosphere to land.

Figure B11 Global averaged carbon flux from atmosphere to ocean.

Figure B12 Global profiles of zooplankton biomass (left y axis) and macroalgal NPP of MOS (right y axis) with and without zooplankton grazing on macroalgae. The “zooplankton” communities do not have large effects, via grazing, on macroalgae NPP or its own biomass.