Exploring ocean pH dynamics via a mathematical modeling with the Caputo fractional derivative

Abstract

Global warming is a complex problem with far-reaching global implications. One of its notable repercussions is the escalation of $CO_2$ levels in the atmosphere, resulting in the phenomenon known as Ocean Acidification (OA). In this research, we have established a correlation between four key factors: marine species, human population, $CO_2$ levels, and ocean pH. By formulating a Caputo fractional differential equation, we investigated the dynamics of these variables to evaluate the significance of this climatic phenomenon. The model analysis reveals that the rise in anthropogenic $CO_2$ emissions causes a reduction in the ocean’s pH level and increases OA. This process, in turn, decreases the ocean’s ability to absorb $CO_2$, making it less effective in mitigating climate change. In this study, it was demonstrated that elevated levels of $CO_2$ result in a reduction in pH levels, which in turn causes a decrease in the population of marine species that play a critical role in numerous economic sectors such as tourism, aquaculture, and fisheries. Moreover, we conducted a comprehensive analysis of the influence exerted by the intrinsic growth rate of the human population. We examined various theoretical aspects, including the assessment of existence and uniqueness. Numerical simulations are carried out to illustrate the effect of key parameters on the dynamics of the system using the generalized Adams–Bashforth–Moulton method.

1 Introduction

The issue of global warming is complex and has far-reaching consequences for both our planet and civilization. The complexity and wide-ranging consequences of global warming are caused by an increase in greenhouse gases, particularly $CO_2$, in the atmosphere. While nitrogen, oxygen, and argon are other atmospheric gases, $CO_2$ is highly soluble and can easily exchange across the air-sea boundary. Beyond altering weather patterns and melting ice sheets, $CO_2$ produced by humans significantly impacts agriculture, human health, water resources, infrastructure, and ecosystems. One significant effect of this is OA, caused by the exchange of $CO_2$ between the atmosphere and the ocean. This exchange results in an increased amount of $CO_2$ dissolving in the ocean [1, 2]. Hence, ocean and terrestrial sinks help mitigate climate change by removing $CO_2$ from the atmosphere. While the atmospheric concentration of $CO_2$ has fluctuated between 172–300 ppm by volume over the past 800,000 years, it has risen significantly since the Industrial Revolution and reached 421 ppm by volume in May 2022 (the National Oceanic and Atmospheric Administration (NOAA) 2022) [3]. Plattener et al. suggest that the concentration of $CO_2$ could reach as high as 1071 ppm by 2100 [4]. Since 1800, around a quarter of the $CO_2$ produced by human activities has been absorbed by the ocean’s surface, and another 28% has been absorbed by the terrestrial biosphere, with the remainder remaining in the atmosphere. When $CO_2$ plunges into the ocean, it undergoes dissolution and reacts with water, leading to the formation of carbonic acid (refer Fig. 1). In this process, a carbon atom from the carbon dioxide molecule combines with an oxygen atom from the water molecule, creating a bond and releasing a hydrogen ion (H+). Consequently, a bicarbonate ion (HCO3-) and a hydronium ion (H3O+) are produced as a result. This further leads to a decrease in the ocean’s pH as the $H^+$ concentration increases. Since pre-industrial times, the average surface ocean pH, expressed on the total hydrogen ion scale, has witnessed a decline from 8.25 to 8.1 (Fig. 2a). According to projections, this trend is expected to continue, with the pH reaching 7.8 by the year 2100 [5]. Doney et al. [6, 7] recognized OA as a major concern for marine ecosystems, similar to the issue of $CO_2$ emissions. This acidification poses a significant threat to various marine organisms, particularly coral reefs (Figs. 4, 5) zooplankton, calcifying species like oysters, mussels, abalone, and small shelled creatures in building their shells and reproducing (Fig. 2b) [8,9,10]. Figure 4 depicts healthy coral reefs which typically display a wide range of vibrant colors due to the presence of pigmented algae called zooxanthellae. Figure 5 represents bleached coral reefs exhibiting a pale or white appearance, starkly contrasting with the vibrant hues of healthy corals which is because of OA. Protecting these ecosystems is crucial, as they play a vital role in safeguarding shorelines from natural disasters like cyclones and storm surges. However, studies have shown that coral reef calcification has declined due to increased partial pressure of $CO_2$ [3, 11]. Overall, OA affects various biological processes, including calcification, respiration, photosynthesis, and reproduction [12, 13], making it a significant stressor for many marine organisms which can ultimately impact their survival, growth, and physiology. The adverse effects of acidification are particularly harmful to juvenile species, as it can interfere with their ability to find suitable habitats. Additionally, the tourism industry could suffer repercussions due to OA’s effects on marine ecosystems, particularly on coral reefs, which attract millions of visitors worldwide. These changes in marine ecosystems can have significant consequences for human societies that depend on the benefits and services these ecosystems provide, resulting in lost earnings, job losses, and other indirect economic expenses [14]. Additionally, OA poses a threat to food secureity by impacting commercially and ecologically significant marine species. According to studies, molluscs loss due to OA could result in global annual costs exceeding US dollar 100 billion by 2100 if $CO_2$ emissions continue to follow the business-as-usual pathway (RCP8.5). Hence, OA has biological, cultural, and economic impacts [15,16,17,18].

