Assessing tipping-point risk in carbon dioxide removal with a network-based framework

INTRODUCTION

The global carbon cycle can settle into more than one steady state. This possibility, supported by a growing body of modeling evidence^[1,2], is consequential because it implies the existence of climate "tipping points": thresholds beyond which the system may shift abruptly and irreversibly into a warmer equilibrium from which recovery is practically impossible^[3,4,5]. Human-induced climate change may be pushing the Earth system closer to such thresholds^[5,6], and the mechanisms governing these transitions remain poorly understood^[7,8].

This uncertainty has direct implications for carbon dioxide removal (CDR), the family of strategies aimed at actively removing CO₂ from the atmosphere as part of long-term carbon neutrality planning^[9,10,11]. Projections of how much CO₂ each strategy can remove over a given horizon (i.e., its mitigation potential) are a core input to national carbon budgets, Intergovernmental Panel on Climate Change (IPCC) scenario pathways, and the carbon accounting that guides long-term climate policy^[9,12,13]. The existence of multiple steady states in a CDR system could have significant consequences for whether that projected potential is actually delivered. If a CDR system settles into a steady state that removes far less CO₂ than projected, the carbon accounting behind its deployment no longer holds, and the broader mitigation pathway built around it may fail to meet its emissions targets. Moreover, abrupt transitions between steady states would make this risk especially difficult to anticipate or quickly recover from, since the shift could occur faster than monitoring, reporting, and verification systems can detect. Yet a critical gap remains in our understanding of when and how CDR systems exhibit multiple steady states under varying environmental conditions^[10,14].

Current approaches to CDR risk assessment address related but different questions. Earth system models project how the climate responds to a given deployment over time, and multi-model intercomparisons have shown that the carbon cycle does not simply retrace its path when emissions turn negative^[15,16]. Integrated assessment models ask which mix of technologies is affordable under land and energy constraints, taking the removal rates themselves as given^[9,17]. Life-cycle, permanence and monitoring frameworks ask whether the carbon stored at a particular site will stay there, and whether that can be verified^[14,18,19]. Deployment-impact studies ask what scaling would cost in resources; recent work finds that removal at the gigatonne scale could strain energy, water and land systems across whole regions^[20,21]. None of them asks a prior question: whether the coupled system of technology and carbon cycle can settle into more than one steady state at all.

Earth system models would seem the natural place to look for multiple steady states, but they are not designed for the task. They are indispensable for projecting climate trajectories and human impacts^[1,22], and they achieve this by integrating the system forward in time under prescribed forcing rather than by locating the states it can settle into. Those states are seldom computed, and the slow processes that would determine them are among the least constrained by observation^[2,16]. The number of steady states a coupled CDR–carbon-cycle system admits is therefore a question posed at the level of smaller and simpler models, where the steady states are well defined and open to direct study^[3,7,8]. Even at that level, locating them by simulation requires sweeping a parameter space that is large and poorly constrained. A more tractable approach is needed, one that can characterize multistationarity, the property of admitting more than one positive steady state, from the structure of the system rather than from exhaustive numerical simulation.

Chemical Reaction Network Theory (CRNT) provides a framework well-suited to this purpose. Originally developed in the context of chemical kinetics, CRNT represents system components as interacting species and transfers between them as reactions governed by rate laws. When formulated with power-law kinetics, it admits a parameter-minimal analysis that relies primarily on the graphical and kinetic structure of the network, rather than precise parameter values, to identify sufficient conditions under which multiple steady states may arise. This makes the framework well-suited to systems in which parameter uncertainty is high. The carbon-cycle box models considered here are such systems. Reaction network modeling has been successfully applied to biochemical networks^[23,24], ecological systems^[25,26], and epidemiology^[27], and its application to the carbon cycle is an active and growing area^{[28,29,30,31]}.

In this study, we apply the Reaction Network Carbon Dioxide Removal (RNCDR) framework^[29] to analyze steady-state multiplicity in four CDR strategies: biochar sequestration, ocean fertilization, soil carbon sequestration, and wetland restoration. Each system is translated into a reaction network with power-law kinetics, in which carbon pools, namely land biota, atmosphere, ocean, total carbon stock, and CDR storage, play the role of interacting species, and carbon transfers between them are modeled as reactions. The resulting network structure is analyzed to identify parameter conditions under which the system may exhibit multistationarity or is guaranteed to be monostationary, admitting at most one positive steady state. This distinction is crucial: a multistationary CDR system could settle into a steady state characterized by low net CO₂ removal performance, undermining its reliability as a climate mitigation tool, while a monostationary system offers more predictable and recoverable dynamics. This approach provides a practical, structure-based way to screen each of the four CDR strategies for the possibility of multistable behavior. The screen comes before any commitment to detailed numerical modeling of a specific deployment.

A screen of this kind establishes less than the language of tipping points suggests, and the distinction is worth drawing at the outset. We use tipping point throughout in a specific sense: an abrupt transition from one positive steady state of the coupled system to a different one. Thus, multiple steady states are necessary since a system admitting only one has no second state to move to. They are not sufficient. A transition also requires that at least two of those states be attractors, that the one currently occupied cease to be available as conditions drift, and that a disturbance of realistic magnitude carry the system across the basin between them. Abrupt behavior of other kinds, which does not end in a second steady state, falls outside this definition and outside what we analyze; Section "Multistationarity and tipping points: what this framework does and does not establish" provides the distinction. Whether that occurs depends on perturbation magnitude, timescale and prevailing conditions, none of which is resolved here, and the other mechanisms by which a system may tip lie outside the analysis altogether.

The two kinds of results reported below therefore differ in force. A monostationarity result is an exclusion: a system admitting a unique positive steady state across a parameter range has no second state with which one could merge and vanish, and so cannot tip by that route however it is forced. A multistationarity result is a flag rather than a diagnosis: the structural precondition is present, and the configuration warrants dynamical study before the strategy is relied upon. We retain the language of tipping-point risk throughout, since it is the concern that makes the structural question worth asking, but no result in this paper establishes that any particular deployment will tip. Section "Multistationarity and tipping points: what this framework does and does not establish" develops the distinction and states its limits.

The remainder of the paper is structured as follows. Section "THE RNCDR FRAMEWORK" presents the RNCDR framework. Section "FOUR CDR STRATEGIES AND THEIR REACTION NETWORK REPRESENTATIONS" describes the four CDR systems and their reaction network representations. Section "CONDITIONS FOR MULTISTATIONARITY IN CDR REACTION NETWORKS" establishes the multistationarity and monostationarity conditions for each system. Section "COMPARATIVE ANALYSIS AND IMPLICATIONS" presents a comparative analysis of the four systems and discusses implications for CDR risk assessment. Section "SUMMARY, CONCLUSION, AND FUTURE RESEARCH" concludes with a summary of the key findings, their broader significance for Earth system science and climate policy, and directions for future research.

