Method for predicting stability of atmosphere-ocean-biosphere phosphorus-oxygen-carbon cycle

Abstract

The present application relates to carbon cycle stability analysis technical field, especially to a kind of method for predicting atmosphere-ocean-biosphere phosphorus-oxygen-carbon cycle stability, comprising: obtaining phosphorus, oxygen, carbon cycle flux parameter and initial state data, construct three coupling dynamics models and obtain simulation data by solving;The simulation data is analyzed by nonlinear dynamics stability;Obtain the isotopic historical data of geological record and correspond with model state variable, form geochemical constraint condition;The statistical characteristics of simulation and observation data are extracted to determine the prior distribution range of key parameters, iterative sampling is solved, whether the difference between the solving result and observation data is judged to be converged, if not converged, adjust the range and resample, until convergence, obtain the optimized model and stability prediction result.The present application verifies the same stability problem from three different angles of analytical calculation, visualization analysis and numerical experiment, which improves the reliability of phosphorus-oxygen-carbon cycle system stability analysis result.

Description

Technical Field

[0001] This invention relates to the field of carbon cycle stability analysis technology, and in particular to a method for predicting the stability of the atmosphere-ocean-biosphere phosphorus-oxygen-carbon cycle. Background Technology

[0002] Research on the evolution of Earth's habitability is a crucial national scientific need and a frontier scientific issue of concern in basic geological and geochemical research. The surface Earth system can be subdivided into the atmosphere, hydrosphere, and biosphere. The complex interactions of matter and energy between these systems, as well as the biogeochemical cycles within them, all require specific numerical simulation methods for characterization. Among the many element cycles, the phosphorus cycle, oxygen cycle, and carbon cycle are the three most core cycles in the surface Earth system: phosphorus is the main limiting nutrient element over long geological timescales (millions of years); the oxygen cycle is one of the most important factors in the formation of a habitable Earth, with atmospheric oxygen content dominating the maximum saturation biomass of the biosphere and determining the rise and fall of life; the carbon cycle affects changes in atmospheric carbon dioxide concentration, thereby regulating global temperature, climate change, and the long-term evolution of biological habitats.

[0003] The three key elements mentioned above are coupled and mutually constrained in their cycles. Whether their coupled system remained stable throughout geological history and whether it would become unstable under disturbances has been a persistent but unresolved issue in the field. Previous studies have largely focused on characterizing the evolution of each cycle itself, lacking numerical simulation methods that integrate the phosphorus, oxygen, and carbon cycles into a unified dynamic framework and systematically analyze the stability of their coupled system. Furthermore, there is a lack of suitable quantitative indicators to effectively compare and constrain model simulation results with geochemical observation data preserved in the geological record. Consequently, the assessment of the stability of the phosphorus, oxygen, and carbon cycle coupled system during geological history lacks sufficient numerical simulation support and observational data verification.

[0004] In the field of numerical simulation and parameter inversion related to the carbon cycle, patent publication number CN114970184A discloses an assimilation method and system for simultaneously inverting high-resolution anthropogenic CO2 emissions and natural CO2 fluxes. This method prepares prior anthropogenic CO2 emissions and biomass combustion data, combines them with atmospheric transport models and observational data for assimilation and inversion, achieving high-resolution simultaneous estimation of carbon fluxes from both anthropogenic and natural sources. While this technical solution achieves inversion constraints on carbon fluxes based on observational data, its inversion object is limited to the single carbon cycle flux itself, and it does not address the issue of determining the overall system stability under the coupling effects of phosphorus, oxygen, and carbon cycles over geological timescales. Summary of the Invention

[0005] In view of this, the present invention provides a method for predicting the stability of the atmosphere-ocean-biosphere phosphorus, oxygen and carbon cycles, in order to solve the technical problems in the prior art of lacking specific numerical simulation methods to study the stability evolution of the coupling of the phosphorus cycle, oxygen cycle and carbon cycle, and lacking suitable statistical features to quantitatively characterize the differences between model simulation results and observation data.

[0006] The technical solution of this invention is achieved as follows: This invention provides a method for predicting the stability of the atmosphere-ocean-biosphere phosphorus, oxygen, and carbon cycle, comprising the following steps: S1. Obtain the mass flux parameters and initial state data of the phosphorus cycle, oxygen cycle and carbon cycle in the atmosphere-ocean-biosphere system. Construct a dynamic model of the phosphorus-oxygen-carbon cycle based on the mass flux parameters and initial state data. The dynamic model of the phosphorus-oxygen-carbon cycle includes a set of differential equations characterizing the coupling relationship between the phosphorus cycle, oxygen cycle and carbon cycle. S2. Solve the phosphorus-oxygen-carbon cycle kinetic model to obtain phosphorus-oxygen-carbon cycle simulation data; S3. Perform nonlinear dynamic stability analysis on the phosphorus-oxygen-carbon cycle simulation data to obtain the stability analysis results corresponding to the phosphorus-oxygen-carbon cycle simulation data. S4. Obtain historical isotope observation data from the geological record, establish a correspondence between the historical isotope observation data and the state variables of the phosphorus-oxygen-carbon cycle dynamics model, and obtain geochemical constraints for verifying the phosphorus-oxygen-carbon cycle simulation data and stability analysis results. S5. Based on geochemical constraints, extract statistical characteristics from the simulated data of the phosphorus-oxygen-carbon cycle and the historical observation data of isotopes, and determine the prior distribution range of key parameters in the dynamic model of the phosphorus-oxygen-carbon cycle based on the statistical characteristics. S6. Iteratively sample key parameters within the prior distribution range, substitute the sampled parameters into the phosphorus-oxygen-carbon cycle dynamics model for solution, and determine whether the preset convergence condition is met based on the difference between the solution result and the historical isotope observation data. If not, adjust the prior distribution range and resample until the preset convergence condition is met, and obtain the phosphorus-oxygen-carbon cycle dynamics model optimized by geochemical constraints and the corresponding prediction results of the atmospheric-oceanic-biosphere phosphorus-oxygen-carbon cycle stability.

[0007] Based on the above technical solution, preferably, in step S1, the differential equations in the phosphorus-oxygen-carbon cycle kinetic model are expressed as follows: ; ; ; Where TP represents the state variable of the marine phosphorus pool, TO represents the state variable of the atmospheric-oceanic oxygen pool, TA represents the state variable of the atmospheric-oceanic carbon pool, and t represents geological time. , , These represent minute changes in TP, TO, and TA, respectively. A quantity that characterizes a small change in geological time t. Input flux to the phosphorus cycle For phosphorus cycle output flux, orgb is organic carbon burial flux, carbb is inorganic carbon burial flux, Oweathering is oxidative weathering flux, Rdegassing is reducing gas consumption flux, Carbweathering is carbonate rock weathering flux, and orgweathering is organic carbon weathering flux. The organic carbon burial flux satisfies Where newp represents primary productivity, , These are preset empirical constants; The oxidation weathering flux satisfies Where n is a preset experience index, A function representing primary productivity newp as the independent variable. It represents an exponential function with the natural constant e as its base.