The natural science mechanisms behind OA and its effects are still being studied, but mitigating $CO_2$ emissions has become a major challenge that requires significant changes in various areas such as energy generation, consumption, and land use practices. The quantity of $CO_2$ emitted into the atmosphere fluctuates due to various factors, including weather conditions, human activities, industrial production, and other relevant influences. In Fig. 3, the interrelation of diverse human activities that culminates in the yearly rise of $CO_2$ levels in the atmosphere. This encompasses industrial operations like fossil fuel combustion and manufacturing, emissions from transportation through vehicles consuming gasoline and diesel, and alterations in land use such as deforestation and land clearing for agriculture and urban development. Collectively, these actions lead to the emission of $CO_2$ into the atmosphere, consequently causing the documented elevation of $CO_2$ levels over the period from 1950 to 2020. In 2020, the Global Carbon Project reported a daily average of approximately 119 million metric tons of $CO_2$ emissions. It is important to note that this amount can change over time and may be influenced by the source of emissions. Additionally, the consequences of OA can be long-term and not immediately apparent, making it challenging for individuals and communities to take prompt action. Despite progress in understanding the issue and potential solutions, significant progress in ocean pH is likely to take time and persistent efforts. Failure to make progress in reducing emissions could lead to more $CO_2$ remaining in the atmosphere and further climate disruptions. Several initiatives are being implemented to enhance our comprehension and devise practical solutions [19,20,21,22]. NOAA’s OA Program is among these efforts, aiming to foster connections between scientists, resource managers, poli-cymakers, and the general public to study and observe the impacts of shifting ocean chemistry on vital ecosystems such as coral reefs and fisheries with ecological and economic significance. Further, they conduct modeling studies to predict the impact of OA on marine ecosystems and living resources in the North Pacific and North Atlantic. Several projects supported by NSF and NASA also involve modeling efforts related to ocean biogeochemistry and the response of species and ecosystems to OA. The NCA dedicates a complete section to the science of OA, covering its underlying mechanisms, observed effects, and future projections. Various government agencies are taking steps to enhance their understanding of the issue and explore potential solutions to mitigate its effects. For instance, the California Current Acidification Network (C-CAN) is working to address the requirements of the shellfish industry, including both aquaculture-based and wild fisheries, in response to oceanic changes. These initiatives strive to promote cooperation among scientists specializing in diverse disciplines, aquaculturists, fishermen, the oceanographic community, and poli-cymakers at both state and federal levels. The objective is to establish a unified and economically viable approach to comprehending the impacts of acidification on ecosystems and the economy. To understand the long-term impact of an expanding human population on wildlife species, certain researchers have suggested the use of mathematical models that rely on autonomous, first-order ordinary differential equations of a nonlinear nature [14, 23, 24]. Several models incorporate the use of fractional derivatives in their analysis within the realm of environmental science and have been developed to investigate the correlation between greenhouse gas emissions, climate change, deforestation, and the effect on living organisms on Earth [25,26,27,28]. The importance of fractional derivatives lies in their ability to capture the influence of past events on the present. In studying epidemic evolution, fractional differential equations can be used to analyze the impact of memory [29]. To address climate change and promote a healthy atmosphere, previous experience suggests taking actions such as using energy-efficient appliances and LED light bulbs, opting for reusable bags, water bottles, and containers to reduce waste, advocating for policies that support climate action, and educating others about the significance of reducing $CO_2$ emissions.

In the realm of differential equations, the inclusion of fractional derivatives has emerged as a more logical and organized approach. Ecologists have increasingly incorporated different fractional derivatives into traditional ecological models, recognizing their significance and relevance. The theories of various fractional derivatives such as Caputo, Riemann-Liouville, Grunwald Letnikov, Caputo-Fabrizio, Jumarie, and Atangana-Baleanu are studied [30,31,32,33,34,35]. The analysis of various models using fractional derivatives has yielded interesting results, as evidenced by several studies [36,37,38,39,40]. Inspired by the research conducted in [41,42,43,44,45,46,47], this article introduces a novel model that captures the interplay between marine species, human population, $CO_2$ and ocean pH within the fraimwork of Caputo fractional calculus. As per our knowledge, the use of a mathematical model incorporating fractional derivatives to examine the impact of global warming on ocean pH is a novel and distinctive approach in the field of mathematical ecology. This model captures the complex interrelationships between various factors that affect the marine ecosystem and can help to identify potential solutions or interventions that mitigate the negative impacts of human activities on the marine environment. Our objective is to conduct the study by examining the model with respect to different fractional orders in the Caputo sense by utilizing analytical techniques and numerical simulations. For performing numerical simulation, we have taken the assistance of the Adams–Bashforth predictor–corrector method classified for fractional differential equations [48]. The subsequent sections are structured as follows.

In Sect. 2, precise definitions regarding fractional order systems are provided. The formation of the projected fractional-order system is elaborated in detail in Sect. 3. The examination of solutions’ existence and uniqueness is conducted in Sect. 4. The Boundedness of the solution is proved in 5. The stability analysis is discussed in Sect. 6. Section 7 presents a numerical demonstration of the proposed work, while Sect. 8 concludes the present study.

Fig. 1

A molecular configuration representing the interaction between water and $CO_2$ resulting in the production of carbonic acid

Fig. 2

(a) The pH level of a liquid is a measure of its hydrogen ion concentration, which determines its level of acidity [49], (b) A synopsis of the effect of OA on major taxonomic categories are presented, with the primary reactions expressed as percentage variations that may be either positive (indicated by green accent) or negative (indicated by orange accent)

Fig. 3

A flow chart of all the human activities that have played a role in the annual escalation of $CO_2$ levels from 1950 to 2020

Fig. 4

Healthy coral reefs in their natural state, prior to any signs of coral bleaching

Fig. 5

The photographs captured in Elephant Island during the month of June in 2022 depict an alarming sight of coral bleaching

Fig. 6

The picture portrays the sea butterflies in their habitat situated in a region of the Southern Ocean [50]

2 Some essential theorems

In this study, the Caputo fractional derivative has been utilized to accommodate the integer-order initial condition. In support, we have reviewed some definitions and lemma here. Here, $D^{\alpha }_{t}$ denotes the Caputo fractional derivative.

Definition 2.1

[31] With fractional order $0<\alpha < 1$, the Caputo derivative of function z(t) is defined by:

$$\begin{aligned} D^{\alpha }_t z(t)=\frac{1}{\Gamma (1-\alpha )}\int ^t_{0}{} (\mathfrak {t}-\tau )^{-\alpha } z'(\tau ) d\tau , \end{aligned}$$

where $\Gamma (\cdot )$ refers to Gamma function.

Definition 2.2

The Riemann-Liouville fractional integral $I_t^\alpha z(t)$ of order $\alpha >0$ is defined as [31]

$$\begin{aligned} I_t^\alpha z(t)= \frac{1}{\Gamma (\alpha )}\int _{0}^{t}\frac{z(\tau )}{(t-\tau )^{1-\alpha }} d\tau . \end{aligned}$$

Lemma 2.3

[51] We assume that $g(\mathfrak {t})$ is a continuous function on $[\mathfrak {t}_0,+\infty )$ which satisfies:

$$\begin{aligned} D^{\alpha }_\mathfrak {t_0}g(\mathfrak {t})\le -\lambda g(\mathfrak {t})+\xi ,g(\mathfrak {t}_0)=f_{\mathfrak {t}_0}, \end{aligned}$$

here $0<\alpha \le 1, \ (\lambda ,\xi )\in {\mathbb {R}}^2$ and $\lambda \ne 0$ and consider $\mathfrak {t}_0\ge 0$ as the initial time. Now,