THE RNCDR FRAMEWORK

The RNCDR framework, formally introduced by Fortun et al.^[29], builds on and extends their earlier work^[31,32]. This approach applies CRNT to model the Earth's carbon cycle and, in the most recent studies^[28,30], to two carbon removal technologies, direct air capture (DAC) and direct ocean capture (DOC). The present study expands this line of work by applying the RNCDR framework to four additional CDR strategies, namely biochar sequestration, ocean fertilization, soil carbon sequestration, and wetland restoration, and by carrying out a comparative steady-state analysis across them to identify structural conditions of multistationarity that are common to, or distinctive of, each strategy.

A chemical reaction network (CRN) is a way of describing a system in terms of three ingredients: the species (the things that interact, in our case the carbon pools), the complexes (the combinations of species that appear on either side of an arrow, such as land plus atmosphere), and the reactions (the directional transfers from one complex to another, each describing a way carbon moves through the system). Together, these define the structure of the network. Each reaction is assigned a rate function describing how fast it proceeds, and the collection of these rates gives the dynamics of the system. In this study, the rate functions take the form of power-law kinetics, a flexible class of rate laws that can approximate a wide range of biogeochemical processes. Precise definitions are provided in Supplementary Section 1.

In this framework, the global carbon cycle is modeled as a CRN in which carbon moves among five interconnected pools: land biota ($$ A_1 $$), representing carbon stored in forests, soil, and living organisms; atmosphere ($$ A_2 $$), representing carbon as atmospheric CO₂; ocean ($$ A_3 $$), representing carbon dissolved in seawater; total carbon stock ($$ A_4 $$), which comprises the geological reserve, including fossil fuel, together with the inorganic carbon pool. It is the reservoir on which the emission reaction $$ R_5 $$ draws, and it is a species distinct from the CDR storage pools $$ A_i $$, which exchange carbon with it through the leak reactions $$ R_{i,6} $$ rather than forming part of it; and CDR storage ($$ A_i $$, $$ i = 8, \ldots, 17 $$), representing carbon captured and stored by a specific negative emissions technology (NET).

Every RNCDR system is assembled from three components: the Anderies pre-industrial subnetwork, the fossil fuel emission reaction, and the CDR carbon storage. Figure 1 provides an overview of how these three components combine to form an RNCDR system. The Anderies pre-industrial subnetwork (Component 1) tracks the natural carbon cycle among land biota, the atmosphere, and the ocean. The fossil fuel emission reaction (Component 2) introduces the total carbon stock and the transfer of carbon from this stock to the atmosphere. The CDR carbon storage (Component 3) introduces a storage pool $$ A_i $$ associated with a specific NET, together with at least one reaction returning carbon to the total carbon stock and one or more technology-specific reactions. Together, these three components are assembled into a single chemical reaction network with power-law kinetics, with five species ($$ A_1 $$, $$ A_2 $$, $$ A_3 $$, $$ A_4 $$, and $$ A_i $$) and at least seven reactions. The remainder of this section describes each component in detail.

Figure 1. Structure of an RNCDR system. Every RNCDR model is assembled from three components: the Anderies pre-industrial subnetwork representing the natural carbon cycle among land biota ($$ A_1 $$), atmosphere ($$ A_2 $$), and ocean ($$ A_3 $$); the fossil fuel emission reaction connecting the total carbon stock ($$ A_4 $$) to the atmosphere; and the CDR carbon storage pool ($$ A_i $$) representing a specific NET. The three components combine into a single chemical reaction network with power-law kinetics.

Classifying RNCDR systems

The dynamic behavior of an RNCDR system is strongly shaped by the properties of its Anderies subnetwork. Following^[31], we classify Anderies systems – and by extension, RNCDR systems^[29]– according to the ratio

(9) \begin{equation} R = \frac{p_2 - p_1}{q_2 - q_1}, \quad q_2 - q_1 \neq 0, \end{equation}

which compares how much respiration and photosynthesis differ in their influence on land biota versus atmosphere. The kinetic orders $$ p_1, p_2 $$ encode the degree to which photosynthesis and respiration are influenced by the size of the land biota pool, while $$ q_1, q_2 $$ encode the degree to which they are influenced by the level of atmospheric CO₂. A value of $$ p_2 > p_1 $$ means that respiration responds more strongly to land biota than photosynthesis does; in other words, as forests grow, they release carbon faster than they absorb it. Similarly, $$ q_2 > q_1 $$ means respiration is more strongly driven by atmospheric CO₂ concentration than photosynthesis is. When both inequalities hold simultaneously, these asymmetries reinforce each other, creating the nonlinear feedback conditions under which multiple equilibria can arise. Four primary classes arise from the ratio R, each with distinct dynamic consequences:

● Positive class ($$ R > 0 $$): both $$ p_2 > p_1 $$ and $$ q_2 > q_1 $$ (or both reversed), so the same process dominates the feedback in both pools: respiration when $$ p_2 > p_1 $$ and $$ q_2 > q_1 $$, photosynthesis when both inequalities are reversed. What makes R positive is that the two inequalities point in the same direction, not the identity of the dominant process; the Anderies subnetwork may admit multiple steady states, meaning the Earth could settle into more than one equilibrium carbon distribution.

● Negative class ($$ R < 0 $$): photosynthesis dominates one feedback while respiration dominates the other, leading to more stabilizing dynamics; the Anderies subnetwork contains monostationary subsets (unique steady states under appropriate kinetic order conditions).

● P-null class ($$ R = 0 $$): the Anderies subnetwork is monostationary.

● Q-null class ($$ q_2 - q_1 = 0,\ p_2 - p_1 \neq 0 $$): defined for the case where the ratio R is undefined but $$ Q = (q_2 - q_1)/(p_2 - p_1) = 0 $$. This class shares the monostationarity of the P-null class in the Anderies subnetwork.

These classes characterize the dynamic behavior of the Anderies subnetwork in isolation. Once the Anderies subnetwork is embedded in an RNCDR system, the addition of the fossil fuel emission reaction and the CDR-specific reactions can change the steady-state behavior: a class that is monostationary at the Anderies level need not remain so for the full RNCDR system, and the propositions in Section "CONDITIONS FOR MULTISTATIONARITY IN CDR REACTION NETWORKS" make this explicit for each of the four CDR strategies analyzed.

Note the asymmetry in these definitions: the P-null class arises when the numerator of R vanishes (i.e., $$ p_2 = p_1 $$), while the Q-null class arises when the denominator vanishes (i.e., $$ q_2 = q_1 $$) with a nonzero numerator. The two classes are therefore structurally distinct, and as Propositions 3 and 4 demonstrate, they can exhibit different multistationarity behavior within the same CDR system.