[0008] Based on the above technical solutions, preferably, the solution to the phosphorus-oxygen-carbon cycle kinetic model includes: The differential equations in the phosphorus-oxygen-carbon cycle dynamics model are solved iteratively using explicit or implicit numerical integration methods to obtain time series data of phosphorus pool state variables, oxygen pool state variables, carbon pool state variables, and the evolution of each cycle flux over geological time, which serve as the simulation data for the phosphorus-oxygen-carbon cycle.

[0009] Based on the above technical solution, preferably, step S3 includes: obtaining the Jacobian matrix at the equilibrium point of the phosphorus-oxygen-carbon cycle kinetic model, and determining the stability of the equilibrium point based on the eigenvalues ​​of the Jacobian matrix. The Jacobian matrix is ​​expressed as: ; In the formula, J represents the phosphorus-oxygen-carbon cycle kinetic model at the equilibrium point. The Jacobian matrix at the equilibrium point; TP represents the state variables of the marine phosphorus pool, TO represents the state variables of the atmospheric-oceanic oxygen pool, and TA represents the state variables of the atmospheric-oceanic carbon pool; each element in the matrix represents the value of the first-order partial derivatives of f, g, and h with respect to TP, TO, and TA, respectively. This indicates that all the above partial derivatives are evaluated at the equilibrium point; The system's operational trajectory in the phase space formed by the phosphorus pool state variable TP, the oxygen pool state variable TO, and the carbon pool state variable TA is visualized to obtain a phase diagram, and the stability state of the system is analyzed based on the phase diagram. A perturbation was applied to the phosphorus-oxygen-carbon cycle kinetic model to test the system's ability to recover from the perturbation. The stability analysis results are obtained by combining at least two of the following analysis methods: Jacobi matrix eigenvalue judgment, phase diagram visualization analysis, and disturbance rejection test.

[0010] Based on the above technical solutions, preferably, in the phase diagram, the position of the equilibrium point is determined by the intersection of three isosurfaces defined by the zero rate of change isosurfaces corresponding to the phosphorus pool state variable TP, the oxygen pool state variable TO, and the carbon pool state variable TA. The stability state is determined by the trend of the system's trajectory relative to the equilibrium point. If the system's trajectory converges directly from the initial state to the equilibrium point, the system is considered to be in a stable state. If the system's trajectory repeatedly intersects three isosurfaces with the equilibrium point as the center and gradually moves away from the equilibrium point, the system is considered to be in a spiral-progressive unstable state. Applying perturbations involves applying white noise perturbation of 5% to 15% of the total carbon pool state variable TA. The recovery time required for the system to return to its initial state before the perturbation is applied is measured after the perturbation ends. The strength of the system's resistance to perturbations is evaluated based on the length of the recovery time.

[0011] Based on the above technical solutions, preferably, step S4 includes: acquiring historical observation data of carbon isotopes, uranium isotopes, lithium isotopes, and oxygen isotopes in the geological record; establishing mass balance equations between the flux of the corresponding circulating reservoir material and the isotope value for carbon isotopes, uranium isotopes, and lithium isotopes, respectively, to obtain the correspondence between carbon isotopes and the oxygen cycle and the carbon cycle, the correspondence between uranium isotopes and the phosphorus cycle and the oxygen cycle, and the correspondence between lithium isotopes and the carbon cycle; For oxygen isotopes, a conversion relationship between isotope observations and sea surface temperature is established, and the correspondence between oxygen isotopes and the carbon cycle is obtained. This forms a multi-isotope orthogonal constraint condition that independently constrains the phosphorus cycle, oxygen cycle, and carbon cycle, serving as a geochemical constraint condition.

[0012] Based on the above technical solutions, preferably, the correspondence between carbon isotopes, uranium isotopes, and lithium isotopes and their corresponding cycles is characterized by isotope mass balance equations of the same form; let X represent the cyclic reservoir material corresponding to any one of the carbon, uranium, and lithium isotopes, and TX be the state variable corresponding to the amount of reservoir material X. The mass balance equation for the state variable TX is expressed as: ; in Let X be the input flux of the stored substance. Let X be the output flux of the stored substance X through an oxygen-deficient environment. The output flux of stored substance X through other environments; The isotopic value of the stored substance X is denoted as Establish the isotopic mass balance equation: ; in , , These are the isotopic characteristic values ​​corresponding to the input end, the hypoxic environment output end, and the other environment output ends, respectively, all of which are preset constants; Substituting the mass balance equation for the state variable TX into the isotope mass balance equation, we obtain the isotope values. Equation of evolution over time: ; in , These are preset constants, representing the isotopic characteristic values ​​and isotopic values ​​at the output terminals of hypoxic environments and other environments, respectively. The differences between them; Will The solution results are used as the model output, and correspondences are established with the historical observation data of carbon isotopes, uranium isotopes, and lithium isotopes, respectively, to obtain the geochemical constraints used to constrain the oxygen cycle and carbon cycle, phosphorus cycle and oxygen cycle, and carbon cycle.

[0013] Based on the above technical solutions, preferably, the correspondence between oxygen isotopes and the carbon cycle is characterized by a temperature conversion formula, which is expressed as follows: ; Where T is the sea surface temperature calculated from oxygen isotopes. These are oxygen isotope observations of apatite. T represents the oxygen isotope value of seawater, and a and b are preset empirical constants. The constraint of oxygen isotopes on the carbon cycle is achieved by comparing T with the temperature output by the phosphorus-oxygen-carbon cycle kinetic model.

[0014] Based on the above technical solutions, preferably, in step S5, the statistical characteristics include the mean, standard deviation, fluctuation amplitude, and fluctuation period of the phosphorus-oxygen-carbon cycle simulation data and the historical observation data of isotopes. For geological records with dating errors and missing data, dating uncertainty is added to the age framework. Historical isotope observation data are transformed into moving averages, maximum moving averages, or minimum moving averages that take into account dating uncertainty. Statistical features are then extracted to determine the prior distribution range.

[0015] Based on the above technical solutions, preferably, step S6 specifically includes: Latin hypercube sampling is performed on key parameters within the prior distribution range to obtain multiple sets of sampled parameter combinations; Each set of sampling parameter combinations is substituted into the phosphorus-oxygen-carbon cycle kinetic model for solution. The Euclidean distance between the solution and the historical observation data of the isotopes is calculated. The sampling parameter combinations with the smallest Euclidean distance are selected, and the prior distribution range is narrowed based on the selected sampling parameter combinations. Based on the reduced prior distribution range, perform Latin hypercube sampling and solve again, and determine whether the preset convergence condition is met. If not, repeat the process of reducing the prior distribution range and resampling until the preset convergence condition is met.