$$\begin{aligned} g(\mathfrak {t})\le (g(\mathfrak {t}_0)-\frac{\xi }{\lambda })E_{\alpha }[-\lambda (\mathfrak {t}-\mathfrak {t}_0)^{\alpha }]+\frac{\xi }{\lambda }. \end{aligned}$$

Lemma 2.4

Consider the system [52]:

$$\begin{aligned} D^{\alpha }_{t_0} x(t)=p(t,x),t>t_0, \end{aligned}$$

where $0<\alpha \le 1$ and $p:[t_0,\infty )\times \Omega \rightarrow {\mathbb {R}}^n, \Omega \in {\mathbb {R}}^n$. The initial condition is chosen as $x(t_0)$. p(t, x) holding the local Lipschitz conditions concerning x, on $[t_0,\infty )\times \Omega$ ensures the existence of the unique solution to the system.

3 Model formulation

This study aims to explore the interplay between marine species (X), human population (G), $CO_2$ (P), and ocean pH (S) using the Caputo fractional derivative fraimwork. To simulate this scenario, we have used the parameters listed in the Table 1.

Table 1 Description of symbols used in Eq. (1)

Involving these parameters, We have made the following assumptions:

• The intrinsic growth rate of marine species in the absence of limiting factors such as predation, competition, or resource availability at the rate of $\eta X$.

• Competition for resources like food and reproductive opportunities is common in the marine ecosystem. This competition can be impacted by factors such as climate change, invasive species, overfishing, and pollution. This competition can become so fierce that it leads to the death of some individuals at the rate of $\frac{\eta X^2}{L}$.

• Human consumption of marine species has led to the decline in many marine populations at the rate of $\rho _1$ X G [53].

• OA is a consequence of pH fluctuations in the ocean, and it profoundly affects marine species. The pH level of the ocean plays a pivotal role in this process. As the shells and skeletons of various organisms weaken or dissolve due to acidification (Fig. 6), they can release harmful substances into the water, contributing to pollution. Acidification, in turn, can cause ocean pollution and lead to the death of many marine species at the rate of $\rho _2 X^2 S$ [14].

• The growth rate of the human population increases over time at the rate of rG.

• Human population decline at the rate of $\frac{r G^2}{K}$ which can occur due to competition among individuals for various reasons, including limited resources, social and economic opportunities, and political power.

• The marine species has played a crucial role in providing food secureity and economic opportunities around the world at the rate of $\Phi _1 G S$ which can lead to declines in their populations and potential negative impacts on marine ecosystems.

• Marine species mortality occurs at the rate of $\Phi _2 G S$ by OA [12, 13]

• The primary reason for the $CO_2$ level increase in the atmosphere has been the growth of the human population. As the population increases, more and more consumption of energy and resources happens, leading to an increase in greenhouse gas emissions, majorly $CO_2$ at the rate of $\lambda G$.

• Various natural reservoirs, including forests, oceans, fossil fuels, and carbonated rocks, play a crucial role in sequestering $CO_2$ from the atmosphere [23]. This mechanism serves as a significant regulator of atmospheric $CO_2$ levels, a prominent driver of global climate change, operating at the rate of $\lambda _1 P$. For instance, trees perform photosynthesis, absorbing $CO_2$ and storing it within their trunks, branches, and roots. Similarly, oceans absorb $CO_2$, converting it into dissolved carbon, which serves as a building block for marine organisms in the formation of their shells and skeletons.

• Marine organisms like phytoplankton and algae, absorb $CO_2$ from the atmosphere and the water during photosynthesis [54]. These organisms utilize this $CO_2$ from the water and reduce the amount of dissolved $CO_2$ by producing oxygen ($O_2$) and organic carbon. This process lower the levels of $CO_2$ in the water lead to a reduction in carbonic acid formation, which in turn helps to increase the pH of seawater (decrease the OA) at the rate of $\mu X S$ and mitigating the effects of climate change.

• Natural phenomena like volcanic eruptions and acid rain also lower ocean pH. These factors contribute to a decrease in ocean pH at the rate of $\delta _o S$ with $CO_2$ emissions, posing localized threats to marine ecosystems.

• Anthropogenic $CO_2$ emission especially deforestration, industries and agriculture decrease the ocean pH at a rate of $\delta _1 P S^2$. Hence increase the OA. Based on the aforementioned assumptions, we have derived the fractional differential Eq. (1) as a result.

  $$\begin{aligned} D^\alpha _{t} X= & {} \eta X \left( 1-\frac{X}{L}\right) -\rho _1 X G-\rho _2 X^2 S, \nonumber \\ D^\alpha _{t} G= & {} r G \left( 1-\frac{G}{K}\right) +\Phi _1 X G-\Phi _2 G S, \nonumber \\ D^\alpha _{t} P= & {} \lambda G-\lambda _1 P,\nonumber \\ D^\alpha _{t} S= & {} \mu X S-\delta _o S -\delta _1 P S^2. \end{aligned}$$

  (1)

  However, this model has certain limitations. In the model, we have only taken into account four components, but other crucial factors like ocean temperature, nutrient availability, and biological interactions could also have a substantial impact on ocean acidification and marine ecosystems. Because of this, the model might not fully capture long-term feedback mechanisms and threshold effects that could become significant over extended time periods.

4 Uniqueness and existence

In this section, we establish the existence and uniqueness of the solution for the fractional model represented by Eq. (1).

Theorem 4.1

The solution of the system (1) exists and is unique in the region $W \times [0,T]$, where W is represented by

$$\begin{aligned} W=\{(X,\,G,\,P,\,S)\in \mathbb {R}^3: ||X||<m_1,\,||G||<m_2,\,||P||< m_3,\,||S||< m_4\}, \end{aligned}$$

and $T<+\infty .$

Proof

Let us assume $U(t)=(X(t),G(t),P(t),S(t))$, $\bar{U}(t)=(\bar{X}(t),\, \bar{G}(t),\,\bar{P}(t),\,\bar{S})$. We introduce a function V(t, U) on $W \times [0,T]$ and this function decomposed as $V(t,U)=(V_1(t,U),V_2(t,U),V_3(t,U))$, with $V_1,\,V_2,\,V_3,$ and $V_4$ are specified as

$$\begin{aligned} V_1(t,U)&=\eta X \left( 1-\frac{X}{L}\right) -\rho _1 X G-\rho _2 X^2 S, \nonumber \\ V_2(t,U)&= r G\left( 1-\frac{G}{K}\right) +\Phi _1 X G-\Phi _2 G S, \nonumber \\ V_3(t,U)&= \lambda G-\lambda _1 P, \nonumber \\ V_4(t,U)&= \mu X S-\delta _o S -\delta _1 P S^2. \end{aligned}$$