Therefore, an RNCDR kinetic system for a single CDR technology consists of five species ($$ A_1 $$, $$ A_2 $$, $$ A_3 $$, $$ A_4 $$, and one CDR storage species $$ A_i $$) and at least seven reactions. The latter include the five reactions $$ R_1 $$–$$ R_5 $$ that are common to all RNCDR models, the common CDR storage reaction $$ R_{i,6}: A_i \to A_4 $$, and at least one CDR-specific reaction $$ R_{i,j} $$ for $$ j \geq 7 $$, whose number and form depend on the particular CDR technology deployed. The dynamics of the system follow power-law kinetics encoded by the matrix F in (8), with the emission kinetic orders $$ e_i $$ and $$ f_i $$ determined by the storage parameters $$ \lambda_i $$ and $$ \mu_i $$ via (5).

Figure 2 summarizes how the RNCDR framework is applied in practice. The analysis proceeds in five steps: (1) selecting a CDR strategy; (2) parametrizing its storage through $$ \lambda_i $$ and $$ \mu_i $$; (3) representing the system as a chemical reaction network together with its kinetic order matrix F; (4) classifying the underlying Anderies subnetwork through the ratio $$ R = (p_2 - p_1)/(q_2 - q_1) $$; and (5) testing for multistationarity or monostationarity through the Deficiency One Algorithm, the Multistationarity Algorithm, and the injectivity test (discussed in Section "CONDITIONS FOR MULTISTATIONARITY IN CDR REACTION NETWORKS"). The outcome of this workflow is a determination of whether a given CDR system can admit multiple steady states, and therefore whether it may be vulnerable to tipping-point behavior under varying conditions. Notably, this characterization is obtained from the graphical and kinetic structure of the network alone, without requiring precise knowledge of rate constants. The four CDR strategies analyzed in this study, presented in Section "FOUR CDR STRATEGIES AND THEIR REACTION NETWORK REPRESENTATIONS", all follow this same workflow, with the differences between them arising from the specific reactions and kinetic orders introduced at steps 2 and 3.

Figure 2. The RNCDR analysis workflow. Each CDR strategy is screened for tipping-point risk through five sequential steps: selecting the strategy, parametrizing its storage, representing it as a chemical reaction network, classifying the Anderies subnetwork, and testing for multistationarity or monostationarity. The framework relies on the graphical and kinetic structure of the network rather than precise rate constants.

The key question the framework is designed to answer is: Under what conditions can a CDR system settle into more than one steady state, and therefore become vulnerable to a tipping point that would compromise the carbon removal it was designed to deliver? As explored in the subsequent sections, the answer depends critically on the class (positive, negative, or null) of the underlying Anderies subnetwork, and on the structural features of the CDR-specific reactions, including their kinetic orders and the resulting network deficiency. This provides a basis for screening CDR strategies for carbon-accounting risk before their projected mitigation potential is committed to in deployment plans or long-term climate policy.

FOUR CDR STRATEGIES AND THEIR REACTION NETWORK REPRESENTATIONS

This section presents the reaction network representations of the four CDR strategies analyzed in this study --- biochar sequestration, ocean fertilization, soil carbon sequestration, and wetland restoration --- together with the box models used to represent each system, their corresponding RNCDR reaction network translation, and the kinetic order matrices encoding their dynamics. Each system extends the RNCDR framework of Section "THE RNCDR FRAMEWORK" by specifying the CDR storage species $$ A_i $$, the storage parameters $$ \lambda_i $$ and $$ \mu_i $$, and the CDR-specific reactions $$ R_{i,6} $$ and $$ R_{i,j} $$ for $$ j \geq 7 $$. Casting all four strategies in a common reaction network form makes their structural differences explicit, and provides the basis for the comparative analysis of multistationarity in Sections "CONDITIONS FOR MULTISTATIONARITY IN CDR REACTION NETWORKS" and "COMPARATIVE ANALYSIS AND IMPLICATIONS". This is one way in which a reaction network approach can complement existing assessments of CDR options: by exposing structural features that may be relevant to choosing among strategies, but that may not be apparent from numerical simulation alone. The full ODE systems for each model are provided in Supplementary Section 3 for reference.

The four strategies were chosen to span the ways in which CDR systems physically hold carbon and the settings in which they are deployed, rather than to sample a single class of technology. Biochar sequestration stores carbon as a solid, thermochemically stabilized product applied to agricultural soil. Ocean fertilization stores it in the deep ocean by way of the marine biological pump. Soil carbon sequestration stores it as soil organic matter maintained by land management practice. Wetland restoration stores it in waterlogged sediment, where anoxia suppresses decomposition^[43]. The four therefore differ in the medium that holds the carbon, in how long it is held, and in how it is lost: biochar by slow oxidation of the applied char^[40,49], ocean carbon by overturning circulation returning it to the surface, soil carbon by tillage or a reversal of land use, and wetland carbon by drainage^[19].

They differ likewise in deployment setting and maturity. Biochar is applied in a distributed way on agricultural land and is already produced commercially. Soil carbon sequestration operates over large areas of existing farmland at low cost per hectare, but its gains saturate as soils approach a new equilibrium and are difficult to verify^[42,50]. Wetland restoration is undertaken at the scale of individual coastal and inland sites, usually alongside biodiversity and flood-protection objectives. Ocean fertilization requires no land at all, but its efficacy is contested and it is restricted under international marine agreements. The set thus ranges from a commercially delivered engineered product, through two forms of land and ecosystem management, to an open-ocean intervention that remains experimental.

This spread is what makes the four informative within the RNCDR framework, since permanence and reversibility enter the model directly through $$ \lambda_i $$, $$ \mu_i $$ and the leak reaction $$ R_{i,6} $$. As a consequence the four occupy near-extremal positions in the $$ (\lambda_i,\mu_i) $$ plane, namely $$ (1,0.5) $$, $$ (1,0) $$, $$ (0.01,0) $$ and $$ (0.01,0.9) $$; they cover both deficiency classes, with biochar and soil carbon sequestration admitting the Deficiency One Algorithm and ocean fertilization and wetland restoration requiring the Multistationarity Algorithm; and they cover the range of coupling to the natural carbon cycle, from CDR reactions that leave $$ R_1 $$–$$ R_4 $$ structurally unchanged to explicit bidirectional atmosphere–storage exchange.

The remaining technologies in Table 2 are excluded for specific reasons. Direct air capture, direct ocean capture and bioenergy with carbon capture and storage have already been analyzed within this framework^[28,29,30], and all three store captured CO₂ by injection into geological formations, so including them here would duplicate both the published analysis and a storage mechanism already covered. Enhanced weathering and ocean alkalinization store carbon geochemically, as dissolved bicarbonate and carbonate rather than as organic matter, and the carbon follows a long multi-step path from mineral grain through soil solution and rivers to the ocean; representing that path faithfully would require pools beyond the five used here, and we therefore treat the two together as the natural next case. Afforestation and reforestation store carbon in standing biomass and soil on a route close to that of biochar, and would add no new coverage of the storage-parameter space, since $$ \lambda_{15} = 0.5 $$ and $$ \mu_{15} = 0 $$ give $$ c_{15} = 0.5 $$, the same withheld fraction as biochar; they are the subject of a separate study in preparation by a group that includes several of the present authors.