[0016] The present invention has the following advantages over the prior art: (1) This invention combines three analytical methods: Jacobi matrix eigenvalue judgment, phase diagram visualization analysis, and anti-disturbance test. By judging the signs of the real and imaginary parts of the Jacobi matrix eigenvalues ​​at the equilibrium point, and combining the convergence or divergence trend of the trajectory near the intersection of the zero rate of change isosurface in the phase space of the state variables, as well as the measurement results of the time required for the system to recover to the initial state after applying a disturbance, the same stability problem is mutually verified from three different perspectives: analytical calculation, visualization analysis, and numerical experiment. This avoids the judgment bias that may be caused by relying on a single analytical method and improves the reliability of the stability analysis results of the phosphorus-oxygen-carbon cycle system.

[0017] (2) This invention utilizes historical observation data of carbon isotopes, uranium isotopes, lithium isotopes and oxygen isotopes in the geological record to establish geochemical constraints corresponding to the phosphorus cycle, oxygen cycle and carbon cycle, respectively, forming multi-isotope orthogonal constraints that independently constrain the phosphorus-oxygen-carbon cycle system, and thereby determining the prior distribution range of the key parameters of the model, so that the model simulation results can be quantitatively compared with the actual geochemical observation data preserved in the geological record, overcoming the problem that previous phosphorus-oxygen-carbon cycle related models lacked effective means of verification with geochemical observation data, and improving the geochemical reliability of the key parameter values ​​and corresponding stability prediction results of the model.

[0018] (3) This invention achieves automated optimization and adjustment of key parameters of the phosphorus-oxygen-carbon cycle dynamics model by iterative sampling and solving within the prior distribution range of key parameters, and gradually narrowing the prior distribution range according to the difference between the solution results and the historical observation data of isotopes until the preset convergence condition is met. This reduces the dependence of parameter values ​​on human experience judgment, and makes the final phosphorus-oxygen-carbon cycle dynamics model and its corresponding stability prediction results simultaneously satisfy the dual constraints of dynamic coupling relationship and geochemical observation data.

[0019] (4) The method of the present invention can make a scientific evaluation of the evolution of the stability of key material cycles in different spheres, calculate the evolution characteristics of cycle stability under different climate-environment boundary conditions in geological history, and make scenario simulations. Attached Figure Description

[0020] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0021] Figure 1 This is a flowchart of the method for predicting the stability of the atmosphere-ocean-biosphere phosphorus, oxygen, and carbon cycles of the present invention. Figure 2 This is a schematic diagram illustrating the biogeochemical processes simulated in this invention. Figure 3 This is a schematic diagram simulating physical transmission according to the present invention; Figure 4 This is a schematic diagram of the statistical feature selection and parameter optimization process of the present invention; Figure 5 The phase diagram analysis results are used to test the stability of the solutions to the system of differential equations in this invention. Figure 6 The results of the perturbation-restoring force test are used to test the stability of the solution of the differential equation system in this invention. Figure 7 This is a schematic diagram comparing the simulation results of this invention with actual geochemical data (carbon isotopes); Figure 8 This is a schematic diagram comparing the simulation results of this invention with actual geochemical data (uranium isotopes); Figure 9 This is a schematic diagram illustrating the parameter optimization and reduction of the model-data Euclidean distance for this invention. Detailed Implementation

[0022] The technical solutions of the present invention will be clearly and completely described below with reference to the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of the present invention.

[0023] like Figure 1-3 As shown, the present invention provides a method for predicting the stability of the atmosphere-ocean-biosphere phosphorus, oxygen, and carbon cycle, characterized by comprising the following steps: S1. Obtain the mass flux parameters and initial state data of the phosphorus cycle, oxygen cycle and carbon cycle in the atmosphere-ocean-biosphere system. Construct a dynamic model of the phosphorus-oxygen-carbon cycle based on the mass flux parameters and initial state data. The dynamic model of the phosphorus-oxygen-carbon cycle includes a set of differential equations characterizing the coupling relationship between the phosphorus cycle, oxygen cycle and carbon cycle. S2. Solve the phosphorus-oxygen-carbon cycle kinetic model to obtain phosphorus-oxygen-carbon cycle simulation data; S3. Perform nonlinear dynamic stability analysis on the phosphorus-oxygen-carbon cycle simulation data to obtain the stability analysis results corresponding to the phosphorus-oxygen-carbon cycle simulation data. S4. Obtain historical isotope observation data from the geological record, establish a correspondence between the historical isotope observation data and the state variables of the phosphorus-oxygen-carbon cycle dynamics model, and obtain geochemical constraints for verifying the phosphorus-oxygen-carbon cycle simulation data and stability analysis results. S5. Based on geochemical constraints, extract statistical characteristics from the simulated data of the phosphorus-oxygen-carbon cycle and the historical observation data of isotopes, and determine the prior distribution range of key parameters in the dynamic model of the phosphorus-oxygen-carbon cycle based on the statistical characteristics. S6. Iteratively sample key parameters within the prior distribution range, substitute the sampled parameters into the phosphorus-oxygen-carbon cycle dynamics model for solution, and determine whether the preset convergence condition is met based on the difference between the solution result and the historical isotope observation data. If not, adjust the prior distribution range and resample until the preset convergence condition is met, and obtain the phosphorus-oxygen-carbon cycle dynamics model optimized by geochemical constraints and the corresponding prediction results of the atmospheric-oceanic-biosphere phosphorus-oxygen-carbon cycle stability.

[0024] In one embodiment of the present invention, step S1 includes: obtaining the mass flux parameters and initial state data of the phosphorus cycle, oxygen cycle and carbon cycle in the atmosphere-ocean-biosphere system, and constructing a phosphorus-oxygen-carbon cycle dynamic model based on the mass flux parameters and initial state data.

[0025] Construct the differential equation for the phosphorus cycle: ; In the formula, For the state variables of the marine phosphorus pool, For atmospheric-oceanic oxygen pool state variables, Here, t represents the state variables of the atmospheric-oceanic carbon pool, and t represents geological time. , , These represent minute changes in TP, TO, and TA, respectively. A quantity that characterizes a small change in geological time t. Input flux to the phosphorus cycle orgb is the phosphorus cycle output flux, orgb is the organic carbon burial flux, Oweathering is the oxidative weathering flux, and Rdegassing is the reducing gas consumption flux.

[0026] In the phosphorus cycle, the phosphorus flux from rivers to the ocean State variables of atmospheric-oceanic oxygen pool State variables of the atmospheric-oceanic carbon pool The function of phosphorus output flux from the ocean State variables of the marine phosphorus pool oxygen pool state variables The function is that the input end of the phosphorus cycle is jointly regulated by the oxygen cycle and the carbon cycle, and the output end is jointly regulated by the phosphorus pool itself and the oxygen cycle, thus establishing the coupling relationship between the phosphorus cycle and the oxygen cycle and the carbon cycle.