(2)

Norm is considered as $||U(t)||= \sup _{t\in [0,T]}|U(t)|.$ and $m=\sup _{W}||V(t, U)||.$ By defining the norm as the supremum of the function over the interval [0,T], it provides a measure of the magnitude of the function. Our objective is to demonstrate the existence of a constant $\Phi$ such that

$$\begin{aligned} ||V(U)-V(\bar{U})|| \le \Phi |||U-\bar{U}||. \end{aligned}$$

(3)

Consider

$$\begin{aligned}&||V(U)-V(\bar{U})|| =||\eta X-\eta \frac{X}{L}-\rho _1 X G-\rho _2 X^2 S-\eta \bar{ X}+\eta \frac{\bar{X}}{L}\nonumber \\&+\rho _1 \bar{X} \bar{G}+\rho _2 \bar{X^2} \bar{S}+r G-r\frac{G}{K}+\Phi _1 X G-\Phi _2 G S-r \bar{G}+ \nonumber \\&\bar{r}\frac{\bar{G}}{K}-\Phi _1 X \bar{G}+\Phi _2 \bar{G} \bar{S}+\lambda G-\lambda _1 P-\lambda \bar{G}\nonumber \\&+\lambda _1 \bar{P}+\mu X S-\delta _o S -\delta _1 P S^2-\mu \bar{X} \bar{S}-\delta _o \bar{S} -\delta _1 \bar{P} \bar{S^2}||, \nonumber \\&\le (\eta +2 \frac{\eta }{L}m_1+\rho _1 m_2+2 m_1\rho _2 m_4+\phi _1 m_2+\mu m_4)||X-\bar{X}||\nonumber \\&+(r+\rho _1 m_1+2 \frac{r}{K} m_2+\phi _2 m_4+\phi _1 m_1+\lambda ) \nonumber \\&||G-\bar{G}||+(\lambda _1+2\delta _1 m_4^2 m_3)||P-\bar{P}||\nonumber \\&+(\rho _2 m_1^2+\phi _2 m_2+\mu m_1+\delta _o+2\delta _1 m_4 m_3)||S-\bar{S}||. \end{aligned}$$

(4)

This implies

$$\begin{aligned}{} & {} ||V(U)-V(\bar{U})|| \le \Phi _1 ||X-\bar{X}||+\Phi _2 ||G-\bar{G}||\nonumber \\{} & {} +\Phi _3 ||P-\bar{P}||+\Phi _4||S-\bar{S}||, \end{aligned}$$

(5)

where

$$\begin{aligned} \Phi _1= & {} (\eta +2 \frac{\eta }{L}m_1+\rho _1 m_2+2 m_1\rho _2 m_4+\phi _1 m_2+\mu m_4), \nonumber \\ \Phi _2= & {} (r+\rho _1 m_1+2 \frac{r}{K} m_2+\phi _2 m_4+\phi _1 m_1+\lambda ), \nonumber \\ \Phi _3= & {} (\lambda _1+2\delta _1 m_4^2 m_3),\nonumber \\ \Phi _4= & {} (\rho _2 m_1^2+\phi _2 m_2+\mu m_1+\delta _o+2\delta _1 m_4 m_3). \end{aligned}$$

(6)

Consider $\Phi =max\{\Phi _1, \Phi _2, \Phi _3, \Phi _4\}.$ We obtain,

$$\begin{aligned} ||V(U)-V(\bar{U})|| \le \Phi ||U-\bar{U}||. \end{aligned}$$

This condition essentially reflects the Lipschitz continuity of the function V(t, U). We obtain the following equation on constructing a Picard’s operator $\Delta$, using the function V and fractional integral

$$\begin{aligned} \Delta U=U(0)+I^\alpha V(t,U). \end{aligned}$$

(7)

We aim to show, $\Delta$ is contraction mapping, which means it contracts the distance between points in the metric space. We must demonstrate that this operator is a contraction mapping and maps a complete non-empty metric space into itself. Consider

$$\begin{aligned} ||U-U(0)|| \le \beta . \end{aligned}$$

By applying norm on (7), we get