Biochar sequestration

Biochar sequestration, hereafter referred to as the BCS system, involves incorporating carbon-rich biochar into soils to enhance long-term carbon storage^[40]. Biochar is produced from biomass through pyrolysis, a process of burning organic material in an oxygen-limited environment, creating a stable, carbon-rich material that resists decomposition. In the RNCDR framework, biochar storage is represented by $$ A_{11} $$, with parameters $$ \lambda_{11} = 1 $$ and $$ \mu_{11} = 0.5 $$: all captured carbon enters long-term storage, half as inorganic biocharcoal in soil and half as organic biofuel in the geological stock.

The box model is shown in Figure 3. The BCS system involves five carbon pools: land biota ($$ A_1 $$), atmosphere ($$ A_2 $$), ocean ($$ A_3 $$), total carbon stock ($$ A_4 $$), and biochar storage ($$ A_{11} $$). Biomass is converted into biochar through pyrolysis ($$ A_1 \to A_{11} $$). When applied to soil, biochar stimulates plant growth and photosynthesis, facilitating further carbon capture ($$ A_{11} \to A_1 $$). Biochar may also undergo fossilization into the total carbon stock ($$ A_{11} \to A_4 $$), from which carbon can be re-emitted to the atmosphere through fossil fuel combustion ($$ A_4 \to A_2 $$). The CDR-specific reactions are:

\begin{align*} R_{11,6} &: A_{11} \to A_4 && \text{(biochar to total carbon stock)}, \\ R_{11,7} &: A_1 \to A_{11} && \text{(land-to-biochar capture)}, \\ R_{11,8} &: A_{11} \to A_1 && \text{(biochar-to-land transfer)}. \end{align*}

Figure 3. Box model of the biochar sequestration system. The green-outlined pools and reactions ($$ A_1 $$, $$ A_2 $$, $$ A_3 $$, and the transfers among them) reproduce the pre-industrial Anderies subnetwork; the black-outlined pools and reactions are specific to biochar and to the fossil-fuel accounting that links the total carbon stock to the atmosphere. The same convention is used in the box models of the remaining CDR systems.

Together with the four Anderies reactions $$ R_1 $$–$$ R_4 $$ and the emission reaction $$ R_5 $$, this gives a BCS system with 5 species and 8 reactions. The dynamics of this system are encoded in the kinetic order matrix:

\begin{equation*} F_{BCS} = \begin{array}{c} A_1 \;\;\; A_2 \;\;\; A_3 \;\;\; A_4 \;\;\; A_{11} \\ \left[\begin{array}{[ccccc]} p_1 & q_1 & 0 & 0 & 0 \\ p_2 & q_2 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & e_{11} & f_{11} \\ 0 & 0 & 0 & 0 & 1 \\ 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ \end{array} \right] \end{array} \quad \begin{array}{l} \\ R_1 \\ R_2 \\ R_3 \\ R_4 \\ R_5 \\ R_{11,6} \\ R_{11,7} \\ R_{11,8} \end{array} \end{equation*}

The row for $$ R_5 $$ encodes the coupled dependence of emission on both the total carbon stock and biochar storage, with $$ e_{11} \geq 1 $$, $$ f_{11} \leq 0 $$, and $$ e_{11} + f_{11} = 1 $$. The rows for $$ R_{11,6} $$, $$ R_{11,7} $$, and $$ R_{11,8} $$ reflect mass-action kinetics for the biochar-specific transfers.

Ocean fertilization

Ocean fertilization, hereafter referred to as the OF system, seeks to stimulate the growth of marine phytoplankton to enhance the oceanic uptake of CO₂^[41]. The process introduces nutrients such as iron to the ocean surface, inducing phytoplankton blooms that absorb atmospheric carbon into their biomass. When the phytoplankton die, their biomass sinks to the deep ocean where the carbon is stored for a long period of time. The efficacy of ocean fertilization is contested. Field experiments have generally found that only a small and highly variable fraction of the carbon fixed in a bloom reaches the deep ocean, and estimates of the cost per tonne removed span roughly two orders of magnitude, driven largely by uncertainty in that export efficiency^[41,51]. Deployment beyond legitimate scientific research has been restricted since 2008 by resolution of the parties to the London Convention and Protocol, an approach codified in a 2013 amendment to the London Protocol that has not yet entered into force^[41]. We retain it here as a modeling case rather than as an endorsement: it continues to appear in CDR portfolio inventories, and it is the analytically distinctive member of our set, being the only one of the four whose emission reaction reduces exactly to linear mass action. In the RNCDR framework, OF storage is represented by the separate pool $$ A_{12} $$, which holds carbon in the deep ocean and is distinct from the total carbon stock $$ A_4 $$. Carbon leaves $$ A_{12} $$ for $$ A_4 $$ only through the leak reaction $$ R_{12,6} $$, which represents fossilization on geological timescales. The parameters $$ \lambda_{12} = 1 $$ and $$ \mu_{12} = 0 $$ describe the disposition of that transferred carbon, namely that it is organic and therefore available for subsequent emission, rather than the physical location of the OF store itself. Since $$ \lambda_{12} = 1 $$ and $$ \mu_{12} = 0 $$, the emission kinetics of the OF system reduce to linear mass action: $$ K_5 = k_5 A_4 $$.

The box model is shown in Figure 4. The OF system involves five carbon pools: land biota ($$ A_1 $$), atmosphere ($$ A_2 $$), ocean ($$ A_3 $$), total carbon stock ($$ A_4 $$), and OF storage ($$ A_{12} $$). The OF process introduces an enhanced photosynthesis pathway from the atmosphere to the ocean, a transfer from the ocean to OF storage representing biomass sinking to the deep ocean ($$ A_3 \to A_{12} $$), and a reverse transfer from OF storage back to the ocean representing ocean water mixing ($$ A_{12} \to A_3 $$). The CDR-specific reactions are:

\begin{align*} R_{12,6} &: A_{12} \to A_4 && \text{(OF to total carbon stock)}, \\ R_{12,7} &: 2A_2 + A_3 \to A_2 + 2A_3 && \text{(atmosphere-to-ocean photosynthesis)}, \\ R_{12,8} &: 2A_3 \to A_3 + A_{12} && \text{(ocean-to-OF capture)}, \\ R_{12,9} &: A_3 + A_{12} \to 2A_3 && \text{(OF-to-ocean release)}. \end{align*}

Figure 4. Box model of the OF system. As in Figure 3, the green-outlined portion reproduces the Anderies subnetwork.