[0027] Construct the differential equation for the oxygen cycle: ; ; ; ; ; In the formula, n is a preset empirical exponent, with a value of 0.5. This equation describes the basic input and output fluxes of the oxygen cycle. Organic carbon burial (orgb) is the main oxygen source, and its flux is a function of TO and TP. Oxygen-involved weathering and reducing gas consumption (Rdegassing) are the main outputs of the oxygen pool. Oweathering is a function of the relative size of the oxygen pool to the power of 0.5, and Rdegassing is a function of the relative size of the carbon pool.

[0028] Construct the differential equation for the carbon cycle: ; ; ; ; ; ; In the formula, carbb represents the inorganic carbon burial flux, and newp represents primary productivity. , As a preset empirical constant, The value is 2.1. The value is 0.7. A function representing primary productivity newp as the independent variable. Characterized by an exponential function with the natural constant e as the base. In the carbon cycle, the reducing gas consumption flux. Carbonate rock weathering flux Organic carbon weathering flux Together, they constitute the input end of the carbon pool, with the exhaust of mantle or deep material being an important source of atmospheric carbon dioxide; both carbonate weathering flux and organic carbon weathering flux are functions of the size of the sedimentary carbon pool, terrestrial runoff, and terrestrial vegetation. Characterizing the effect of sedimentary carbon pool size on flux regulation, Characterizing the regulatory effect of land runoff on flux, Characterizing the regulatory effect of terrestrial plant biomass on flux; the output end of the carbon pool is determined by organic carbon burial flux. The inorganic carbon burial flux (carbb) is composed of organic carbon burial flux, which is considered a function of oxygen pool size and primary productivity. The inorganic carbon burial flux and carbonate weathering flux maintain a dynamic equilibrium, i.e., carbb ≈ This ensures that the carbon cycle system is completely closed, with no additional carbon loss or increase.

[0029] Material flux parameters and initial state data can be determined based on the specific geological history period studied, referring to the range of values ​​in published geochemical databases, sedimentary record inversion results, or existing biogeochemical cycle models (such as the GEOCARB model, COPSE model, etc.). This invention does not limit the specific values.

[0030] In one embodiment of the present invention, step S2 includes: solving the phosphorus-oxygen-carbon cycle kinetic model by iteratively solving the differential equations in the model using explicit or implicit numerical integration methods to obtain time series data of the phosphorus pool state variable TP, oxygen pool state variable TO, carbon pool state variable TA, and the evolution of each cycle flux with geological time t, which serve as the phosphorus-oxygen-carbon cycle simulation data. Explicit numerical integration methods can employ methods such as the fourth-order Runge-Kutta method, while implicit numerical integration methods can employ implicit methods such as the backward difference formula. The specific solution can be implemented using general scientific computing software packages (e.g., the scipy.integrate module in Python, the ode solver in MATLAB, or the DifferentialEquations.jl package in Julia). The present invention does not limit the specific solution tools.

[0031] In one embodiment of the present invention, step S3 includes: obtaining the Jacobian matrix at the equilibrium point of the phosphorus-oxygen-carbon cycle kinetic model, and judging the stability of the equilibrium point based on the eigenvalues ​​of the Jacobian matrix; visualizing the system's trajectory in the phase space formed by the phosphorus pool state variable TP, the oxygen pool state variable TO, and the carbon pool state variable TA to obtain a phase diagram, and analyzing the stability state of the system based on the phase diagram; applying a perturbation to the phosphorus-oxygen-carbon cycle kinetic model and testing the system's resistance to perturbation recovery; and combining at least two of the analytical methods in the Jacobian matrix eigenvalue judgment, phase diagram visualization analysis, and perturbation resistance test to obtain the stability analysis result.

[0032] Specifically, the differential equations of the phosphorus-oxygen-carbon cycle kinetic model are denoted as follows: , , Assume that the system of differential equations has an equilibrium point. ,satisfy We perform substitution on the states near the equilibrium point, letting... , , This represents applying a small perturbation to the equilibrium point. Performing Taylor expansions on the substituted variables u, v, and w, retaining the first-order terms, and taking the equation corresponding to u as an example, the expansion yields... ; In the formula, , , The functions representing the time-varying rates of TP, TO, and TA with the three state variables are respectively expressed as the net mass flux expressions for the phosphorus cycle, oxygen cycle, and carbon cycle. , , TP, TO, and TA represent the values ​​of TP, TO, and TA at the equilibrium point, respectively; u, v, and w represent the small deviations of TP, TO, and TA from the equilibrium point, respectively, which represent a small perturbation applied to the equilibrium point. , , These represent the first-order partial derivatives of function f with respect to TP, TO, and TA at the equilibrium point, respectively, which are the sensitivity of f to each state variable.

[0033] because And higher-order terms Near the equilibrium point, the coefficients can be ignored, and the above equation can be approximated as a first-order linear equation; the equations corresponding to v and w can be derived in the same way. Arranging the first-order coefficients of the three equations u, v, and w row-wise yields the Jacobian matrix at the equilibrium point: ; In the formula, J represents the phosphorus-oxygen-carbon cycle kinetic model at the equilibrium point. The Jacobian matrix at the equilibrium point; each element in the matrix represents the value of the first-order partial derivatives of f, g, and h with respect to TP, TO, and TA, respectively, with subscripts... This indicates that all the partial derivatives mentioned above are evaluated at the equilibrium point.

[0034] The eigenvalues ​​of the Jacobian matrix J at the equilibrium point are obtained by computer program. The stability of the equilibrium point can be determined by the signs of the real and imaginary parts of the eigenvalues: if the real parts of all eigenvalues ​​are negative, the disturbance decays with time, the system converges back to the equilibrium point, and it is determined to be stable; if there is an eigenvalue with a positive real part, the disturbance amplifies with time and is rejected by the equilibrium point, and it is determined to be unstable; if there is a non-zero imaginary part among the eigenvalues, it indicates that the system oscillates during the process of approaching or deviating from the equilibrium point.

[0035] In phase diagram visualization analysis, isosurfaces with zero rate of change (dTP / dt=0, dTO / dt=0, dTA / dt=0) are plotted for the phosphorus pool state variable TP, the oxygen pool state variable TO, and the carbon pool state variable TA, respectively. The unique intersection of the three isosurfaces is the location of the equilibrium point. If the system's trajectory falls on a certain isosurface, it means that the corresponding state variable has reached equilibrium; if the trajectory enters the intersection of two isosurfaces, it means that both state variables have reached equilibrium simultaneously; if the trajectory finally reaches the intersection of all three isosurfaces, it means that all three state variables are in complete equilibrium. If the system's trajectory gradually approaches and converges directly to the equilibrium point from the initial state along the isosurfaces, the system is considered to be in a stable state; if the system's trajectory repeatedly intersects the three isosurfaces around the equilibrium point and gradually moves away from the equilibrium point, the system is considered to be in a spirally asymptotically unstable state.