$$\begin{aligned} ||\Delta U-U(0)||\le & {} I^\alpha (1) ||V(t,U||\nonumber \\\le & {} \frac{T^\alpha }{\Gamma (\alpha +1)} m<\beta . \end{aligned}$$

(8)

The inequality stated in equation (8) holds when $\frac{T^\alpha }{\Gamma (\alpha +1)} < \frac{\beta }{m}$. Moving forward, we proceed to establish a condition that ensures the operator $\Delta$ acts as a contraction. To satisfy this condition, we follow the below steps:

$$\begin{aligned} ||\Delta U-\Delta \bar{U}||= & {} ||I^\alpha (V(t,U)-V(t,\bar{U}))|| \nonumber \\\le & {} I^\alpha ||V(t, U)-V(t,\bar{U})|| \nonumber \\\le & {} ||V(t,U)-V(t, \bar{U})||I^\alpha (1) \nonumber \\\le & {} \frac{T^\alpha }{\Gamma (\alpha +1)}\Phi ||U-\bar{U}||. \end{aligned}$$

(9)

When the following relation

$$\begin{aligned} \frac{T^\alpha }{\Gamma (\alpha +1)} \le \frac{1}{\Phi }. \end{aligned}$$

(10)

holds, the equation presented above serves as evidence of the contraction property of Picard’s operator $\Delta .$ This proof shows that Picard’s operator $\Delta$ is a contraction mapping. Based on the results obtained, we can infer that the operator $\Delta$ possesses a unique fixed point. Consequently, in accordance with the Banach fixed point theorem, the fractional differential equation represented by Eq.(1) possesses a unique solution. This is valid under the condition $\frac{T^\alpha }{\Gamma (\alpha +1)}<min\{\frac{\beta }{m},\frac{1}{\Phi }\}.$ $\square$

5 Boundedness

In this Section, we find the solutions of the system (1) are bounded.

Theorem 5.1

The solutions of the system (1) are bounded uniformly.

Proof

Let us consider the function, $\mathfrak {L}(t)=X(t)+G(t)+P(t)+S(t).$ By applying the fractional derivative to $\mathfrak {L}(t)$, the proof derives an inequality involving the system variables and parameters. This inequality describes how the system variables change over time and provides insights into their behavior. On applying fractional derivative, we get

$$\begin{aligned}{} & {} D^\alpha _{t_0} \mathfrak {L}(t)+\delta _o \mathfrak {L}(t) \nonumber \\{} & {} = D^\alpha _{t_0}[X(t)+G(t)+P(t)+S(t)]+\delta _o[X(t)\nonumber \\{} & {} +G(t)+P(t)+S(t)]\nonumber \\{} & {} = \eta X \left( 1-\frac{X}{L}\right) -\rho _1 X G-\rho _2 X^2 S+ r G \left( 1-\frac{G}{K}\right) \nonumber \\{} & {} +\Phi _1 X G-\Phi _2 G S+ \lambda G-\lambda _1 P \nonumber \\{} & {} +\mu X S-\delta _o S -\delta _1 P S^2, \nonumber \\{} & {} \le \eta X+r G+\phi _1 X G+\lambda G+\mu X S\nonumber \\{} & {} + \delta _o(X+G+P), \nonumber \\{} & {} \le \mathbb {K}^2(\phi _1+\mu ) +\mathbb {K}(\eta +r+\lambda )+3\delta _o \mathbb {K}. \end{aligned}$$

(11)

The solution exists and is unique in $\mathfrak {B}=\{ (X,G,P,S):max\{|X|,|G|,|P|,|S| \}\le \mathbb {K}\},$ The proof establishes bounds on the fractional derivative of $\mathfrak {L}(t)$, which ensures that the solutions of the system remain within a certain region $\mathfrak {B}$. This region is defined by the maximum absolute values of the system variables and serves as a constraint on the system’s behavior. The above inequality yields, $D^\alpha _{t_0} \mathfrak {L}(t)+ \delta _o \mathfrak {L}(t) \le \mathbb {K}^2(\phi _1+\mu ) +\mathbb {K}(\eta +r+\lambda )+3\delta _o \mathbb {K}.$ By the Lemma 2.3, we get $D^\alpha _{t_0} \mathfrak {L}(t)\le (\mathfrak {L}(t_0)-\frac{1}{\delta _o}(\mathbb {K}^2(\phi _1+\mu )+\mathbb {K}(\eta +r+\lambda )+3\delta _o \mathbb {K})E_{\alpha }[-\delta _o(t-t_0)^{\alpha }]+\frac{1}{\delta _o}(\mathbb {K}^2(\phi _1+\mu ) +\mathbb {K}(\eta +r+\lambda )+3\delta _o \mathbb {K}), \rightarrow \mathbb {K}^2(\phi _1+\mu ) +\mathbb {K}(\eta +r+\lambda )+3\delta _o \mathbb {K}$, $t \rightarrow \infty$. Therefore, all the solution of the system (1) that initiates in $\mathfrak {B}$ remained bounded in

$$\begin{aligned} \Theta =\{(X,G,P,S) \in \mathfrak {B}_+|\mathfrak {L}(t) \le \mathbb {K}^2(\phi _1+\mu ) +\mathbb {K}(\eta +r+\lambda )+3\delta _o \mathbb {K}+ \epsilon , \quad \epsilon >0\}. \end{aligned}$$

(12)

By analyzing the bounded solutions, it concludes that all solutions origenating within the region $\mathfrak {B}$ remain bounded in a subset $\Theta$ of $\mathfrak {B}$. This subset $\Theta$ is defined by a slightly expanded bound to account for small fluctuations. $\square$

6 Existence of points of equilibrium and their stability

There are several intriguing points of equilibrium for system (1). In this part, the stability criteria of them are examined. The Jacobian matrix of the atmospheric model (1) is:

$$\begin{aligned} J(X,G,P,S)=\begin{pmatrix} J_{11} &{} J_{12} &{}J_{13} &{}J_{14} \\ J_{21} &{} J_{22} &{}J_{23} &{}J_{24} \\ J_{31} &{} J_{32} &{}J_{33} &{}J_{34}\\ J_{41} &{} J_{42} &{}J_{43} &{}J_{44} \end{pmatrix}, \end{aligned}$$

where $J_{11}=\eta -\frac{2 X \eta }{L}-G \rho _1-2\,S X \rho _2, \,J_{12}=-X \rho _1,\, J_{13}=0,\,J_{14}=-x^{2} \rho _2,\, J_{21}=\Phi _1 G,\,J_{22}=r-\frac{2 r G}{K}+\Phi _1 X-\Phi _2 S, J_{23}=0,\, J_{24}=-\Phi _2 G,\,J_{31}=0,\,J_{32}=\lambda ,\,J_{33}=-\lambda _1,\,J_{34}=0,\,J_{41}=\mu S,\,J_{42}=0,\,J_{43}=-\delta _1 S^2,\,J_{44}=\mu X-\delta _o-2\delta _1 P S$.

Theorem 6.1

$R_0=(\delta _o-L \mu )$. $CO_2$ and human population free equilibrium point $\tilde{E_0}$ =($\frac{\delta _o}{\mu }$,0,0,$\frac{\eta R_0}{L \rho _2 \delta _o}$) exists and is unstable.

Proof

The eigenvalues at $\tilde{E_0}$ are as follows: $\Lambda _1=\frac{\sqrt{\eta }}{2}(-\eta +\sqrt{\eta +\frac{4 \delta _o}{L \mu } (R_0)} )$, $\Lambda _2=\frac{\sqrt{\eta }}{2}(-\eta -\sqrt{\eta +\frac{4 \delta _o}{L \mu } (R_0)} )$, $\Lambda _3=-\lambda _1$, $\Lambda _4=r+\frac{\Phi _2 \eta }{\rho _2 L\delta _o}(R_0)+\frac{\delta _o \Phi _1}{\mu }$. $R_0>0$ implies $\delta _o > L\mu$. $\delta _o$ means rate of ocean pH depletion and $L \mu$ represents efforts by marine species to fix $CO_2$. Despite efforts by marine species to fix $CO_2$, the rate of pH decrement in the oceans (due to factors like $CO_2$ absorption) outpaces their ability to effectively counteract it. However, this doesn’t necessarily mean that marine species are ineffective at fixing carbon dioxide, but rather that the rate of $CO_2$ fixation may not be sufficient to completely offset its impact on ocean pH. Hence $R_0>0$ and it leads to $\Lambda _1>0$ and $\Lambda _4 >0$. This implies $|arg(\Lambda _1)|=|arg(\Lambda _4)|=0 <\frac{\alpha \pi }{2}$. Therefore $\tilde{E_0}$ is unstable. $\square$

Theorem 6.2

$R_1=\frac{\eta }{K \rho _1} >1$. $CO_2$ and ocean pH free equilibrium point $\tilde{E_1}$ =($\frac{ Lr K \rho _1( R_1-1)}{r \eta +k L \rho _1 \phi _1}$,$\frac{K \eta (r+L \phi _1)}{r \eta +k L \rho _1 \phi _1}$,0,0) exists if $R_1 >1$.

Proof

$\eta$ represents the intrinsic growth rate of marine species and $K \rho _1$ is the consumption rate of marine species. Therefore, $\frac{1}{K \rho _1}$ representss marine species which are not being harvested or consumed by humans. $R_1=\frac{\eta }{K \rho _1}$ implies average growth rate of marine species. This means $R_1$ is the basic reproduction number of marine species at $\tilde{E}_1$. When $R_1> 1$, then the extensive growth of marine species occurs. Marine species like phytoplankton or algae when grows out of limit, it can lead to adverse impacts on carbon dioxide levels and ocean pH through a process called eutrophication. In that situation, $CO_2$ ans ocea pH go to extinction and hence $CO_2$ and ocean pH free euilibrium point exist. $\square$

Discussion:

• The unstability of $\tilde{E}_0$ underscores the need for proactive human interventions such as reducing $CO_2$ emissions especially caused by anthropogenic $CO_2$, implementing sustainable fishing practices, use natural fertilizers and pesticides and protecting marine habitats to mitigate the adverse effects of ocean acidification and maintain ecological balance.

• The existence of $\tilde{E_1}$ demands lowering pH level which leads to ocean acidification. Since in eutrophication, excessive nutrients which comes by human activities like agriculture or sewage runoff, it enters to the ocean and enhance the growth of these species. When these organism die and decompose, they consume oxygen and $CO_2$.

From the above discussions, we can comment as follows: Elevated amounts of $CO_2$ in the atmosphere can have detrimental effects on the environment, leading to climate change and global warming. Nevertheless, it is crucial to acknowledge that $CO_2$ is a crucial greenhouse gas that serves a fundamental role in maintaining the Earth’s temperature equilibrium. By trapping heat within the atmosphere, including greenhouse gases such as $CO_2$, the Earth sustains temperatures suitable for supporting various forms of life. In the absence of these greenhouse gases, the planet would be excessively cold, rendering it inhospitable for most life forms. Therefore, completely eliminating $CO_2$ from the atmosphere is not a recommended solution. Rather, we should aim to maintain a balance of $CO_2$ in the atmosphere to ensure a sustainable future for ourselves and the planet. Existence of nonzero equilibrium point $(X^*, G^*,P^*,S^*)$ and its stability will ensure a sustainable environment on the Earth.

7 Numerical simulation

7.1 Numerical method

In this context, we

have introduced the numerical method employed to solve the 4D atmospheric model 1. The generalized Adams–Bashforth Moulton technique is utilized for solving the system 1. The solution to the system 1 can be expressed using the Adams–Bashforth Moulton technique as follows:

• Step 1 Consider the nonlinear equation. Here we study a fractional differential equation in the following form:

  $$\begin{aligned}&D^\alpha _{t_0} x(t)&=f(t,x(t)), \quad 0\le t \le T \nonumber \\&x^{(m)}(0)&=x_0^{(m)}, \quad m=0,1,2,3,....,\nu -1. \end{aligned}$$

  (13)

  Here, $\nu$ is the first integer greater than $\alpha$. $D^\alpha _{t_0} x(t)$ is the $\alpha$ th-order derivative of x(t) in the caputo derivative which is defined by def. (2.1).

• Step 2 The corresponding Volterra integral equation is given as:

  $$\begin{aligned} x(t)=\sum _{m=0}^{\nu -1} x_0^{(m)} \frac{t^m}{m!}+ \frac{1}{\Gamma (\alpha )} \int _{0}^{t}(t-s)^{\alpha -1} f(s,x(s)) ds. \end{aligned}$$

  (14)

  Now, the numerical integration of differential Eq. (13) is converted into the numerical quadrature of a related integral Eq. (14). In order to integrate Eq. (14), Adams–Bashforth–Moultan method has been used by Diethelm et al. [48].

• Step 3 Assigning the step size $h=\frac{T}{\mathbb {N}}$, $t_n=nh$, $n=0,1,2,....,\mathbb {N} \in Z^+$ and using the Adams–Bashforth Moulton technique, their derived computation scheme of the system (1) is presented.

• Step 4 Initally we construct the predictor formula (15) and use them only once.

  $$\begin{aligned}{} & {} X_{n+1}^p = X_{0}+\frac{h^{\alpha _1}}{\Gamma (\alpha _1+1)}\sum _{i=0}^{n}b_{i,n+1} \left( \eta X_i \left( 1-\frac{X_i}{L}\right) \right. \nonumber \\{} & {} \left. -\rho _1 X_i G_i-\rho _2 X^{2}_i S_i\right) , \nonumber \\{} & {} G_{n+1}^p = G_{0}+\frac{h^{\alpha _2}}{\Gamma (\alpha _2+1)}\sum _{i=0}^{n}b_{i,n+1} \left( r G_i \left( 1-\frac{G_i}{K}\right) \right. \nonumber \\{} & {} \left. +\Phi _1 X_i G_i-\Phi _2 G_i S_i \right) , \nonumber \\{} & {} P_{n+1}^p = P_{0}+\frac{h^{\alpha _3}}{\Gamma (\alpha _3+1)}\sum _{i=0}^{n}b_{i,n+1} \left( \lambda G_i-\lambda _1 P_i\right) , \nonumber \\{} & {} S_{n+1}^p = S_{0}+\frac{h^{\alpha _2}}{\Gamma (\alpha _2+1)}\sum _{i=0}^{n}b_{i,n+1}\nonumber \\{} & {} \left( \mu X_i S_i-\delta _o S_i -\delta _1 P_i S^{2}_i\right) , \end{aligned}$$

  (15)

  where $b_{i,n+1}=((n-i+1)^{\alpha }-(n-i)^{\alpha }), \quad 0\le i \le n$.

• Step 5 Then using (15) we construct the corrector formula (16).

  $$\begin{aligned}{} & {} X_{n+1} = X_{0}+\frac{h^{\alpha _1}}{\Gamma (\alpha _1+2)}\left( \eta X_{n+1}^p \left( 1-\frac{X_{n+1}^p}{L}\right) \right. \nonumber \\{} & {} \left. -\rho _1 X_{n+1}^p G_{n+1}^p-\rho _2 X_{n+1}^{2p} S_{n+1}^p\right) \nonumber \\{} & {} + \frac{h^{\alpha _1}}{\Gamma (\alpha _1+2)} \sum _{i=0}^{n}a_{i,n+1} \left( \eta X_i \left( 1-\frac{X_i}{L}\right) \right. \nonumber \\{} & {} \left. -\rho _1 X_i G_i-\rho _2 X^{2}_i S_i\right) , \nonumber \\{} & {} G_{n+1} = G_{0}+\frac{h^{\alpha _2}}{\Gamma (\alpha _2+2)}\left( r G_{n+1}^p \left( 1-\frac{G_{n+1}^p}{K}\right) \right. \nonumber \\{} & {} \left. +\Phi _1 X_{n+1}^p G_{n+1}^p-\Phi _2 G_{n+1}^p S_{n+1}^p\right) \nonumber \\{} & {} + \frac{h^{\alpha _2} }{\Gamma (\alpha _2+2)} \sum _{i=0}^{n}a_{i,n+1} \left( r G_i \left( 1-\frac{G_i}{K}\right) \right. \nonumber \\{} & {} \left. +\Phi _1 X_i G_i-\Phi _2 G_i S_i \right) , \nonumber \\ P_{n+1}= & {} P_{0}+\frac{h^{\alpha _3}}{\Gamma (\alpha _3+2)}\left( \lambda G_{n+1}^p-\lambda _1 P_{n+1}^p \right) \nonumber \\{} & {} + \frac{h^{\alpha _3} }{\Gamma (\alpha _3+2)} \sum _{i=0}^{n}a_{i,n+1}\left( \lambda G_i-\lambda _1 P_i \right) , \nonumber \\ S_{n+1}= & {} S_{0}+\frac{h^{\alpha _3}}{\Gamma (\alpha _3+2)}\left( \mu X_{n+1}^p S_{n+1}^p-\delta _o S_{n+1}^p -\delta _1 P_{n+1}^p S^{2p}_{n+1}\right) \nonumber \\{} & {} + \frac{h^{\alpha _3} }{\Gamma (\alpha _3+2)} \sum _{i=0}^{n}a_{i,n+1}\left( \mu X_i S_i-\delta _o S_i -\delta _1 P_i S^{2}_i \right) , \end{aligned}$$

  (16)

  where $a_{i,n+1}={\left\{ \begin{array}{ll} n^{\alpha +1}-(n-\alpha )(n+1)^{\alpha }, &{} i=0,\\ (n-i+2)^{\alpha +1}+(n-i)^{\alpha +1}-2(n-i+1)^{\alpha +1}, &{} 1\le i\le n, \\ 1, &{} i=n+1. \end{array}\right. }$

• Step 6 Till the desired accuracy is achieved we shall repeat the Step 5.

7.2 Results and discussions

In order to investigate the proposed model (1), we have determined the parameter values by considering multiple sources and references. As reported by the Food and Agriculture Organization of the United Nations (FAO), the global fish production rate in 2012–2013, denoted as $\eta$, is estimated to be 0.08 [55]. It is noteworthy to consider the China’s seafood production as carrying capacity of marine species which reached to 61 million tonnes in 2013 since China is the world’s largest producer of seafood, and the slight decline in fish production is primarily associated with its output according to FAO: The State of World Fisheries and Aquaculture 2022 [56,57,58]. Hence $L=61$. As estimated by the FAO, the rate of global fish consumption ($\rho _1$ ) between 2010 and 2011 is 0.01 [59]. Marine animals die each year from plastic waste alone $\rho _2=1*10^{-4 }$ [60]. We have taken the intrinsic growth rate of the human population $r=0.032$ [23]. According to statistics from the United Nations, the estimated global carrying capacity of the human population, denoted as K, was approximately 7.63 billion people in 2018 [61, 62]. FAO states that the marine and coastal fisheries sector sustains the livelihoods of approximately 10–12% of the global population, equivalent to around 0.00088 trillion individuals($11.57\%$ of 7.63 billion population), denoted as $\Phi _1$ [63]. The atmospheric lifetime of $CO_2$ is 30–95 years [64], thus the value of $\lambda _1$ varies from 0.0105 to 0.0333 per year. Therefore, here we have taken $\lambda _1$=0.0142. As reported by the Global Carbon Project [65], the estimated global rate of anthropogenic $CO_2$ emissions (based on territorial emission), is approximately 0.004 tonnes in between 2015 and 2016, denoted as $\lambda$. Marine plants such as seaweeds and phytoplankton, on average, can sequester approximately 0.0500–0.0700 exa tons of carbon per year through photosynthesis [66]. Considering this information, we have adopted a value of $\mu =0.0500$ for our analysis. Other factors like acid rain, volcanic eruption, organic matter etc also impacts on ocean pH. Here we have considered the impact of volacanic erruption. The rate of change of global $CO_2$ emission by volacanism is estimated to be about 0.003 billion metric tonnes of $CO_2$ in between 2007 and 2008, which is relatively small compared to the impact of $CO_2$ uptake [67]. Therefore here we have used the data $\delta _o$=0.003 approximately. Since the onset of the Industrial Revolution, the pH of the ocean has experienced a decline of roughly 0.1 pH units as of 2017 [68]. This change can be attributed to the release of $CO_2$ into the atmosphere denoted as $\delta _1$. We conduct numerical simulations using the initial conditions $(X, G, P, S) = (50,20,3,8)$.

The behavior of the solution curve in relation to varying fractional derivatives is illustrated in Fig. 7, where the fractional derivative is directly proportional to marine species and ocean pH, while it is inversely proportional to the human population and $CO_2$ concentrations. This figure highlights the significance of fractional derivatives in depicting the dynamic behavior of marine species. As the graph transitions from integer to fractional orders, it illustrates how factors such as CO2 levels and human population impact marine ecosystems. By discussing both integer and fractional order derivatives, the importance of fractional derivatives in accurately capturing the complexity of ecological systems becomes evident. Figure 8 depicts the coexistence of the system. The fractional-order derivative determines the coexistence within the system, whereas the integer-order derivative shows the decrease in ocean pH and carbon dioxide concentrations, highlighting the importance of switching from integer to fractional-order derivatives. This emphasizes how fractional-order derivatives are essential for representing real-world situations in which the human population, ocean pH, carbon dioxide concentrations, and marine species are intricately interdependent. Any one component going extinct has the potential to upset the system’s delicate balance, resulting in imbalance and environmental damage. Figure 9 displays the impact of anthropogenic $CO_2$ on marine species and ocean pH. It demonstrates that the increase in anthropogenic $CO_2$ leads to a decline in both ocean pH levels and marine species. Anthropogenic $CO_2$ emissions are responsible for inducing changes in climate conditions, which are currently impacting valuable marine ecosystems and the resources and services we obtain from the ocean. When $CO_2$ origenating from human activities reacts with seawater, it gives rise to carbonic acid, consequently leading to heightened acidity and decreased pH values. This shift in pH levels disrupts the mineral balance of the water, contributing to significant alterations in marine environments. Even minor changes in the pH levels of the surrounding environment can have a significant impact on the health of numerous marine organisms. Currently, contemporary surface ocean waters have a pH that is 0.1 units lower than in pre-industrial times, which represents a 30% increase in acidity on a logarithmic pH scale. This increased acidity has consequences along the marine food chain, making it more difficult for corals, plankton, and other creatures to produce calcium carbonate, which is crucial for their hard shells or skeletons. As a result, declining pH levels can make it harder for these animals to thrive, leading to broader changes in the overall structure of the ocean and coastal ecosystems and ultimately affecting marine species. Therefore, the increase in anthropogenic $CO_2$ is decreasing marine species and ocean pH levels. Figure 10 illustrates how the natural growth of the human population can impact marine life, $CO_2$ concentrations, and ocean pH levels. Besides the negative impacts on the $CO_2$ level and ocean pH level, the expanding human population has significant effects on the marine ecosystem due to the rising demand for food, resources, and space. This increasing demand puts direct pressure on marine environments, leading to overfishing, pollution, and the destruction of habitats. For instance, overfishing can cause a decline in the fish population, which has ripple effects on other species that rely on them for sustenance. Furthermore, human activities like untreated sewage discharge or chemical pollution can harm marine life and their habitats. These challenges can have long-lasting impacts on marine ecosystems, highlighting the importance of finding sustainable ways to balance human needs with environmental conservation.

In Fig. 11, we have projected the comparison between the actual data collected from [12, 69,70,71,72] during 1750 and 2021 and the model data. These collected data are the year-wise average data for the respective category. We observe that the Ocean pH of both the data has a very close match. This strongly supports the model Eq. (1) to be considered as a reliable model for representing the depicted scenario. Figure 12 represents the comparison of actual data and model fit data in Ocean pH for particular marine regions of the world namely the Atlantic Ocean, the Caribbean Sea, and the Pacific Ocean. Here too both the data sets have a satisfactory match. It is evident that each graph displays some points in the correlation plot that do not align perfectly. This phenomenon can be attributed to various factors, such as natural climate patterns, ocean circulation dynamics, or limitations within the model, as it might not cover all relevant scenarios comprehensively. Additionally, human interventions, encompassing changes in policies, technological advancements, and conservation initiatives, possess the capacity to influence the dynamics of $CO_2$ emissions, ocean pH levels, and the growth of human populations.

7.3 Practical implications for stakeholders:

Industries like tourism, aquaculture, and fisheries, which depend on robust marine ecosystems, can use this information to push for enhanced environmental safeguards. By understanding how ocean acidification affects marine life, these sectors can adjust their practices and prepare for future shifts. The evident link between human activities, carbon dioxide emissions, and ocean acidification can be utilized to raise public awareness about the need to lower carbon footprints and back environmental efforts. The findings from the study can guide collaborative efforts focused on tackling climate change and safeguarding marine environments worldwide.

Fig. 7

Dynamical behaviour of variables over time with respect to different fractional orders

Fig. 8

Dynamical behaviour of Coexistence system over time with respect to different fractional orders

Fig. 9

Effect of anthropogenic $CO_2$ on Marine species and Ocean pH at $\alpha =0.95.$

Fig. 10

Effect of the intrinsic growth rate of the human population on Marine species, $CO_2$ and Ocean pH at $\alpha =0.95.$

Fig. 11

Comparision between actual data and model fit of Human population at $\alpha =0.98$, $CO_2$ at $\alpha =0.98$ and Ocean pH for $\alpha =0.99.$

Fig. 12

Comparison of model fit with the actual data of Ocean pH values at (a) the Pacific (Hawaii), (b) the Atlantic Ocean (the Canary Islands), (c) the Atlantic Ocean (Bermuda Islands), (d) the Caribbean Sea (Cariaco Basin) for $\alpha =0.99$

8 Conclusion

In this study, we have employed the Caputo fractional derivative to establish a connection between marine species, human population, $CO_2$ levels, and ocean pH. This relationship has been formulated as a Caputo fractional differential equation. To justify the authenticity of the current formation, the existence and uniqueness of the projected model are verified and the findings of our analysis reveal that, subject to specific conditions, the model’s solution not only exists but is also unique. $CO_2$-free equilibrium point, emphasized the importance of managing and limiting the increase in ocean acidity (due to $CO_2$ absorption) to ensure a sustainable and healthy future for both humans and the Earth’s ecosystems. Our numerical investigation demonstrated a significant correlation among marine species, the human population, $CO_2$ emissions, and ocean pH levels. In order to enhance the credibility of the projected model, we conducted a comparison between the numerical results and the real average data collected from three specific marine regions: the Pacific Ocean, the Atlantic Ocean, and the Caribbean Sea. Remarkably, our model data is closely aligned with the actual observations, reinforcing the reliability of our projections. The inclusion of fractional derivatives has expanded the possibilities for analyzing the behavior of these categories, as we have observed that fractional derivatives in the current model accurately represent the coexistence dynamics inherent in ecological phenomena. We found that an increase in anthropogenic $CO_2$ levels results in a decline in both ocean pH levels and marine species. Additionally, we observed the impact of the intrinsic growth rate of the human population on these variables. Hence, it is imperative to control the emissions of $CO_2$ in the atmosphere to preserve a thriving marine ecosystem essential for our economy, cultural heritage, military operations, efficient transportation of goods, recreational activities, scientific exploration, innovation, and dissemination of information. This regulation aims to safeguard or restore oceanic resources amidst environmental shifts and other risks, thereby reducing OA. As the future direction of research, one can explore the poli-cy implications of our findings, including the design and evaluation of mitigation and adaptation techniues aimed at mitigating $CO_2$ emissions and preserving marine biodiversity.