Together with $$ R_1 $$–$$ R_4 $$ and $$ R_5 $$, this gives the OF system with 5 species and 9 reactions. The dynamics of the OF system are encoded in the kinetic order matrix:

\begin{equation*} F_{OF} = \begin{array}{c} A_1 \;\;\; A_2 \;\;\; A_3 \;\;\; A_4 \;\;\; A_{12} \\ \left[\begin{array}{[ccccc]} p_1 & q_1 & 0 & 0 & 0 \\ p_2 & q_2 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ 0 & p_3 & q_3 & 0 & 0 \\ 0 & 0 & p_4 & 0 & 0 \\ 0 & 0 & p_5 & 0 & q_5 \\ \end{array}\right] \end{array} \quad \begin{array}{l} \\ R_1 \\ R_2 \\ R_3 \\ R_4 \\ R_5 \\ R_{12,6} \\ R_{12,7} \\ R_{12,8} \\ R_{12,9} \end{array} \end{equation*}

The row for $$ R_5 $$ reflects linear mass action kinetics ($$ e_{12} = 1 $$, $$ f_{12} = 0 $$) since $$ \lambda_{12} = 1 $$ and $$ \mu_{12} = 0 $$. The rows for $$ R_{12,7} $$, $$R_{12,8} $$, and $$ R_{12,9} $$ introduce additional kinetic orders $$ p_3, q_3, p_4, p_5 $$, and $$ q_5 $$ that encode the influence of the atmosphere and ocean pools on the OF processes.

Soil carbon sequestration

Soil carbon sequestration, hereafter referred to as the SCS system, aims to increase the organic carbon content of agricultural and natural soils through various land management practices^[42]. Carbon from the atmosphere is absorbed by plants through photosynthesis and stored in the soil as soil organic carbon. Practices such as cover cropping, crop rotation, and reduced tillage are used to increase the rate at which carbon is added to the soil while decreasing the rate at which it is lost through decomposition. In the RNCDR framework, SCS storage is represented by $$ A_{13} $$, with parameters $$ \lambda_{13} = 0.01 $$ and $$ \mu_{13} = 0 $$: only a very small fraction of sequestered carbon enters long-term geological storage, and this fraction is entirely organic.

The box model is shown in Figure 5. The SCS system involves five carbon pools: land biota ($$ A_1 $$), atmosphere ($$ A_2 $$), ocean ($$ A_3 $$), total carbon stock ($$ A_4 $$), and SCS storage ($$ A_{13} $$). Carbon decays from land biota into SCS storage ($$ A_1 \to A_{13} $$), from which it may enter the total carbon stock through fossilization ($$ A_{13} \to A_4 $$). From there, carbon can be re-emitted to the atmosphere through fossil fuel combustion ($$A_4 \to A_2 $$). The CDR-specific reactions are:

\begin{align*} R_{13,6} &: A_{13} \to A_4 && \text{(SCS to total carbon stock)}, \\ R_{13,7} &: A_1 \to A_{13} && \text{(land-to-SCS capture)}. \end{align*}

Figure 5. Box model of the SCS system. As in Figure 3, the green-outlined portion reproduces the Anderies subnetwork.

Together with $$ R_1 $$–$$ R_4 $$ and $$ R_5 $$, this gives the SCS system with 5 species and 7 reactions. The dynamics of the SCS system are encoded in the kinetic order matrix:

\begin{equation*} F_{SCS} = \begin{array}{c} A_1 \;\;\; A_2 \;\;\; A_3 \;\;\; A_4 \;\;\; A_{13} \\ \left[\begin{array}{[ccccc]} p_1 & q_1 & 0 & 0 & 0 \\ p_2 & q_2 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & e_{13} & f_{13} \\ 0 & 0 & 0 & 0 & 1 \\ 1 & 0 & 0 & 0 & 0 \\ \end{array}\right] \end{array} \quad \begin{array}{l} \\ R_1 \\ R_2 \\ R_3 \\ R_4 \\ R_5 \\ R_{13,6} \\ R_{13,7} \end{array} \end{equation*}

The row for $$ R_5 $$ encodes the coupled dependence of emission on both the total carbon stock and SCS storage, with $$ e_{13} \geq 1 $$, $$ f_{13} \leq 0 $$, and $$ e_{13} + f_{13} = 1 $$. The rows for $$ R_{13,6} $$ and $$ R_{13,7} $$ reflect mass-action kinetics for the SCS-specific transfers. As with the BCS system, the CDR storage species $$ A_{13} $$ does not appear in the rows corresponding to the natural carbon cycle reactions $$ R_1 $$–$$ R_4 $$, so the SCS-specific reactions extend the natural carbon cycle without directly modulating its rates.

Wetland restoration

Wetland restoration, hereafter referred to as the WR system, leverages the natural carbon sequestration capacity of wetland ecosystems^[43]. Wetlands are ecosystems where plants and soil are submerged in water; they can capture large amounts of CO₂ from the atmosphere and store it in soils for hundreds or even thousands of years. In the RNCDR framework, WR storage is represented by $$ A_{14} $$, with parameters $$ \lambda_{14} = 0.01 $$ and $$ \mu_{14} = 0.9 $$: only a very small fraction of sequestered carbon enters long-term geological storage, and most of this fraction is inorganic.

The box model is shown in Figure 6. The WR system involves five carbon pools: land biota ($$ A_1 $$), atmosphere ($$ A_2 $$), ocean ($$ A_3 $$), total carbon stock ($$ A_4 $$), and WR storage ($$ A_{14} $$). The WR system is the most complex of the four, with carbon transfers between the atmosphere and WR storage in both directions, mediated by respiration and photosynthesis processes. Carbon also transfers between land biota and WR storage in both directions ($$ A_1 \to A_{14} $$ and $$ A_{14} \to A_1 $$), reflecting the interaction between land biota and wetland ecosystems. The CDR-specific reactions are:

\begin{align*} R_{14,6} &: A_{14} \to A_4 && \text{(wetland to total carbon stock)}, \\ R_{14,7} &: A_2 + 2A_{14} \to 2A_2 + A_{14} && \text{(wetland respiration)}, \\ R_{14,8} &: 2A_2 + A_{14} \to A_2 + 2A_{14} && \text{(wetland photosynthesis)}, \\ R_{14,9} &: A_1 \to A_{14} && \text{(land-to-wetland capture)}, \\ R_{14,10} &: A_{14} \to A_1 && \text{(wetland-to-land transfer)}. \end{align*}

Figure 6. Box model of the WR system. As in Figure 3, the green-outlined portion reproduces the Anderies subnetwork.

Together with $$ R_1 $$–$$ R_4 $$ and $$ R_5 $$, this gives the WR system with 5 species and 10 reactions. The dynamics of the WR system are encoded in the kinetic order matrix:

\begin{equation*} F_{WR} = \begin{array}{c} A_1 \;\;\; A_2 \;\;\; A_3 \;\;\; A_4 \;\;\; A_{14} \\ \left[\begin{array}{[ccccc]} p_1 & q_1 & 0 & 0 & 0 \\ p_2 & q_2 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & e_{14} & f_{14} \\ 0 & 0 & 0 & 0 & 1 \\ 0 & q_4 & 0 & 0 & p_4 \\ 0 & q_5 & 0 & 0 & p_5 \\ 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ \end{array}\right] \end{array} \quad \begin{array}{l} \\ R_1 \\ R_2 \\ R_3 \\ R_4 \\ R_5 \\ R_{14,6} \\ R_{14,7} \\ R_{14,8} \\ R_{14,9} \\ R_{14,10} \end{array} \end{equation*}

The WR system has the largest kinetic order matrix among the four systems, reflecting its greater complexity. The rows for $$ R_{14,7} $$ and $$ R_{14,8} $$ introduce additional kinetic orders $$ q_4, p_4, q_5 $$, and $$ p_5 $$ that encode the influence of both the atmosphere ($$ A_2 $$) and WR storage ($$ A_{14} $$) on the respiration and photosynthesis processes within the wetland ecosystem. This bidirectional coupling between $$ A_2 $$ and $$ A_{14} $$ is a structural feature unique to the WR system.

CONDITIONS FOR MULTISTATIONARITY IN CDR REACTION NETWORKS

Having established the RNCDR representations of the four CDR systems, we now present the conditions under which each system may exhibit multistationarity, the capacity to admit more than one positive steady state, or is guaranteed to be monostationary, settling into at most one. Recall from Section "THE RNCDR FRAMEWORK" that the dynamic behavior of the underlying Anderies subnetwork is organized by the ratio $$ R = (p_2 - p_1)/(q_2 - q_1) $$. It was shown in^[31] that a positive class ($$ R > 0 $$) is associated with potential multistationarity at the level of the Anderies subnetwork, while negative, P-null, and Q-null classes are associated with monostationarity under appropriate conditions. When the Anderies subnetwork is embedded in a full RNCDR system, however, the addition of CDR-specific reactions can either preserve or change this baseline behavior. The propositions below establish, for each of the four CDR systems, the precise conditions under which the full RNCDR system inherits, departs from, or qualitatively rewrites the Anderies-level behavior.

The analysis proceeds by applying two complementary tools to each CRN representation. A key structural index that determines which tools are applicable is the deficiency of the network. The deficiency is a non-negative measure that quantifies the degree of "linear independence" among a network's reactions^[52]. A deficiency of zero represents the highest level of linear independence consistent with the network's directed graph structure. Conversely, a higher deficiency indicates a lower degree of linear independence, indicating greater complexity in the network structure. As shown in Supplementary Section 3, the BCS and SCS systems have deficiency $$ 1 $$, while the OF and WR systems have deficiency $$ 2 $$.

To establish monostationarity, we apply the injectivity test of Wiuf and Feliu^[53,54], which uses a computational approach and Maple script to determine whether the determinant of the matrix $$ M^* $$ has all positive or all negative coefficients; if so, the network is injective and hence monostationary. To establish multistationarity, we apply either the Deficiency One Algorithm (DOA) for power-law kinetic systems ^[32] for the BCS and SCS systems (deficiency 1), or the Multistationarity Algorithm (MSA) of Hernandez et al.^[55] for the OF and WR systems (deficiency 2). Both the DOA and MSA translate the question of whether multiple positive steady states can exist, a nonlinear problem, into a tractable linear system of equations and inequalities, whose solvability is equivalent to the existence of multiple steady states. The detailed computations and proofs are provided in the Supplementary Materials; here we state and interpret the results.

Proposition 1 (Biochar Sequestration). The following hold for the BCS system:

(i) Positive class. If $$ p_2 > p_1 $$ and $$ q_2 > q_1 $$, or $$ p_1 > p_2 $$ and $$ q_1 > q_2 $$, then the system is multistationary.

(ii) Negative class. If $$ p_1 < 0 $$, $$ p_2 > 0 $$, $$ q_1 > 0 $$, $$ q_2 < 0 $$, $$ e_{11} \geq 1 $$, and $$ f_{11} \leq 0 $$, then the system is monostationary.

(iii) P-null class. The system is multistationary.

(iv) Q-null class. The system is multistationary.

The conditions $$ p_2 > p_1 $$ and $$ q_2 > q_1 $$ mean that respiration is more strongly influenced by both land biota and atmosphere than photosynthesis is. When these inequalities hold, the nonlinear feedback between the carbon pools is strong enough to sustain multiple equilibria. Physically, this corresponds to a scenario in which the Earth's natural carbon cycle is itself in a potentially multistationary regime, and BCS does not eliminate this underlying instability. Beyond the positive class, multistationary subsets have also been identified in the P-null and Q-null classes, indicating that BCS can admit multiple steady states even when the Anderies subnetwork is in a structurally degenerate regime (i.e., when the kinetic-order ratio R is zero or undefined). The negative class remains the one Anderies regime in which the BCS system is guaranteed to be monostationary, and therefore the one in which the carbon dynamics admit at most one positive equilibrium distribution, a condition for predictable and recoverable long-term behavior.

Proposition 2 (Ocean Fertilization). The following hold for the OF system:

(i) Positive class. If $$ q_1 > q_2 > 0 $$ and $$ p_1 > p_2 > 0 $$ (or $$ q_2 > q_1 > 0 $$ and $$ p_2 > p_1 > 0 $$), together with $$ p_3, q_3, p_4, p_5, q_5 > 0 $$ and $$ p_3 \mu_{A_2} > M_3 $$, then the system is multistationary.

(ii) Negative class. If $$ q_2 > q_1 > 0 $$ and $$ p_1 > 0 > p_2 $$, together with $$ p_3, q_3, p_5, q_5 > 0 $$, $$ p_4 < 0 $$ and $$ p_3 \mu_{A_2} > M_3 $$, then the system is monostationary.

(iii) Q-null class. If $$ q_2 = q_1 $$ and $$ p_1 > 0 > p_2 $$, together with $$ p_3, q_3, p_5, q_5 > 0 $$, $$ p_4 < 0 $$ and $$ p_3 \mu_{A_2} > M_3 $$, then the system is monostationary.

(iv) P-null class. If $$ p_2 = p_1 $$, $$ q_2 \neq q_1 $$, and $$ q_3 > 0 > p_3 $$, then the system is monostationary.