[0036] In the disturbance immunity test, a white noise disturbance of 5% to 15% of the total carbon pool state variable TA is applied. The recovery time required for the system to return to its initial state before the disturbance is measured. A longer recovery time indicates a weaker disturbance immunity, while a shorter recovery time indicates a stronger disturbance immunity. The amplitude of this disturbance can be analogous to the random fluctuations in the eruption rate of volcanic activity during geological history, for example, approximately 2 × 10⁻⁶. 10 Flux fluctuations on the order of mol / yr simulate background eruption rate noise. Within the perturbation amplitude range of 5% to 15%, lower amplitude perturbations usually do not cause a fundamental change in system stability. However, when the perturbation amplitude increases to the upper part of this range (e.g., close to or exceeding 10%), it may lead to a collapse in system stability. Based on this, the system's ability to resist perturbations varies with the intensity of the perturbation by applying perturbations of different amplitudes in stages.

[0037] It should be noted that the judgment criteria for the three methods—Jacobi matrix eigenvalue judgment, phase diagram visualization analysis, and disturbance rejection test—are consistent (all based on the behavior of the system deviating from or converging to the equilibrium point after being disturbed). The combination of the three methods can mutually verify each other from three perspectives: analysis, visualization, and numerical experiment, avoiding misjudgments that may occur with a single method and improving the reliability of stability analysis results.

[0038] In one embodiment of the present invention, step S4 includes: acquiring historical observation data of carbon isotopes, uranium isotopes, lithium isotopes, and oxygen isotopes in the geological record; establishing mass balance equations between the flux of corresponding circulating reservoir materials and isotope values ​​for carbon isotopes, uranium isotopes, and lithium isotopes, respectively, to obtain the correspondence between carbon isotopes and the oxygen cycle and carbon cycle, the correspondence between uranium isotopes and phosphorus cycle, and the correspondence between lithium isotopes and carbon cycle; establishing a conversion relationship between isotope observation values ​​and sea surface temperature for oxygen isotopes, to obtain the correspondence between oxygen isotopes and carbon cycle; thereby forming multi-isotope orthogonal constraint conditions that independently constrain the phosphorus cycle, oxygen cycle, and carbon cycle, which serve as geochemical constraint conditions.

[0039] Specifically, the correspondences between carbon isotopes, uranium isotopes, and lithium isotopes and their corresponding cycles are all characterized by isotope mass balance equations of the same form. Let X represent the cyclic reservoir material corresponding to any one of the carbon, uranium, or lithium isotopes. It is important to note that when different elements are substituted in the isotope mass balance calculation, the calculation of input and output fluxes will differ, and the types of input and output fluxes will also change. However, the principle of isotope mass balance calculation remains the same. Here, we simply assume a basic case with one input flux and two output fluxes. TX is the state variable corresponding to the amount of reservoir material X, and the mass balance equation for the state variable TX is expressed as: ; in, Let X be the input flux of the stored substance. Let X be the output flux of the stored substance X through an oxygen-deficient environment. Let X be the flux of the reservoir material X through other environments. The isotopic value of the reservoir material X is denoted as... Establish the isotopic mass balance equation: ; in, , , These are the isotopic eigenvalues ​​corresponding to the input, hypoxic environment output, and other environment output, respectively, all of which are preset constants. Expanding the left side of the above equation using the product rule, we get: ; Substituting the mass balance equation for the state variable TX into the above equation, expanding it using the product rule, and combining like terms, we obtain the isotope values. Equation of evolution over time: ; in, , These are preset constants, representing the isotopic characteristic values ​​of the hypoxic environment output terminal, other environmental output terminals, and the (atmosphere-)oceanic reservoir isotopic values, respectively. The differences between them. The solution results are used as the model output and compared with the historical observation data of carbon isotopes, uranium isotopes and lithium isotopes respectively, so as to achieve constraints on the oxygen cycle and the carbon cycle, phosphorus cycle and carbon cycle.

[0040] Taking uranium isotopes as an example, the above general formula can be further explained: the output of uranium in the environment is a function of the combined effects of the phosphorus cycle and the oxygen cycle, that is, the output flux in an oxygen-deficient environment. The output fluxes are functions of the phosphorus pool state variable TP and the oxygen pool state variable TO, and other environmental output fluxes. This can be considered a fixed value, thus establishing a functional relationship between the uranium cycle and the phosphorus and oxygen cycles. Uranium isotopes can then be used as geochemical indicators of the phosphorus cycle and changes in ocean redox states. It should be noted that, unlike the differential equations in the phosphorus-oxygen-carbon cycle kinetic model, the uranium cycle is controlled by the phosphorus and oxygen cycles, while none of the cycles in the phosphorus-oxygen-carbon cycle is inversely affected by the uranium cycle. Therefore, the uranium cycle is not a state variable in the phosphorus-oxygen-carbon cycle kinetic model, but only serves as an indicative isotope system for the coupling relationship between the phosphorus and oxygen cycles. The uranium reservoir material is denoted as TU, and the isotope value is denoted as... Substituting the aforementioned general formula yields the equation for the evolution of uranium isotope values ​​over time, which can be compared with historical uranium isotope observation data from marine carbonate rocks. The specific functional dependencies between carbon isotopes, lithium isotopes, and their corresponding cycles can be determined by analogy to the uranium isotope example above, based on the actual input-output flux of their respective cycles. The derivation process of their mass balance equations is completely consistent with that of uranium isotopes, and will not be elaborated upon further in this invention.

[0041] The correspondence between oxygen isotopes and the carbon cycle is characterized by a temperature conversion formula, which is expressed as: ; Where T is the sea surface temperature calculated from oxygen isotopes. These are oxygen isotope observations of apatite. Here, T represents the oxygen isotope value of seawater, and a and b are preset empirical constants. In one embodiment of the present invention, a can be 117.4 and b can be 4.50. The specific values ​​should be determined according to the type of apatite mineral selected and the applicable temperature conversion standard. By comparing T with the temperature output by the phosphorus-oxygen-carbon cycle kinetic model, the constraint of oxygen isotopes on the carbon cycle is realized.

[0042] The Euclidean distances between the model outputs of the four isotopes (carbon, oxygen, uranium, and lithium) and the corresponding historical observation data are calculated separately and then directly added together to obtain the overall distance between the model and the actual observation data. This serves as a comprehensive index to measure the coupling degree between the phosphorus-oxygen-carbon cycle dynamics model and the geochemical observation data. Since the four isotopes (carbon, oxygen, uranium, and lithium) independently constrain the oxygen cycle and carbon cycle, and the phosphorus cycle and oxygen cycle, the orthogonal constraints of multiple isotope systems can reduce the ambiguity of single geochemical indicators, reduce the degrees of freedom of key model parameters, and improve the reliability of subsequent parameter optimization steps.

[0043] In one embodiment of the present invention, such as Figure 4 As shown, step S5 includes: extracting statistical features from the phosphorus-oxygen-carbon cycle simulation data and isotope historical observation data based on geochemical constraints. The statistical features include the mean, standard deviation, fluctuation amplitude, and fluctuation period of the phosphorus-oxygen-carbon cycle simulation data and isotope historical observation data. For geological records with dating errors and missing data, dating uncertainty is added to the age framework, and the isotope historical observation data is converted into a moving average, moving maximum, or moving minimum value after considering the dating uncertainty. Then, statistical features are extracted to determine the prior distribution range.

[0044] After clarifying the combination of statistical characteristics, the key parameters in the phosphorus-oxygen-carbon cycle kinetic model are classified into two categories: the first category consists of parameters with strong geochemical constraints, such as empirical constants controlling the input and output fluxes and isotopic characteristic values ​​at the empirical input and output ends. The prior distribution range of these parameters can be directly constrained based on published geochemical data. The second category consists of initial hypothesis parameters with relatively weak constraints that require further verification through model-data comparison, such as the magnitude of changes in terrestrial plant biomass and the magnitude of changes in ocean anoxic area. For these parameters, a relatively wide prior uncertainty (e.g., 15%–30%) can be given as the initial prior distribution range. For geological records with dating errors, the dating uncertainty interval can be superimposed on the sample age, and the moving average, moving maximum, or moving minimum value within this interval can be used to replace the single-point observation value. Then, statistical characteristics such as mean, standard deviation, fluctuation amplitude, and fluctuation period can be extracted from this interval to avoid distortion in the comparison between simulated and observed data caused by dating errors.

[0045] In one embodiment of the present invention, step S6 includes: performing Latin hypercube sampling on key parameters within the prior distribution range to obtain multiple sets of sampling parameter combinations; substituting each set of sampling parameter combinations into the phosphorus-oxygen-carbon cycle kinetic model for solution, calculating the Euclidean distance between the solution result and the historical observation data of the isotopes, selecting the sampling parameter combinations with the smallest Euclidean distance, and narrowing the prior distribution range based on the selected sampling parameter combinations; re-performing Latin hypercube sampling and solving based on the narrowed prior distribution range, determining whether a preset convergence condition is met, and repeating the process of narrowing the prior distribution range and re-sampling if not met, until the preset convergence condition is met.

[0046] Specifically, considering that the time taken for a single model run is typically on the order of 100 to 1000 milliseconds, the number of Latin hypercube sampling operations can be set to around 5000 to 10000 times with the help of parallel computing. After each round of sampling, the top 10% of the sampled parameter combinations with the smallest Euclidean distance are selected to reduce the uncertainty range of the corresponding parameters (for example, after the first round of selection, the uncertainty can be reduced from 15% to 20% to 10% to 12%), and the reduced upper and lower limits are used as the prior distribution range for the next round of Latin hypercube sampling. This process is repeated iteratively. The preset convergence condition can be set to the uncertainty reduction of the prior distribution range obtained from the previous and next iterations being less than 10%, or the weighted uncertainty of all key parameters being less than 5%. If either of the above conditions is met, the iteration stops, and the final result of parameter optimization is obtained.

[0047] The phosphorus-oxygen-carbon cycle dynamic model constructed in steps S1-S6 above couples the phosphorus, oxygen, and carbon cycles into a unified set of differential equations for modeling. It comprehensively evaluates the nonlinear dynamic stability of the system by combining multiple methods such as Jacobi matrix eigenvalue judgment, phase diagram visualization analysis, and anti-disturbance testing. This changes the previous situation where most studies focused on single cycle evolution processes and lacked coupled stability analysis methods. At the same time, by introducing orthogonal constraints of carbon, oxygen, uranium, and lithium multi-isotopes, and combining iterative parameter optimization processes of statistical feature extraction and Latin hypercube sampling, it achieves quantitative comparison and parameter verification between model simulation results and geochemical observation data in the geological record, thereby improving the geochemical reliability of stability prediction results.

[0048] The technical solution of the present invention will be described in more detail below with reference to specific embodiments.

[0049] Example 1 This embodiment provides a method for predicting the stability of the atmosphere-ocean-biosphere phosphorus, oxygen, and carbon cycle. The main body comprises a system of differential equations consisting of three differential equations corresponding to the phosphorus pool state variable TP, the oxygen pool state variable TO, and the carbon pool state variable TA. In this embodiment, inorganic carbon isotopes from marine carbonate rocks are introduced (…). As a geochemical indicator of the temporal changes in the carbon and oxygen pools, this is established through the isotopic mass balance equation between carbon isotopes and corresponding circulating reservoir materials in step S4 above. The correspondence between the state variables and the phosphorus-oxygen-carbon cycle kinetic model is established in the same way as the isotope mass balance equation described in step S6.

[0050] After solving the system of differential equations consisting of the three differential equations mentioned above, a three-dimensional phase diagram consisting of the phosphorus pool state variable TP, the oxygen pool state variable TO, and the carbon pool state variable TA can be drawn, as shown below. Figure 5 As shown, the three-dimensional phase diagram presents the zero-rate-of-change isosurfaces corresponding to the three differential equations. These three isosurfaces intersect at a unique point, which is the equilibrium solution of the system of differential equations. The stability of the equilibrium solution can be determined based on the geometric differences in the zero-rate-of-change isosurfaces: if the system's trajectory gradually approaches and converges to the intersection point along the isosurfaces, it is considered stable; if the trajectory repeatedly intersects the isosurfaces and gradually moves away from the intersection point, it is considered unstable. Furthermore, a disturbance recovery test can be performed on the system, measuring the time it takes for the phosphorus-oxygen-carbon cycle system to return to its initial stable state after being disturbed. The longer the time, the worse the system's disturbance recovery stability. Figure 6As shown in the figure, within a certain time interval, the stacking flux drops significantly and even reaches zero, then gradually recovers to the level before the disturbance and remains stable. This process of drop and recovery reflects the dynamic trajectory of the system gradually recovering to the initial stable state after being disturbed. The length of time taken for the drop and recovery is the recovery time index for measuring the system's ability to resist disturbances. The longer the recovery time, the weaker the system's ability to resist disturbances.

[0051] Example 2 This embodiment provides a method for predicting the stability of the atmosphere-ocean-biosphere phosphorus, oxygen, and carbon cycles. The main body also includes a system of differential equations consisting of three differential equations corresponding to the phosphorus pool state variable TP, the oxygen pool state variable TO, and the carbon pool state variable TA. Building upon Embodiment 1, this embodiment further adds other geochemical indicators to the original system of three differential equations, such as introducing uranium isotopes from marine carbonate rocks to indicate changes in the phosphorus cycle and the ocean's redox state. The specific method for establishing this system is consistent with the uranium isotope mass balance equation described in step S6.

[0052] In this embodiment, the phosphorus-oxygen-carbon cycle dynamics model must simultaneously calculate the Euclidean distance between the carbon, oxygen, uranium, and lithium isotope results output by the model and the corresponding historical observation data in the actual geological record. Compared to using only a single isotope index for constraint, introducing multiple isotope systems for coupled analysis provides stricter constraints on the parameters of the phosphorus-oxygen-carbon cycle dynamics model, reduces the ambiguity of single geochemical indices, and reduces the degrees of freedom of key model parameters. Specifically, for the application of geochemical observation data in the prediction of phosphorus-oxygen-carbon cycle stability, a combination method suitable for the time series statistical characteristics of isotope geochemical data is proposed to achieve more effective parameter optimization.