A notable feature of the OF system is that, since $$ \lambda_{12} = 1 $$ and $$ \mu_{12} = 0 $$, the emission reaction reduces to linear mass action ($$ K_5 = k_5 A_4 $$), the simplest possible form. This is the one system of the four for which the withheld fraction $$ c_i $$ vanishes, $$ c_{12} = 0 $$, and the reduction is contingent on that: any departure from these values, meaning $$ \lambda_{12} < 1 $$ or $$ \mu_{12} > 0 $$, would give $$ f_{12} < 0 $$ and move the system off mass action. This means the total carbon stock contributes to atmospheric CO₂ at a rate directly proportional to its size, with no modulating effect from the OF storage pool. Despite this simplification, the OF system retains the same potential for multistationarity as the more complex systems. The additional kinetic orders $$ p_3, q_3, p_4, p_5 $$, and $$ q_5 $$ introduced by the OF-specific reactions must also satisfy appropriate sign conditions for multistationarity to arise, reflecting the richer kinetic structure of the OF system relative to the other three. In addition to these sign conditions, the positive, negative, and Q-null cases each carry the requirement $$ p_3 \mu_{A_2} > M_3 $$, where $$ \mu_{A_2} $$ is the equilibrium value of the atmospheric pool $$ A_2 $$ and $$ M_3 $$ is a parameter arising from the multistationarity analysis (see Supplementary Materials). Intuitively, this requirement constrains how strongly the OF-mediated atmosphere-to-ocean photosynthesis pathway depends on the atmospheric pool relative to the rest of the network; it must be satisfied for the parametric conditions identified in Proposition 2 to apply. In real-world terms, this means that the steady-state behavior of an OF deployment depends not only on which Anderies class the natural carbon cycle is in, but also on whether the enhanced photosynthesis pathway is sensitive enough to atmospheric CO₂ for the framework's sufficient conditions to be in effect.

Proposition 3 (Soil Carbon Sequestration). The following hold for the SCS system:

(i) Positive class. If $$ p_2>p_1>0 $$ and $$ q_2>q_1>0 $$, or $$ p_1>p_2>0 $$ and $$ q_1>q_2>0 $$, then the system is multistationary.

(ii) Negative class. If $$ p_1,q_2<0 $$ and $$ p_2,q_1>0 $$, together with $$ e_{13}\ge1 $$ and $$ f_{13}\le0 $$, then the system is monostationary.

(iii) Q-null class. If $$ q_2 = q_1 $$ and $$ p_1>p_2>0 $$, then the system is multistationary.

(iv) P-null class. If $$ p_2 = p_1 $$ and $$ q_1\neq q_2 $$ (with $$ q_1\neq0 $$), then the system is multistationary.

Despite differences in their underlying physical processes, the SCS and BCS systems exhibit the same multistationarity profile across Anderies classes: both admit multistationary subsets in the positive, P-null, and Q-null classes, with monostationarity guaranteed only in the negative class under appropriate sign conditions. The two systems therefore share the broadest reach of multistationarity across Anderies classes among the four strategies analyzed.

Both BCS and SCS admit the DOA, which identifies parameter conditions under which multiple steady states exist in the positive, P-null, and Q-null classes for these networks. The DOA's applicability follows from the deficiency-one structure of both networks together with the specific kinetic-order pattern in which the CDR storage species ($$ A_{13} $$ for SCS, $$ A_{11} $$ for BCS) appears only in rows corresponding to CDR-specific reactions in the kinetic order matrix, leaving the natural carbon cycle reactions $$ R_1 $$–$$ R_4 $$ structurally unchanged. Under this combined structural setting, multistationarity can arise in the P-null and Q-null classes, classes that are monostationary at the level of the Anderies subnetwork in isolation. From a policy standpoint, this means that SCS deployments may face tipping-point risks not only when the natural carbon cycle is in a regime that would produce them on its own, but also under regimes that would otherwise be classified as structurally degenerate (i.e., when the kinetic-order ratio R is zero or undefined). Analyses that classify CDR risk solely on the basis of the Anderies regime of the natural carbon cycle therefore risk missing the additional multistationarity pathways that deficiency-one CDR systems of the kind shared by BCS and SCS can introduce.

Proposition 4 (Wetland Restoration) The following hold for the WR system:

(i) Positive class. If $$ p_1 \neq p_2 $$ or $$ p_4 \neq p_5 $$, together with $$ e_{14} > 0 $$ and $$ f_{14} < 0 $$, then the system is monostationary.

(ii) Negative class. If $$ p_1 \neq p_2 $$ or $$ p_4 \neq p_5 $$, then the system is multistationary.

(iii) Q-null class.

(a) If $$ (q_4 - q_5)/ (p_5 - p_4) > 0 $$ and $$ p_4 \neq p_5 $$, then the system is monostationary.

(b) If $$ (q_4 - q_5)/ (p_5 - p_4) < 0 $$ and $$ p_4 \neq p_5 $$, then the system is multistationary.

The WR system exhibits the most nuanced multistationarity behavior among the four systems. Unlike the BCS, OF, and SCS systems, where multistationarity arises in the positive class, the WR system is monostationary in the positive class and multistationary in the negative class. This is a consequence of the richer network structure of the WR system, which involves bidirectional transfers between the atmosphere and WR storage mediated by additional kinetic orders $$ q_4, p_4, q_5 $$, and $$ p_5 $$. These additional parameters introduce a new ratio $$ (q_4 - q_5)/(p_5 - p_4) $$ that governs the behavior of the Q-null class, independent of the Anderies subnetwork class. This makes the WR system the most structurally complex of the four systems, with multistationarity behavior that depends on the largest number of kinetic parameters among the four systems analyzed.

This result also depends on the storage parameterization in a way the others do not. Proposition 4(i) requires $$ f_{14} < 0 $$ strictly, which holds only when $$ c_{14} \neq 0 $$. Table 3 assigns wetland restoration $$ \lambda_{14} = 0.01 $$ and $$ \mu_{14} = 0.9 $$, giving $$ c_{14} = 0.999 $$, so the condition is comfortably met; under a parameterization with $$ \lambda_{14} = 1 $$ and $$ \mu_{14} = 0 $$ one would have $$ f_{14} = 0 $$ and Proposition 4(i) would not apply. Monostationarity of the wetland system in the positive Anderies class is therefore contingent on those values, and Supplementary Section 2 sets out the dependence in full.

This reversal of the positive- and negative-class behavior carries an implication for how WR might be represented in larger climate models. Models that capture the natural carbon cycle but treat the wetland–atmosphere interaction as a one-way flux to a storage pool would correspond, in the RNCDR framework, to a system without the bidirectional coupling encoded in $$ q_4, p_4, q_5, $$ and $$ p_5 $$, and therefore to a deficiency-one system without atmosphere–storage interaction, structurally closer to the BCS and SCS cases. Under such a simplified representation, the wetland system would be expected to inherit the multistationarity pattern that BCS and SCS exhibit, in particular multistationarity in the positive Anderies class. The reversal observed here, with monostationarity in the positive class once the two-way exchange between the atmosphere and the wetland is explicitly represented, indicates that this inherited behavior may not be the right baseline expectation for WR. Whether the same reversal appears in higher-dimensional climate models is a question this framework cannot settle on its own, since the kinetic orders $$ q_4, p_4, q_5, $$ and $$ p_5 $$ would need to be estimated from observational or modeling data before the ratio $$ (q_4 - q_5)/(p_5 - p_4) $$ could be evaluated for a specific deployment context. What the framework does offer is a structural reason to expect that the level of detail at which wetland–atmosphere exchange is represented can change the qualitative steady-state behavior of the coupled system, not just the quantitative size of the carbon flux.