[0053] Specifically, the key to selecting suitable geochemical index time series for parameter optimization lies in finding geochemical indices that can independently describe the phosphorus, oxygen, and carbon cycles, forming orthogonal constraints on the respective subsystems of the phosphorus, oxygen, and carbon cycles. This ensures that the Euclidean distance calculated by combining these indices represents relatively independent and separate components in phase space, thereby avoiding reduced constraint effectiveness due to coupling redundancy between indices. In this embodiment, the specific method for calculating the Euclidean distance is as follows: uranium isotopes indicate the phosphorus cycle, carbon isotopes indicate the oxygen and carbon cycles, lithium isotopes indicate the weathering process mainly controlled by the carbon cycle, and oxygen isotopes indicate temperature changes mainly controlled by the carbon cycle. The Euclidean distances for carbon, oxygen, and lithium isotopes were calculated separately for each of the four isotopes (carbon, oxygen, uranium, and lithium) in the phosphorus-oxygen-carbon cycle dynamics model and the corresponding historical observation data in the actual geological record, forming a set of four Euclidean distances. These four Euclidean distances were then normalized by dividing each isotopic distance by one-thousandth (‰) of the corresponding isotopic value in the modern ocean, converting them to dimensionless values. The four normalized Euclidean distances were then summed to obtain the final composite distance between the phosphorus-oxygen-carbon cycle dynamics model and the geochemical observation data. This composite distance serves as an indicator of the coupling degree between the model parameters and the geochemical observation data. Figure 7-9 As shown, where Figure 7 The figure shows a comparison between carbon isotope simulation results and observational data. The simulation curves in the figure generally exhibit an oscillation pattern similar to that of the observational data (including several obvious negative shifts and subsequent recovery processes), indicating that the model can reproduce the carbon cycle perturbation events reflected by the carbon isotope record well. The branches that appear at the end of the curve (the difference between the solid and dashed lines) reflect the divergence of the model output results under different parameter combinations or different prediction scenarios. Figure 8 The figure shows a comparison between uranium isotope simulation results and observational data. The simulation curve reproduces several significant negative shifts and step-like changes in the observational data. These features usually indicate redox state changes such as the expansion of ocean anoxic area. This shows that the simulation curve can reflect the coupling changes between the phosphorus cycle and ocean redox state indicated by the uranium isotope record. The solid and dashed branches at the end of the curve also reflect the differences in model output results under different parameter combinations. Figure 9The figure illustrates the parameter optimization process. The solid black line and the dashed black line represent the simulated curves of the model output under two different parameter combinations (simulation result 1 and simulation result 2), respectively. Both of them have certain deviations from the observed data points. The red curve is the optimal simulation result obtained after parameter optimization and screening. Compared with simulation result 1 and simulation result 2, it fits the distribution trend of the observed data points better. This indicates that the Euclidean distance corresponding to this parameter combination is the smallest, that is, the coupling degree between the model output and the geochemical observation data is optimal. The figure intuitively shows the optimization process of obtaining the optimal parameter combination by iteratively sampling and screening parameters and gradually reducing the model-data Euclidean distance.

[0054] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for predicting the stability of the atmosphere-ocean-biosphere phosphorus, oxygen, and carbon cycle, characterized in that, Includes the following steps: S1. Obtain the mass flux parameters and initial state data of the phosphorus cycle, oxygen cycle and carbon cycle in the atmosphere-ocean-biosphere system. Construct a dynamic model of the phosphorus-oxygen-carbon cycle based on the mass flux parameters and initial state data. The dynamic model of the phosphorus-oxygen-carbon cycle includes a set of differential equations characterizing the coupling relationship between the phosphorus cycle, oxygen cycle and carbon cycle. S2. Solve the phosphorus-oxygen-carbon cycle kinetic model to obtain phosphorus-oxygen-carbon cycle simulation data; S3. Perform nonlinear dynamic stability analysis on the phosphorus-oxygen-carbon cycle simulation data to obtain the stability analysis results corresponding to the phosphorus-oxygen-carbon cycle simulation data. S4. Obtain historical isotope observation data from the geological record, establish a correspondence between the historical isotope observation data and the state variables of the phosphorus-oxygen-carbon cycle dynamics model, and obtain geochemical constraints for verifying the phosphorus-oxygen-carbon cycle simulation data and stability analysis results. S5. Based on geochemical constraints, extract statistical characteristics from the simulated data of the phosphorus-oxygen-carbon cycle and the historical observation data of isotopes, and determine the prior distribution range of key parameters in the dynamic model of the phosphorus-oxygen-carbon cycle based on the statistical characteristics. S6. Iteratively sample key parameters within the prior distribution range, substitute the sampled parameters into the phosphorus-oxygen-carbon cycle dynamics model for solution, and determine whether the preset convergence condition is met based on the difference between the solution result and the historical isotope observation data. If not, adjust the prior distribution range and resample until the preset convergence condition is met, and obtain the optimized phosphorus-oxygen-carbon cycle dynamics model and the corresponding prediction results of the atmospheric-oceanic-biosphere phosphorus-oxygen-carbon cycle stability.

2. The method for predicting the stability of the atmosphere-ocean-biosphere phosphorus, oxygen, and carbon cycle as described in claim 1, characterized in that, In step S1, the differential equations in the phosphorus-oxygen-carbon cycle kinetic model are expressed as follows: ； ； ； Where TP represents the state variable of the marine phosphorus pool, TO represents the state variable of the atmospheric-oceanic oxygen pool, TA represents the state variable of the atmospheric-oceanic carbon pool, and t represents geological time. , , These represent minute changes in TP, TO, and TA, respectively. A quantity that characterizes a small change in geological time t. Input flux to the phosphorus cycle For phosphorus cycle output flux, orgb is organic carbon burial flux, carbb is inorganic carbon burial flux, Oweathering is oxidative weathering flux, Rdegassing is reducing gas consumption flux, Carbweathering is carbonate rock weathering flux, and orgweathering is organic carbon weathering flux. The organic carbon burial flux satisfies Where newp represents primary productivity, , As a preset empirical constant, A function representing primary productivity newp as the independent variable. Characterized by an exponential function with the natural constant e as its base; The oxidation weathering flux satisfies , where n is the preset experience index.

3. The method for predicting the stability of the atmosphere-ocean-biosphere phosphorus, oxygen, and carbon cycle as described in claim 2, characterized in that, The solution to the phosphorus-oxygen-carbon cycle kinetic model includes: The differential equations in the phosphorus-oxygen-carbon cycle dynamics model are solved iteratively using explicit or implicit numerical integration methods to obtain time series data of phosphorus pool state variables, oxygen pool state variables, carbon pool state variables, and the evolution of each cycle input flux with geological time, which are used as the simulation data of the phosphorus-oxygen-carbon cycle.