COMPARATIVE ANALYSIS AND IMPLICATIONS

SUMMARY, CONCLUSION, AND FUTURE RESEARCH

This study applied the RNCDR framework to analyze the steady-state multiplicity of four CDR strategies, namely biochar sequestration, ocean fertilization, soil carbon sequestration, and wetland restoration, by translating each system into a chemical reaction network with power-law kinetics and examining the graphical and kinetic structure of the resulting networks. The analysis identified sufficient conditions under which each system may exhibit multistationarity or is guaranteed to be monostationary, without requiring precise knowledge of rate constants. The key findings are summarized below.

● Among the four systems, OF is the one whose multistationary subsets are confined to a single Anderies class, the positive class ($$ R > 0 $$), with monostationary subsets identified in the negative, P-null, and Q-null classes. This confinement to one class does not mean the Anderies class alone determines the outcome, however: every case in Proposition 2, including the positive class, additionally requires sign conditions specific to the OF reactions, so the technology's own parameters remain part of the determination for OF exactly as they do for BCS, SCS, and WR.

● The BCS and SCS systems share the broadest reach of multistationarity across Anderies classes among the four strategies analyzed. Both exhibit multistationary subsets in the positive, P-null, and Q-null classes, with monostationarity guaranteed only in the negative class under appropriate sign conditions. This shared profile is a consequence of shared structural properties: both BCS and SCS are deficiency-one networks in which the CDR-specific reactions extend the natural carbon cycle without directly modulating its rates, and the DOA identifies parameter conditions under which multistationarity arises in classes that would be monostationary at the level of the Anderies subnetwork in isolation.

● The WR system exhibits the most parameter-dependent multistationarity behavior among the four systems analyzed. In contrast to the BCS, OF, and SCS systems, the WR system is monostationary in the positive Anderies class and multistationary in the negative class. The Q-null class contains both monostationary and multistationary subsets, depending on the sign of $$ (q_4 - q_5)/(p_5 - p_4) $$. This behavior is a consequence of the bidirectional coupling between the atmosphere and WR storage introduced by the kinetic orders $$ q_4, p_4, q_5 $$, and $$ p_5 $$, which gives rise to a new ratio governing the Q-null class behavior independently of the Anderies subnetwork class.

● The structural design of a CDR system is itself a risk-relevant parameter, and it can take qualitatively different forms across CDR strategies. The four systems analyzed here span a range from a deficiency-two network with no direct modulation of the natural cycle (OF), to deficiency-one networks whose CDR-specific reactions extend the natural cycle without modulating its rates while admitting multistationarity in additional classes through the conditions identified by the DOA (BCS and SCS), to a deficiency-two network with bidirectional coupling between the storage pool and the atmosphere (WR). Each of these forms produces a distinct multistationarity profile across Anderies classes, demonstrating that the network-level structure of a CDR system can either introduce or suppress tipping-point risk in ways that are not predictable from the Anderies subnetwork class alone.

● The RNCDR framework provides a structure-based route to screening for multistationarity, which is a necessary condition for an abrupt transition between steady states. By relying on the graphical and kinetic structure of the network rather than precise parameter values, the framework can characterize multistationarity conditions across all four CDR systems without exhaustive numerical simulation. The sign conditions on $$ p_1, p_2, q_1, q_2 $$ identified in Propositions 1–4 provide a concrete set of criteria that could guide the empirical estimation of kinetic orders from observed carbon flux data in future work.

For Earth system science, this study demonstrates that CRNT offers a tractable and systematic complement to conceptual and numerical carbon-cycle modeling for the analysis of steady-state multiplicity, and hence of one structural precondition for tipping-point behavior. In this framework, carbon pools are represented as interacting components and carbon transfers as reactions, and the question of whether the system can settle into more than one equilibrium is answered from the structure of these interactions rather than from numerical simulation. The RNCDR framework identifies structural drivers of multistationarity from the graphical and kinetic structure of simple heuristic box models, complementing the insights available from high-dimensional computational models. The framework is designed to accommodate other CDR technologies and portfolios of simultaneously deployed strategies, as noted in Section "Limitations", though demonstrating this extensibility remains a direction for future work.

For climate policy, the findings underscore that the reliability of a CDR strategy as a climate mitigation tool cannot be assessed independently of the dynamic regime of the natural carbon cycle into which it is deployed. The natural carbon cycle can operate in qualitatively different regimes, characterized here by the Anderies class, which determines whether the system is capable of settling into more than one equilibrium carbon distribution. A CDR system deployed into a multistationary regime may settle into a steady state characterized by low net CO₂ removal performance, undermining the carbon accounting assumptions underlying its deployment. The results suggest that a basic characterization of this background regime, through estimation of how photosynthesis and respiration respond to changes in carbon pool sizes, should be treated as a prerequisite screening step in CDR risk assessment, prior to the evaluation of any specific technology. Furthermore, the structural design of CDR technologies, specifically the nature of the feedback between the CDR storage pool and the natural carbon cycle, should be recognized as a risk-relevant design parameter in the development and evaluation of CDR deployment strategies. Aligning these findings with IPCC-style scenario analyses could provide further insights into the role of each CDR strategy in achieving long-term mitigation targets and inform the prioritization of technologies in national and global climate action plans^[12].

Several directions for future research emerge from this study. Empirical validation of the kinetic order conditions identified here will be essential to strengthen the credibility and applicability of the results. In particular, estimating the kinetic orders $$ p_1, p_2, q_1, q_2 $$ for the natural carbon cycle and the CDR-specific kinetic orders from real-world data would allow the Anderies class and the multistationarity conditions to be evaluated for specific geographical regions or deployment scenarios.

Extending the framework to portfolios of multiple CDR technologies deployed simultaneously is another important direction. The RNCDR framework is designed to accommodate such combinations, and analyzing the multistationarity behavior of CDR portfolios could reveal synergies or trade-offs that are not apparent from single-technology analyses. One case is of particular interest. BCS and SCS have multistationarity broadly distributed across Anderies classes, whereas that of OF is anchored more narrowly to a single class. Whether combining the two kinds reduces the overall tipping-point risk of a portfolio is an open question.

Finally, the model outcomes could be aligned with IPCC-style scenario analyses^[12]. Doing so would give insight into the role of each CDR strategy in meeting long-term mitigation targets, and would inform how technologies are prioritized in national and global climate action plans.