4. The method for predicting the stability of the atmosphere-ocean-biosphere phosphorus, oxygen, and carbon cycle as described in claim 1, characterized in that, Step S3 includes: obtaining the Jacobian matrix at the equilibrium point of the phosphorus-oxygen-carbon cycle kinetic model, and determining the stability of the equilibrium point based on the eigenvalues ​​of the Jacobian matrix. The Jacobian matrix is ​​expressed as: ； In the formula, J represents the phosphorus-oxygen-carbon cycle kinetic model at the equilibrium point. The Jacobian matrix at the equilibrium point; TP represents the state variables of the marine phosphorus pool, TO represents the state variables of the atmospheric-oceanic oxygen pool, and TA represents the state variables of the atmospheric-oceanic carbon pool; each element in the matrix represents the value of the first-order partial derivatives of f, g, and h with respect to TP, TO, and TA, respectively. This indicates that all partial derivatives are evaluated at the equilibrium point; The system's operational trajectory in the phase space formed by the phosphorus pool state variable TP, the oxygen pool state variable TO, and the carbon pool state variable TA is visualized to obtain a phase diagram, and the stability state of the system is analyzed based on the phase diagram. A perturbation was applied to the phosphorus-oxygen-carbon cycle kinetic model to test the system's ability to recover from the perturbation. The stability analysis results are obtained by combining at least two of the following analysis methods: Jacobi matrix eigenvalue judgment, phase diagram visualization analysis, and disturbance rejection test.

5. The method for predicting the stability of the atmosphere-ocean-biosphere phosphorus, oxygen, and carbon cycle as described in claim 4, characterized in that, In the phase diagram, the position of the equilibrium point is determined by the intersection of three isosurfaces defined by the zero rate of change isosurfaces corresponding to the phosphorus pool state variable TP, the oxygen pool state variable TO, and the carbon pool state variable TA. The stability state is determined by the trend of the system's trajectory relative to the equilibrium point. If the system's trajectory converges directly from the initial state to the equilibrium point, the system is considered to be in a stable state. If the system's trajectory repeatedly intersects three isosurfaces with the equilibrium point as the center and gradually moves away from the equilibrium point, the system is considered to be in a spiral-progressive unstable state. Applying perturbations involves applying white noise perturbation of 5% to 15% of the total carbon pool state variable TA. The recovery time required for the system to return to its initial state before the perturbation is applied is measured after the perturbation ends. The strength of the system's resistance to perturbations is evaluated based on the length of the recovery time.

6. The method for predicting the stability of the atmosphere-ocean-biosphere phosphorus, oxygen, and carbon cycle as described in claim 1, characterized in that, Step S4 includes: acquiring historical observation data of carbon isotopes, uranium isotopes, lithium isotopes, and oxygen isotopes in the geological record; establishing mass balance equations between the flux of the corresponding circulating reservoir material and the isotope value for carbon isotopes, uranium isotopes, and lithium isotopes, respectively, to obtain the correspondence between carbon isotopes and the oxygen cycle and the carbon cycle, the correspondence between uranium isotopes and the phosphorus cycle and the oxygen cycle, and the correspondence between lithium isotopes and the carbon cycle; For oxygen isotopes, a conversion relationship between isotope observations and sea surface temperature was established, and the correspondence between oxygen isotopes and the carbon cycle was obtained. This results in multi-isotope orthogonal constraint conditions that independently constrain the phosphorus cycle, oxygen cycle, and carbon cycle, serving as geochemical constraint conditions.

7. The method for predicting the stability of the atmosphere-ocean-biosphere phosphorus, oxygen, and carbon cycle as described in claim 6, characterized in that, The correspondences between carbon isotopes, uranium isotopes, and lithium isotopes and their corresponding cycles are all characterized by isotope mass balance equations of the same form; let X represent the cyclic reservoir material corresponding to any one of the carbon, uranium, or lithium isotopes, and TX be the state variable corresponding to the amount of cyclic reservoir material X. The mass balance equation for the state variable TX is expressed as: ； in Let X be the input flux of the stored substance. Let X be the output flux of the stored substance X through an oxygen-deficient environment. The output flux of stored substance X through other environments; The isotopic value of the stored substance X is denoted as Establish the isotopic mass balance equation: ； in , , These are the isotopic characteristic values ​​corresponding to the input end, the hypoxic environment output end, and the other environment output ends, respectively, all of which are preset constants; Substituting the mass balance equation for the state variable TX into the isotope mass balance equation, we obtain the isotope values. Equation of evolution over time: ； in , These are preset constants, representing the isotopic characteristic values ​​and isotopic values ​​at the output terminals of hypoxic environments and other environments, respectively. The differences between them; Will The solution results are used as the model output, and correspondences are established with the historical observation data of carbon isotopes, uranium isotopes, and lithium isotopes, respectively, to obtain the geochemical constraints used to constrain the oxygen cycle and carbon cycle, phosphorus cycle and oxygen cycle, and carbon cycle.

8. The method for predicting the stability of the atmosphere-ocean-biosphere phosphorus, oxygen, and carbon cycle as described in claim 6, characterized in that, The correspondence between oxygen isotopes and the carbon cycle is characterized by a temperature conversion formula, which is expressed as: ； Where T is the sea surface temperature calculated from oxygen isotopes. These are oxygen isotope observations of apatite. T represents the oxygen isotope value of seawater, and a and b are preset empirical constants. The constraint of oxygen isotopes on the carbon cycle is achieved by comparing T with the temperature output by the phosphorus-oxygen-carbon cycle kinetic model.

9. The method for predicting the stability of the atmosphere-ocean-biosphere phosphorus, oxygen, and carbon cycle as described in claim 1, characterized in that, In step S5, the statistical characteristics include the mean, standard deviation, fluctuation amplitude, and fluctuation period of the phosphorus-oxygen-carbon cycle simulation data and the historical isotope observation data; For geological records with dating errors and missing data, dating uncertainty is added to the age framework. Historical isotope observation data are transformed into moving averages, maximum moving averages, or minimum moving averages that take into account dating uncertainty. Statistical features are then extracted to determine the prior distribution range.

10. The method for predicting the stability of the atmosphere-ocean-biosphere phosphorus, oxygen, and carbon cycle as described in claim 1, characterized in that, Step S6 specifically includes: Latin hypercube sampling is performed on key parameters within the prior distribution range to obtain multiple sets of sampled parameter combinations; Each set of sampling parameter combinations is substituted into the phosphorus-oxygen-carbon cycle kinetic model for solution. The Euclidean distance between the solution and the historical observation data of the isotopes is calculated. The sampling parameter combinations with the smallest Euclidean distance are selected, and the prior distribution range is narrowed based on the selected sampling parameter combinations. Based on the reduced prior distribution range, perform Latin hypercube sampling and solve again, and determine whether the preset convergence condition is met. If not, repeat the process of reducing the prior distribution range and resampling until the preset convergence condition is met.

Images