The classical NPZD approach is an improvement of the NPZ (Evans et al., 1985; Franks, 2002) which considers a single variable for each compartment (nutrients - phytoplankton - zooplankton - detritus, sometimes including a different compartment for bacteria) neglecting the differences between the different biological types (Fasham et al., 1990; Fasham, 1993; McGillicuddy et al., 1995b). In this approach, nitrogen is typically considered a limiting nutrient because it is a necessary component of chlorophyll, nucleic acids, and proteins, all of which are important for the growth of phytoplankton and photosynthesis (Liefer et al., 2019). Also, compared to other nutrients, nitrogen availability is frequently limited, which can have an adverse effect on the entire food web by limiting phytoplankton development (Bristow et al., 2017). The detritus component accounts for the organic matter pool, which is derived from faecal materials, the fraction non-assimilated by grazing zooplankton and phytoplankton decay. The detritus is recycled by bacteria and degraded to dissolved organic nitrogen, or digested by zooplankton (Leles et al., 2016). The process through which autotrophic organisms, including phytoplankton, produce organic matter by converting sunlight and nutrients into biomass referred as primary production. Contrarily, in secondary production, organic matter is consumed by heterotrophic creatures like zooplankton, which graze on primary producers or consume other organisms.

Modifications have been made in the classical NPZD approaches to improve the estimation of nitrogen and carbon flux (Fennel et al., 2006), with the addition of variable Chl-a and introducing a nitrogen-based nutrient-phytoplankton-heterotroph model which is of intermediate complexity with respect to Fasham’s and McGillicuddy’s models (Fasham et al., 1990; McGillicuddy et al., 1995a). Although nitrate is considered one of the most important limiting nutrients for the growth of phytoplankton biomass and primary production in open ecosystems; in some locations, such as the equatorial Pacific, North Pacific, and the Southern Ocean, nitrate (NO₃) levels may exceed 4 µM (central Pacific) and phytoplankton biomass is low (Chai et al., 1996). The high nitrate and low Chl-a (HNLC) conditions observed in these regions are commonly associated with limited iron availability (Leonard et al., 1999; Pitchford and Brindley, 1999; Venables and Moore, 2010).

Nevertheless, in high-latitude areas, other factors besides the iron deficiency may also inhibit phytoplankton growth. Low temperatures (Rhee and Gotham, 1981), imbalanced nutritional stoichiometry (Gerhard et al., 2022), short days and prolonged cloud cover (ChrEilertsen and Degerlund, 2010), and grazing pressure (Dale et al., 1999) from higher trophic levels are some examples of these constraints. Leonard et al., 1999 accounted for iron limitation (nutrient compartment) in his extended NPZD model with the employment of nine compartments. The model incorporated iron partitioning into new iron and regenerated iron, as well as dissolved iron bioavailability in the marine environment. Following these developments, Galbraith developed the Biogeochemistry with Light Iron Nutrients and Gases (BLING) model, also based on NPZD (Galbraith et al., 2009). This model employs a specific approach to isolate the global impact of iron on two major features of photosynthesis: maximum light-saturated photosynthesis rates and photosynthetic efficiency. The model uses sophisticated parameterizations or algorithms to separate and quantify the impact of iron on these distinct components. These sophisticated formulations take into account various factors and processes, such as iron uptake mechanisms, iron limitation effects, and biochemical reactions, to provide a detailed understanding of how iron affects photosynthetic process. The model also considers an implicit representation (without explicitly modeling it as a separate parameter or variable) of the growth rate of phytoplankton.

The proportion of the main elements in most biogeochemical models is commonly described with a constant elemental ratio (stoichiometry): 106 C (carbon):16 N (nitrogen):1 P (phosphorous) (Redfield, 1933). This ratio is the most widely used stoichiometric reference for planktonic nutrient limitation (Ptacnik et al., 2010). Hence, biogeochemical models may consider simple growth formulations with constant C: Nutrients stoichiometry or neglect photo-acclimation by considering a constant Chl:C ratio. The classical Michaelis Menten equations are the simplest representation of phytoplankton growth (Monod, 1949) assuming a constant stoichiometric ratio between nitrogen, phosphorous, and carbon (Redfield, 1933; Ho et al., 2003; Quigg et al., 2003; Quigg et al., 2011; Finkel et al., 2006; Finkel et al., 2006). The Redfield ratio has been modified to take into account observed differences from the conventional ratio. As study by (Anderson and Sarmiento, 1994) have been conducted to analyse nutrient data and evaluate the variability of nutrient ratios under different conditions in several regions. Their research emphasizes the significance of taking into account context-specific elements that can affect nutrient stoichiometry, such as nutrient availability, microbial community composition, and ecological interactions Biogeochemical models that incorporate these modified formulations of the Redfield ratio use specific units (nitrogen, phosphorus, or carbon currency) to represent phytoplankton biomass and these models can generally be utilized for global scale studies (Aumont and Bopp, 2006; Follows et al., 2007a; Dutkiewicz et al., 2009). On the other hand (Droop, 1983), has more advanced and complex growth formulations which account for the dynamics of phytoplankton cell internal quota (Klausmeier et al., 2004; Bougaran et al., 2010; Bernard, 2011; Mairet et al., 2011; Wang et al., 2022). Thus, phytoplankton biomass in Droop’s model (Droop, 1983) is represented by at least two state variables, typically carbon and nitrogen currency. Different versions of Droop’s model have been embedded and utilized in 1D (Lancelot et al., 2000; Allen et al., 2002; Lefe et al., 2003; Mongin et al., 2003; Blackford et al., 2004; Salihoglu et al., 2008) and 3D (Tagliabue and Arrigo, 2005; Vichi et al., 2007; Vichi and Masina, 2009; Vogt et al., 2010) marine biogeochemical models. In these models, Chl-a dynamics can be also defined with various levels of complexity. The ratio of Chl: C can be constant (non-photo-acclimation), or prognostic with a dynamic behaviour (Flynn and Flynn, 1998; Geider et al., 1998; Ross and Geider, 2009) or diagnostic with an empirical (Bernard, 2011) and mechanistic relationships (Geider and Platt, 1986; Doney et al., 1996; Bissett et al., 1999). These approaches vary in terms of their level of complexity, assumptions, and underlying mechanisms employed. The constant Chl:C ratio method is the simplest, assuming no acclimation or variation in Chl-a concentration. Prognostic techniques enable dynamic modifications of the Chl-a:C ratio, capturing phytoplankton’s adaptation to changing environmental conditions. Both empirical and mechanistic diagnostic techniques estimate Chl-a content based on observational or mechanistic correlations, respectively, without explicitly modelling the physiological processes. Each method provides information about the regulation of Chl-a dynamics and their ecological implications in marine ecosystems.

Other complexities have been added to the simple NPZD model, including 1) dissolved organic nitrogen in order to allow for a comprehensive representation of a significant nitrogen pool, particularly in regions where high primary productivity is observed (Hood et al., 2001; Druon et al., 2010), 2) oxygen cycle which is known to be an essential component to represent the oxygen minimum zone caused by excessive denitrification (Naqvi et al., 2000; Codispoti et al., 2001).

Unlike the simple NPZD, more complex models include different phytoplankton functional types (PFTs) and at least 15 state variables (e.g., Gregg, 2000; Gregg et al., 2003; Moore et al., 2004; Le Quéré et al., 2005; Dutkiewicz et al., 2009; Stock et al., 2014b; Wright et al., 2021; Luo et al., 2022). Considering a range of plankton groups and increasing the number of variables allows for a more thorough assessment of ecological interactions. Table 2 provides a summary of the models presented in this section, emphasizing their key characteristics and improvements. Here, we expand on the interpretation of each model improvement and discuss the relevance of increasing the number of variables.

The incorporation of Plankton Functional Types (PFTs) in biogeochemical models has been pursued to improve the accuracy of marine ecosystem representations (Anderson, 2005). Notably, the integration of additional complexity beyond NPZD models presents challenges, including limited understanding of ecology, data constraints, and the need to consolidate diversity within functional groups (Anderson, 2005). This approach is exemplified in biogeochemical models such as ERSEM (European Regional Seas Ecosystem Model) versions 1 and 2 (Baretta et al., 1995; Baretta-Bekker et al., 1997; Blackford et al., 2004). ERSEM adopts a functional role-based classification of biotic groups, categorizing phytoplankton as producers, bacteria as decomposers, and zooplankton as consumers, employing a general lower trophic approach. This approach is similar to the functional roles in an NPZD model, where each component has specific functions within the nutrient dynamics and planktonic community interactions. However, what distinguishes ERSEM from other NPZD models is its comprehensive integration of ecosystem-specific processes and interactions, allowing for a more tailored representation of biotic groups and their functions. The division of the food web into three functional categories is based on traits, including silica uptake, maximum growth rate, and nutrient affinity. The ERSEM model was initially applied to the North Sea to study the seasonal cycling of nutrients (i.e., nitrogen, phosphorous, silica, carbon). ERSEM was later modified to produce the Biogeochemical Flux Model (BFM), which accounts for the Chemical Functional Families (CFFs) in the state variables. The CFFs is split into living, non-living, and inorganic states (Vichi et al., 2007).

The Pelagic Interactions Scheme for Carbon and Ecosystem Studies (PISCES) is another example of a simple PFT model adapted from HAMOCC (Aumont et al., 2003) and extensively and widely applied in more than a hundred studies directly or indirectly (Aumont et al., 2015) The NASA Ocean Biogeochemical Model (NOBM) is a PFT model coupled with the Ocean-Atmosphere Spectral Irradiance Model (OASIM) (Gregg, 2001; Gregg et al., 2009) and is capable of simulating global-scale biogeochemical cycles. Additionally, the PlankTOM model describes the lower trophic level of marine ecosystems and has various extensions with different numbers of PFTs (Le Quéré et al., 2005; Kwiatkowski et al., 2014).

The Model of Ecosystem Dynamics, nutrient Utilization, Sequestration and Acidification (MEDUSA) includes the biogeochemical cycles of iron, silicon, and nitrogen as well as small and large plankton size classes (Yool et al., 2011; 2013). MEDUSA focuses on the biological sequestration of carbon in the deep ocean. Similarly, the Tracers of Phytoplankton with Allometric Zooplankton (TOPAZ), developed by Dunne et al., 2010, focuses on simulating the interactions among phytoplankton and allometric zooplankton, emphasizing their roles within the ecosystem dynamics.This model has been implemented in several studies (Jung, et al., 2019a; Bronselaer et al., 2020; Sharada et al., 2020) and extended in the Carbon, Ocean Biogeochemistry and Lower Trophics (COBALT) model, which improves the planktonic food web dynamics resolution to examine the influence of climate on the flow of energy from phytoplankton to fish (Stock et al., 2014a). Since the planktonic food web representation in TOPAZ is highly idealized due to an implicit representation of zooplankton and bacteria, COBALT addressed these limitations by including three explicit zooplankton groups and bacteria with a mechanistic parametrization of the impacts of these groups on the biogeochemistry (Stock et al., 2014b).

The DARWIN biogeochemical model is a more complex PFT model allowing for several different configurations and number of PFTs. It is initially developed to study Prochlorococcus species (Follows et al., 2007b) and later applied to global distributions of phytoplankton (Dutkiewicz et al., 2015; Lo et al., 2019). Nevertheless, the Regulated Ecosystem Model (REcoM) is based on a functional group approach (Schourup-Kristensen et al., 2014) and has also been commonly used in the Southern Ocean studies (Taylor et al., 2013; Losch et al., 2014; Hauck and Völker, 2015).

The complexity of the aforementioned PFT models depends on the number of the independent elements, processes parameterization (e.g., shape of functional responses, number of parameters, linear vs. nonlinear) along with the number of PFTs considered (Ford et al., 2018). As the number of PFTs increases, the complexity of the model increases as well. Generally, adding complexity to the model does not necessarily improve the model’s performance as the addition of complexity will require validation of each component of the ecosystem. Even under ideal circumstances, biogeochemical observations are lacking for full spatiotemporal coverage. As a result, the validity of additional processes is sometimes simply determined by how much they contribute to certain quantities like the total chlorophyll content (Ford et al., 2018). These studies comparing models of varying complexity have typically concluded that increasing model complexity does not always result in an increase in model skill (Friedrichs et al., 2007; Kriest et al., 2010; Ward et al., 2013; Kwiatkowski et al., 2014; Xiao and Friedrichs, 2014).

Bacteria play a key role in their contribution in the microbial loop, particularly in the recycling of dissolved organic matter (DOM) and provide a linkage between the higher trophic level and the dynamics of DOM (Pauer et al., 2014). In biogeochemical models studying dissolved organic matter (DOM) dynamics, various classes of organic matter, along with chemical elements and Plankton Functional Types (PFTs), have been incorporated. DOM is divided into three distinct pools based on its lability: labile DOM, semi-labile DOM, and refractory DOM. Labile DOM consists of organic molecules with a lifespan ranging from minutes to days. Semi-labile DOM undergoes decomposition over a period of 1–20 years and typically accumulates above a seasonal or permanent pycnocline. Refractory DOM, the most abundant component, exhibits the longest residence time in the water column, persisting for tens of thousands of years (Carlson and Hansell, 2015; Osterholz et al., 2015).

The ratios of DOC:DON vary among different DOM pools. The refractory DOM pool has a ratio of 3511:202, whereas the labile DOM pool has a ratio of 199:20 (Hopkinson and Vallino, 2005). These ratios indicate the relative amounts of dissolved organic carbon to dissolved organic nitrogen within each pool. It is essential to consider these ratios when incorporating DOM dynamics into biogeochemical models.

DOM in the oceans is an important input to the biogeochemical cycle, with the ocean’s DOC pool serving as the biosphere’s largest storage facility for reduced organic carbon (Hansell et al., 2009). Labile and, to a lesser extent, semi-refractory pools are driven by microbial reactions. In the microbial loop, the DOC-bacteria pathway constitutes a primary mechanism through which primary production is eitherrespired to CO₂, directed to higher trophic levels, or transformed into refractory compounds. However, the refractory DOC is transported deep into the ocean (Anderson, 1999; Kahler and Koeve, 2001).

The aforementioned different forms of organic matter have been incorporated into some biogeochemical models. For instance, the semi-labile DOC and semi-labile dissolved organic nitrogen (DON) have been included in biogeochemical models that focus on regions known of their high primary productivity, particularly coastal and continental shelves regions (Druon et al., 2010; Fennel, 2010; Xue et al., 2013). These areas are potential sites for sediment denitrification and organic carbon transport to the open ocean (Devol and Christensen, 1993; Rao et al., 2007; Fennel et al., 2008; Druon et al., 2010). Another study by Schartau et al. (2007) considered the labile DOM exudation and extracellular POC formation. This process occurs during the growth of phytoplankton in which the assimilation of dissolved inorganic carbon (DIC) by phytoplankton is exuded as DOC, which is simultaneously converted into extracellular POC. A significant proportion of extracellular POC is available as transparent exopolymer particles (TEP), which consist mainly of dissolved polysaccharides (PCHO). PCHO exudation is associated with high consumption of DIC that is not directly measurable from the use of dissolved inorganic nitrogen (DIN), a phenomenon known as carbon overconsumption (Toggweiler, 1993).

Sediment processes have been also integrated into the pelagic biogeochemical models to investigate the impact of atmospheric conditions on the ecosystem dynamics, the distribution of sediments processes in the water column, and the nitrogen budget derivation (Soetaert et al., 2001; 2000; Arthur et al., 2016; Ehrnsten et al., 2019; Kearney et al., 2020). The lack of sediment regeneration of nutrients in the biogeochemical model may lead to underestimating primary productivity (Liu et al., 2007). Adding a representation of marine sediment in biogeochemical models is important, particularly in shallow water regions like shelves, where strong vertical mixing occurs. This inclusion helps to accurately account for nutrient recycling from sediments to surface waters, preventing the underestimation of productivity (Yool et al., 2013). However, it is crucial to emphasize that in many models, sediment is represented as a boundary condition with simplistic representation, and even global models that address diagenesis frequently lack a thorough representation of the benthic ecosystem.

While incorporating sediment interactions in biogeochemical models is important to understand the linkages between sediment diagenetic processes and the global carbon cycle (Middelburg, 2018; Snelgrove et al., 2018), many models use an oversimplified representation of sediment. Several regional models have been developed to examine specific areas of the ocean at different timescales, yielding important insights into biogeochemical processes. These models such as (Soetaert et al., 2000) and the Ecological ReGional Ocean Model with vertically resolved sediments (ERGOM SED 1.0) (Radtke et al., 2019) also add to the discussion of sediment representation and its effects on oceanic and coastal ecosystems. For example, in a coastal environment (Reed et al., 2011), investigated the dynamics of sedimentary phosphorus and the development of bottom-water hypoxia. Furthermore, the relationship between sediment dynamics and biogeochemical cycling was examined by (Radtke et al., 2019).

Marine sediments representation in the biogeochemical model can be parametrized using a meta-model that is based on the downward fluxes of organic matter (Middelburg et al., 1996), a sediment box module which uses simple 1D or 2D submodules in the pelagic model to represent a simple interaction of the water column with the sedimentary pools (Yool et al., 2013), or multilayer sediment modules with chemically active sediment layers (Heinze and Maier-Reimer, 1999). Most of the biogeochemical models utilize sediment box modules or meta-models to examine eutrophication (Cerco et al., 2006; Soetaert and Middelburg, 2009; Fennel et al., 2011) or hypoxia in several regions, such as Tokyo Bay (Sohma et al., 2008), the Baltic Sea (Reed et al., 2011) and the North Sea (Lancelot et al., 2005; Meire et al., 2013).

The aforementioned biogeochemical models have been applied in various tropical seas due to their high dynamism and significant biogeochemical activity, making them among the most active zones in the ocean (Lønborg et al., 2021). This section encompasses multiple regions, including the tropical Indian Ocean (specifically the Arabian Sea), the South China Sea, the tropical Pacific Ocean, and the oxygen minimum zone in the Arabian Sea. One prime example of a tropical ecosystem is the Indian Ocean, which experiences the world’s most significant monsoon system caused by the seasonal warming and cooling of massive air masses (Wiggert et al., 2005). In winter and summer, the wind patterns near the ocean change, with the southwest monsoon (SWM) providing half of the year’s wind and the northeast monsoon (NEM) providing the other half (Wiggert et al., 2005). The dynamics of the Arabian Sea’s (the most extensively studied ecosystem in the Indian Ocean) phytoplankton bloom clearly show a temporal alignment with the southwest monsoon and northeast monsoon (Madhupratap et al., 1996; Wiggert et al., 2005). Strong south-westerly winds during the summer propel coastal upwelling and offshore movement of entrained nutrients, which promote the Arabian Sea’s observed phytoplankton blooms (Brock and Mcclain, 1992; Hitchcock et al., 2000; Lierheimer and Banse, 2002). Modest north-easterly winds during the winter carry cool, dry air from the Tibetan Plateau over the northern Arabian Sea (Huang, 2010). Additionally, convective mixing during NEM entrains nutrients and promotes biological activity north of 15°N (Madhupratap et al., 1996). Thus, primary productivity and nutrients budgets in the Arabian Sea, spanning a temporal scale of seasonal variations and a spatial scale ranging from 5°S to 27°N are highly impacted by the seasonal cycle and mesoscale impacts (Resplandy et al., 2011). The primary production in terms of surface Chl-a production in the Arabian Sea is widely studied using several types of ocean biogeochemical models. For example, the estimated surface Chl-a by PISCES coupled to the Nucleus for European Modelling of the Ocean (NEMO) is compared with the weekly climatology of SeaWiFS (Sea-viewingWide Field-of-view Sensor Data). The PISCES-NEMO reproduced surface Chl-a levels of more than 0.2 mg m^-3 during the NEM. It also predicted the maximum intensity of surface Chl-a, which exceeds 0.5 mg m^-3, but underpredicted surface Chl-a by 30%–50% north of 20°N. This bias was found in other studies and may be attributed to the lack of diurnal layer mixing (Kawamiya and Oschlies, 2003; Wiggert et al., 2006; Koné et al., 2009). The bias in the bloom dynamics is explained by the lack of deepening of the mixed layer and nutrient feed (Wiggert et al., 2006), and light limitation (Koné et al., 2009). Another cause of inconsistencies in the PISCES model (Resplandy et al., 2011) is partially linked to the underprediction of silicate input from the Sea of Oman. PISCES captured the observed high surface Chl-a along the coasts of Somalia and southwestern India, reaching values higher than 0.4 mg m^-3 during the SWM. On the other hand, it underestimated Chl-a in Oman. These inconsistencies can be attributed to the overprediction of the offshore advection by the Great Whirl and the Socotra Passage, resulting in an excessive outflow of warmer and oligotrophic waters from the Aden Gulf. Therefore, PISCES-NEMO showed better agreement with observations (seasonal SeaWiFS data) for the surface Chl-a distribution than in (Koné et al., 2009) with a correlation coefficient (r) around 0.6. These comparisons encompass a spatial scale ranging from 30°S to 30°N and from 30°E to 120°E, with a focus on seasonal variations. In addition to PISCES, the PFT model NOBM coupled with the oceanic general circulation and radiative forcing model was able to capture the intra-seasonal variability of the Indian summer monsoon, specifically focusing on the spatial scale of the Bay of Bengal (BOB) (6°–22° N; 80°–100° E) and Arabian Sea (AS) (3°–17° N; 55°–73.5° E) in the North Indian Ocean. The study showed important differences between the two regions indicating that the Arabian Sea is richer in nutrients than the Bay of Bengal with a coefficient of determination (R²) around 0.78 for the monthly mean surface Chl-a (Das et al., 2019). Furthermore, surface Chl-a is reproduced by Fasham-ROMS (Regional Ocean Modeling System), in the study region, which extends in latitude from 5°S to 30°N and in longitude from 34°E to 78°E. The performance of the model is evaluated by comparing its outputs with seasonal SeaWiFs observations, resulting in correlation coefficient ranging from 0.69 in summer to 0.74 in winter (Lachkar et al., 2017).

In addition to surface Chl-a estimation, surface nitrogen content in the Arabian Sea is estimated by the European Regional Seas Ecosystem Model (ERSEM) coupled with a 1D Princeton/Mellor–Yamada (Blackford and Burkill (2002). The model performance was assessed at three stations located at 19°N 59°E, 16°N 62°E, and 8°N 67°E, covering a seasonal timescale. The coupled model has a bias for surface nitrate concentration (3 mmol m^-3) and lacks horizontal advective processes, which might be enhanced by including diurnal physical processes (horizontal advection, vertical mixing, or tidal currents). Another study by using Fasham-ROMS indicated a correlation coefficient (r) of 0.88 with observation for surface nitrate in summer. Nitrogen is not the only limiting nutrient affecting phytoplankton growth in the Arabian Sea and the ignorance of other components such as iron may underestimate nitrogen budget in the region. Iron, as an important limiting nutrient in the Arabian Sea, plays a significant role in facilitating nitrogen fixation, which is a crucial process driving biological productivity. However, nitrogen fixation is not taken into account in the NPZD, which may amplify the impact of denitrification on the nitrogen budgets in the Arabian Sea (Lachkar et al., 2017). A few studies have been conducted to examine the effect of iron limitation on global biogeochemical processes, including those in the Indian Ocean. For instance, studies using PFT approach, such as PISCES-ROMS and TOPAZ-MOM (Jung et al., 2019b) have explored the implications of iron limitation on various aspects, including macronutrient availability and ecosystem dynamics. Unlike the previous version of PISCES in (Aumont and Bopp, 2006), this version incorporated new parametrizations for iron chemistry and nitrogen fixation. PISCES exhibited a better iron estimate than TOPAZ, whereas TOPAZ showed larger bias with observation (SeaWiFS) (1.5 nMol m^-3) (Jung et al., 2019b).

Another prominent feature of the Arabian Sea is its oxygen minimum zone which is noted for its low oxygen concentrations (<4 mmol/m^-3) at depths of 150–1,250 m in the central part of the Arabian Sea. Thus, the dynamics of the oxygen minimum zone (OMZ) at spatial coordinates (65°E, 13°N) is studied by applying ERSEM coupled with General Ocean Turbulence Model (GOTM). The model incorporates most of the biogeochemical processes driving the dynamics of OMZ with the inclusion of denitrification process which was not considered in Butenschön et al. (2016). ERSEM-GOTM ignores the episodic intrusion of oxygen within the OMZ which may play a role in maintaining aerobic microorganisms. Therefore, the aerobic biomass within the OMZ is maintained by the residual initial conditions (Sankar et al., 2018). In addition, the model is unable to reproduce meso-zooplankton at depth, and the discrepancy between the estimates and the climatological seasonal cycles of nutrients and oxygen at depth might be due to the absence of lateral circulation in the model. Further studies on the OMZs in the Arabian Sea were investigated by Lachkar using a simple NPZD approach with the addition of the oxygen cycle (Fasham- ROMS) (Lachkar et al., 2017; Lachkar et al., 2019; Lachkar et al., 2020). Fasham-ROMS indicated that the primary productivity and OMZ are highly impacted by monsoon wind intensification, with an overall high correlation coefficient (r) 0.93 for oxygen (both seasons) in the upper layer. While the model’s horizontal grid resolution of 1/12° and vertical grid consisting of 32 sigma layers allow for some representation of eddy dynamics, the resolution may still be too coarse to fully capture all eddy processes.The OMZ intensity (oxygen concentration) has also been overestimated because of the underestimation of the northern Arabian Sea ventilation. Studies on OMZ in the Arabian Gulf are still limited.

In addition to the Indian Ocean, some parts of the Pacific Ocean are also located in the tropical zone, and south China Sea, which is the largest marginal seas of the West Pacific Ocean (Qi-zhou et al., 1994), is an example of an ecosystem that belongs to this region. Similar to the Arabian Sea, this region is largely controlled by seasonal monsoon circulation. Strong northeast (NE) winter monsoon and weaker southeast (SE) summer monsoon predominate at the sea surface of the South China Sea (SCS) (Qi-zhou et al., 1994; Li et al., 2021). Studies have shown a clear response of Chl-a and temperature of the upper oceanic conditions to atmospheric forcing in the SCS (Sun et al., 2009; Sun et al., 2010; Yang et al., 2010; Zhao et al., 2014). Estimates of the mean surface Chl-a has been reported with varying performances. A correlation coefficient (r) of 0.85 for the annual mean surface Chl-a is obtained by (Jung et al., 2019b) using both global PISCES and TOPAZ coupled to NEMO representing surface Chl-a concentration conditions (i.e., phytoplankton blooms) in the tropical Pacific. Furthermore, the monthly mean surface Chl-a was determined by PISCES-NEMO and TOPAZ-NEMO (Jung et al., 2019b) with r equal to 0.6 and 0.75, respectively. However, both models overpredicted the zonal averaged Chl-a in the Pacific Ocean east of Japan by 67%, primarily because of errors in the Kuroshio current path, as seen in low-resolution models (Gnanadesikan et al., 2002; Dutkiewicz et al., 2005; Aumont and Bopp, 2006). Nevertheless, the regional Chai’s NPZD model (Chai et al., 1996) configured for the central South China Sea was modified by adding phosphate, silicate, dissolved oxygen, DIC, two phytoplankton groups (diatoms, small phytoplankton) and their chlorophyll components (Chl1 and Chl2) has shown a correlation coefficient of approximately 0.62 for the monthly mean surface Chl-a in the SCS compared to SeaWiFS data (Ma et al., 2019). Other NPZD models were also applied in the tropical Pacific seas (SCS) to study the physical processes of mesoscale eddies’ impacts on biology such as in (Xiu and Chai, 2011). The study covered the entire Pacific Ocean, spanning from 45°S to 65°N and from 99°E to 70°W. The model exhibited a bias of approximately 0.05 mg m^-3 (April) and ∼0.1 mg m^-3 (September) for mean surface Chl-a concentrations when compared to SeaWiFS data for the respective months.

Other biogeochemical components such as nutrients including nitrate, phosphate, silicate, and iron have been estimated using global ocean biogeochemical models in the Pacific region. Higher positive bias for nutrients is obtained using the global TOPAZ-NEMO model compared to PISCES-NEMO in the Southern Pacific Ocean using the same physical model setup. With bias of 4.5 μmol kg⁻¹,0.32 μmol kg⁻¹, and 16 μmol kg⁻¹ for nitrate, phosphate, and silicate respectively, with a higher positive bias in the central and southern Pacific for silicate (8 μmol kg⁻¹). With regards to a NPZD based approach, Fasham-ROMS exhibited the impact of physical processes, such as lateral and vertical fluxes, as well as vertical diffusion in the shelf water, on nutrient distribution in the SCS (Ji et al., 2017; Lu et al., 2020) showing a reasonable linearity between the model estimates and field observations (JGOFS: a global, multidisciplinary program with participants from over 20 countries) for surface nitrate with R²: 0.95 (Ji et al., 2017). In (Ji et al., 2017), the model domain focused on the northwestern Pacific (NWP) region (3°–52°N, 98°–158°E), while in (Lu et al., 2020), the model domain extended from approximately 0.95°N, 99°E in the southwest corner to around 50°N, 147°E in the northeast corner.

With regards to iron, low surface concentrations of iron were predicted in the North and North Central Pacific, while high concentrations were predicted in the North and Equatorial Indian Ocean using global NOBM (Gregg et al., 2003). The model domain considered spans from approximately −84°S to 72°N latitude in increments of 1.25° longitude by 2/3° latitude. The overall correlation coefficient for the annual mean iron with in situ data is 0.86 and the reason for the iron overprediction is attributed to the lack of scavenging, excessive remineralization, and slow detritus sinking rate.

Compared to the tropical zones, the temperate zone is characterized by moderate climate with mean annual temperatures. Average monthly air temperatures in the temperate zone typically exceed 10°C during the warmest months while water temperatures generally remain above −3°C, with a minimum around ∼-2°C. (Trewartha and Horn, 1980). In contrast, the subtropical zones, located between the tropical and temperate zones, are characterized by hot summers and moderate winters. Some parts of the North Atlantic Ocean comprises both temperate and varying tropical to subtropical zones. Biogeochemical studies in this region have focused on understanding primary productivity in continental shelves and coastal areas, which are highly productive ecosystems and are significantly influenced by human activities (Economidou, 1982; Fennel, 2010). Additionally, air-sea CO₂ fluxes in coastal waters exhibit greater variability and larger amplitudes compared to fluxes in the open ocean (Tsunogai et al., 1999; Cai et al., 2003; Thomas et al., 2004). For instance, in the Mid Atlantic Bight, located in the central region of the eastern U.S. continental shelf, high primary productivity rates, sediment denitrification, and strong residual circulation contribute to the region’s significance as a potential source of organic carbon export to the open ocean (Druon et al., 2010). In the Northwest North Atlantic and other ecosystems within the region, many studies have utilized simple NPZD models to investigate the influence of physical processes on biological dynamics, rather than aiming to capture the full complexity of ecological interactions. Therefore, this section covers both coastal areas and open sea regions of the Mid Atlantic Bight, northeastern US continental shelves, Georges Bank, Gulf of Mexico, and Mediterranean waters. An example of such modeling approach is demonstrated in (Fennel et al., 2008), where the distribution of Chl-a, an indicator of productivity in the Atlantic regions, is predicted using a NPZD model. The model is coupled with a high-resolution physical model ROMS, which incorporates accurate descriptions of physical processes and an explicit representation of inorganic carbon chemistry and air-sea gas exchange of CO₂ (Fennel et al., 2008). The model’s inclusion of tides in ROMS improves the seasonal mean surface Chl-a representations, showing good agreement with SeaWiFS climatology in Georges Bank (correlation coefficients: 0.75 in winter, 0.72 in spring, 0.85 in summer, and 0.83 in fall) (Fennel et al., 2008). In another study by Druon et al., 2010, Fennel’s model (Fennel et al., 2008) was further enhanced by incorporating a dissolved organic matter (DOM) module to investigate organic carbon dynamics in the North-eastern U.S. continental shelves. The model demonstrated good agreement with monthly satellite data, particularly in the inner shelf and on Georges Bank, attributing this agreement to tidal mixing and continuous nutrient suppl (Bias: Chl-a (with DOM) = 4 mg m^-3, Chl-a (without DOM) = 6 mg m^-3). The inclusion of the DOM module improved the model by capturing the role of dissolved organic matter in nutrient cycling and its influence on primary production in the North-eastern U.S. continental shelf. In a separate study by (Xue et al., 2013), phosphate cycling was omitted from Fennel’s model (Fennel et al., 2011) because previous studies (Rabalais et al., 2002) showed that the primary production is generally nitrogen limited in the Louisiana-Texas shelf in the Gulf of Mexico. Therefore, the correlation coefficient of the seasonal mean of simulated and observed (MODIS: Moderate Resolution Imaging Spectroradiometer) surface Chl-a for the spring, summer, fall and winter were 0.6, 0.65, 0.53, 0.45 respectively. The simulated data successfully captured the high concentration of surface Chl-a in coastal areas close to main rivers, such as the Bay of Campeche, Campeche Bank, and the Louisiana-Texas shelf, whereas Chl-a concentrations were very low in the deep ocean. It is important to note that while Louisiana/Texas shelf in the Gulf of Mexico, as discussed in Xue et al. (2013), falls outside the tropics, the Gulf of Mexico circulation does include tropical latitudes. This reference to the Gulf of Mexico was mentioned to highlight the omission of phosphate cycling in Fennel’s model due to nitrogen limitation in the primary production observed in the Louisiana-Texas shelf (Rabalais et al., 2002). In contrast to continental shelves and coastal areas of the Mid Atlantic Bight, some areas of the Atlantic Ocean and the Mediterranean Sea shows limitation of phosphate rather than nitrogen. The Mediterranean Sea encompasses the most oligotrophic ecosystems in the global ocean, with low nutrient and low chlorophyll content (Krom et al., 1991; Durrieu et al., 2011; Thingstad et al., 2014), being characterised by strong seasonality as well as highly dynamic system in terms of physics and biology (Ortenzio et al., 2009; Christaki et al., 2011). This limitation of phosphate was modelled using ERSEM which accounts for variable stoichiometric formulation in phytoplankton cells allowing the adaptation of internal N:P quota to different limiting environments (Lazzari et al., 2016). ERSEM reproduced the seasonal cycle and the observed patterns of surface Chl-a and inorganic nutrients in the Mediterranean Sea (MS). The interannual variability of Chl-a showed an agreement of r: 0.75 with SeaWiFS observation for the MS. However, the model underestimates phytoplankton blooms in the western basin of the MS while overestimates the blooms in the eastern basin during spring chlorophyll peaks. This discrepancy is attributed to the overestimation/underestimation of the winter mixing intensity by the physical model (Kalaroni et al., 2019). The spatial variability of the simulated surface Chl-a was assessed by calculating the correlation coefficient (r) between the model outputs and SeaWiFS observations. The correlation coefficient (r) was found to be 0.64, indicating a moderate positive linear association between the two datasets. The standard deviation of the spatial model variability was computed to be lower than 1 (normalised), suggesting relatively low dispersion of Chl-a values across the Mediterranean Sea. The temporal variability, on the other hand, was evaluated using unbiased root mean square differences (RMSDs) for the annual mean surface spatial and temporal variability. The RMSD value for the annual mean surface spatial variability was 0.78, indicating the average magnitude of differences between the model simulations and the observed data. Likewise, the RMSD value for the annual mean surface temporal variability was 0.68, suggesting the average magnitude of differences between the model simulations and the observed data over time. Lower RMSD values indicate better agreement in terms of both spatial and temporal variations.

The effects of climate change, such as a considerable decrease in ice extent and rising sea levels, can exacerbate stratification in polar and subpolar regions (Farmer et al., 2021). In addition, changes to the global climate system are anticipated to have a significant impact on polar and subpolar marine ecosystems through their effects on water temperature and sea ice cover hence affecting marine production and carbon cycle (Gibson et al., 2020). This section covers the Arctic Ocean, Antarctic Peninsula, Drake Passage, Scotia Sea, northern Weddell Sea, South Georgia, Barents Sea, Hudson Bay, Beaufort Gyre, Labrador Sea, Baffin Bay, Canadian Arctic Archipelago, Russian shelves (East Siberian Sea), and the Dotson Ice Shelf in West Antarctica, along with the adjacent Amundsen Sea Polynya. Several ocean biogeochemical studies in the polar region focused on the effect of climate change on the carbon cycle (MCGUIRE et al., 2009; Terhaar et al., 2021; 2019; 2020; Negrete-garcía et al., 2019; Bruhwiler et al., 2021) and nutrients dynamics (Hátún et al., 2017; Randelhoff et al., 2018; 2019; Tuerena et al., 2021). Therefore, carbon is estimated using biogeochemical models with the inclusion of an explicit representation of the DIC component, such as PISCES-MIT general circulation model (MITgcm), which have been used to qualitatively examine the carbon cycle in the Arctic Ocean and have demonstrated their capacity to capture observed seasonal and regional changes in dissolved pCO₂ (Manizza et al., 2011). PISCES-MITgcm represented summer levels of surface pCO₂ in the Canadian Archipelago (200–250 µatm) correctly, but underestimated pCO₂ spring values (∼300 µatm) relative to observations (Carbon Dioxide in the Atlantic (CARINA) data set (Jutterstrom et al., 2010)) (400–450 µatm). The lack of riverine POC and the impact of terrestrial carbon resulting from coastal erosion in the model resulted in the underprediction of the observed carbon levels. Additionally, sedimentation and resuspension processes, which may be important to enrich the water column with carbon hence impacting air-sea gas exchange, were neglected in the model. When comparing PISCES-NEMO and TOPAZ-NEMO (Jung et al., 2019b) for annual mean (averaged over all depths) DIC, TOPAZ-NEMO’s results agreed with observation (r > 0.95), and the zonally averaged surface content was better represented by TOPAZ-NEMO in the Northern Hemisphere (bias < 10 μmol kg⁻¹). However, surface mean DIC reproduced by TOPAZ and PISCES showed a bias of −90 and 90 respectively with GLobal Ocean Data Analysis Project version 2 (GLODAPv2 (Lauvset et al., 2016)).

In addition, a modified version of Chai’s model (Chai et al., 1996) coupled to ROMS is applied over a model domain that spans the Antarctic Peninsula, Drake Passage, Scotia Sea, northern Weddell Sea, and South Georgia, indicating that dominant iron sources in the Scotia Sea are derived from sediments in the Antarctic Peninsula shelf along with the South Orkney Plateau (Jiang et al., 2019). In addition to these sources, the Antarctic Circumpolar Current, the northern side of the Weddell Gyre, upwelling, atmospheric dust deposition, and icebergs are the common sources of iron in the Southern Ocean (Jiang et al., 2019). However, the monthly means of iron levels estimated by the modified Chai’s model have shown an average overestimation by 0.26 nM, deviating from the observed average value of 0.35 nM, resulting in r: 0.76.

On the other hand, the representation of the monthly observed nitrate, mixed layer depth, as well as euphotic layer depth is reproduced in the Arctic region using a skill assessment of 21 regional and global coupled biogeochemical models based on a phytoplankton functional types, including PISCES, PlankTOM, COBALT, TOPAZ, HAMOCC, BIOMASS, MEDUSA, ERSEM, PELAGOS, PISCES and NOBM (Babin et al., 2016). Most of these models have shown a positive bias for depth-averaged nitrates, explaining the overestimation of nitrates in the upper layer (r ≤ 0.68), and none of these models were able to represent accurately the variability in field measurements (NOAA/NODC ArcNut database (Codispoti et al., 2013)). However, REcoM with its variable stoichiometry for phytoplankton’s growth and nutrient cycling (Geider et al., 1998), coupled with a high-resolution physical model: Finite Element Sea Ice-Ocean Model (FESOM), exhibited a correlation coefficient of (r = 0.75) for spatial seasonal surface values of both dissolved inorganic nitrogen (DIN) and silicate (Schourup-Kristensen et al., 2018). The simulated DIN concentration indicated low levels in the Russian shelves (East Siberian Sea) with the exception of the Lena Delta in which the nutrients are added to the water by the riverine supply.

As for the representation of Chl-a in polar ecosystems, the original version of BLING coupled with NEMO was used to study monthly high surface Chl-a levels (i.e., blooms) in the high-latitude eastern North Atlantic and the central Arctic (Castro de la Guardia et al., 2019). This model is compared with MODIS observations and showed R²: 0.76, 0.42, 0.81, 0.9, 0.1, and 0.93 in the Barents Sea (BB), Hudson Bay (HB), Beaufort Gyre (BG), Labrador Sea (LS), Baffin Bay (BB), and Canadian Arctic Archipelago (CAA) respectively. However, it underestimated the spring bloom by 1.7 mg m^-3 and the fall bloom by 0.7 mg m^-3 in the BB, likely because it did not take into account nutrients from riverine input. Furthermore, the model exhibited a tendency to overpredict the spring blooms in Hudson Bay (HB) and Baffin Bay (BB) by approximately 0.5 mg m^-3 and 0.3 mg m^-3 respectively, occurring about a month earlier than expected (March instead of April). This discrepancy might be attributed to the underprediction of sea ice concentration. Additionally, the REcoM model exhibited a bias of 1.1 mg m^-3 in predicting widespread subsurface chlorophyll maxima in the central Arctic Ocean, as compared to the mean summer surface Chl-a data obtained from Globcolor observations (Schourup-Kristensen et al., 2018).

Moreover, the MITgcm-REcoM model is used to investigate the impact of sea ice on phytoplankton blooms and surface dissolved inorganic carbon (DIC) in the ocean coastal-shelf-slope ecosystem west of the Antarctic Peninsula (Schultz et al., 2019). The study area encompasses a vast region ranging from 74.7°S, 95°W in the southwest to 55°S, 55.6°W in the northeast. Simulation results were compared with ocean color observations, and research ship surveys from the Palmer Long-Term Ecological Research (LTER) program. The Palmer LTER data and MODIS/SeaWifs data showed significant discrepancies in the study of the interannual variability of the phytoplankton bloom, with correlations for surface Chl-a in different sub-regions ranging from 0.52 to 0.73. While late summer retreat years displayed positive or low negative anomalies, early sea-ice retreat years primarily displayed negative or neutral anomalies in Chl-a. The model demonstrated the influence of sea-ice dynamics on phytoplankton dynamics by indicating an earlier bloom and lower Chl-a concentrations (relative to climatology) in years with early retreat of sea ice, resulting in lower chlorophyll concentrations at the end of the summer. Conversely, in years with late retreat of sea ice, the model showed greater Chl-a concentrations later in the season. Despite discrepancies between the model and observations, the model captured the temporal shifts in bloom timing during years characterized by high or low ice extent. Additionally, the model reproduced overall seasonal anomalies in dissolved inorganic carbon (DIC), which are related to seasonal net community production (NCP).

DARWIN is another model that has been used to study regional phytoplankton distribution in the Southern Ocean (Lo et al., 2019). The model was customized to improve its accuracy by increasing the affinity of phytoplankton for nutrients, adding two distinct size classes of diatoms, and considering two different life stages for Phaeocystis (single-cell vs. colonial). These improvements increased the agreement between the simulated coccolithophore and diatom levels with the in-situ observations. However, the model inaccurately simulated monthly diatoms and haptophytes, including small non-silicified phytoplankton, in the Ross Sea. Additionally, the model overestimated small non-silicified phytoplankton in this region, resulting in a general mean absolute error for diatoms and haptophytes of 0.74 mg m⁻³ and 0.22 mg m⁻³, respectively. This inconsistency can be partially attributed to inaccuracies in representing PFT phenology and distribution, and the potential impact of iron limitation should also be considered. The representation of co-existing coccolithophores and Phaeocystis remains a challenge, and any small changes in DARWIN physiological parameters can lead to either Phaeocystis or coccolithophores loss. Furthermore, the explicit representation of sea-ice algae, which is important for accurately modeling ice-covered regions, has not been adequately incorporated.

In addition to studying Arctic and Antarctic phytoplankton dynamics, researchers have also been investigating the role of polynyas in polar ecosystems. Polynyas, which refer to expansive openings in the sea ice cover, play a crucial role in facilitating the exchange of heat, moisture, and gases between the atmosphere and the ocean (Maqueda et al., 2004; Dieckmann and Hellmer, 2009). Moreover, they drive patterns of ocean circulation, contribute to nutrient cycling, and influence ecosystem dynamics in polar regions (Maqueda et al., 2004). The rapid melting beneath the Dotson Ice Shelf (DIS) in West Antarctica, caused by warm, saline Circumpolar Deep Water (CDW) intrusions, has a significant impact on the adjacent Amundsen Sea Polynya (ASP) and its biology. Further, studies focusing on polynyas has garnered significant attention due to their substantial impact (Fichefet and Goosse, 1999; Hunke and Ackley, 2001; Wu et al., 2003). A recent study by (Twelves et al., 2020) using a coupled MITgcm-BLING model examined the effect of ice shelf melting on net primary production (NPP). The results show that melting ice shelves enhances upper ocean iron concentrations, which increases NPP. The study highlights how phytoplankton self-shading delays the bloom and reduces peak NPP while simultaneously emphasizing that changes in CDW intrusion regulate the horizontal dispersion of iron-rich meltwater, thereby impacting NPP levels.

Biogeochemical models applied in this zone are mostly based on PFT approach which account for at least two plankton species. The section covers the Subarctic northeast Pacific at Ocean Station Papa (OSP), HNLC areas of the Southern Ocean, and the Southern Ocean (Antarctic Circumpolar Current region). Surface Chl-a concentrations are low in this region caused by iron limitation (Edwards et al., 2004). More recently (Zhang et al., 2021), used an extended version of the Chai’s model (CoSiNE-Fe) that incorporated iron cycling to investigate phytoplankton dynamics at the Ocean Station Papa (OSP) in the subarctic northeast Pacific (50°N, 145°W). The model was able to demonstrate a strong correlation between nitrate and the averaged depth-integrated Chl-a in the upper 200 m of the water column, with correlation coefficients of 0.8 and 0.69, respectively, compared to observations from a Biogeochemical-Argo float (BGC-Argo). The model’s performance was also evaluated using root-mean-square differences (RMSDs), which were 0.65 and 0.98 for nitrate and averaged depth-integrated Chl-a (in the upper 200 m), respectively.

However, the global PFT-based model PlankTOM10 coupled to NEMO was used to evaluate the role of grazing versus iron limitation in HNLC areas of the Southern Ocean (Le Quéré et al., 2016). PlankTOM10 was able to reproduce annual mean surface Chl-a estimates with a correlation coefficient (r) around 0.8 in the summer season when macro-zooplankton grazing is more intense. While this suggests that microzooplankton grazing may be a significant factor contributing to the low Chl-a concentrations observed in HNLC areas, it is important to note that high growth of phytoplankton can coexist with grazing (Lu¨rling, 2021). While both PlankTOM10 and PlankTOM6 share similar formulations, the former encompasses a greater number of plankton groups. These models yielded comparable outcomes for the annual mean surface Chl-a concentrations, demonstrating a correlation coefficient (r) of approximately 0.8 when compared with SeaWiFS satellite data. PlankTOM10, in particular, displayed slightly better performance to PlankTOM6 in terms of seasonal surface Chl-a distribution, with a bias of only 1.2%.

In contrast, DARWIN-MITgcm (Dutkiewicz et al., 2015) captured significant spatial variability in the annual mean surface Chl-a across the global oceans. Despite this, the model exhibited a relatively lower correlation coefficient of around 0.55 and tended to overestimate surface Chl-a concentration, particularly in the Southern Ocean.

Surface nitrate concentration, on the other hand, is high relative to surface Chl-a production in the HNLC zone (Edwards et al., 2004) and the global PlankTOM10 model predicted slightly lower (by 5%) of annual mean surface nitrate concentrations than observations World Ocean Atlas 2009 (Garcia et al., 2006; Le Quéré et al., 2016). Biases in the PlankTOM10 model can indeed be attributed to various limitations, including simplified overwintering strategies for zooplankton, coarse representation of iron (Fe) dynamics, the absence of certain organic matter types such as semi-refractory organic matter, and the exclusion of ecosystem pathways like viral lysis. Additionally, the model’s omission of specific zooplankton groups like salps, pteropods, and auto- and mixotrophic dinoflagellates can also contribute to biases in its predictions.

Although iron is required at low concentrations for phytoplankton growth, the lack of its sources in the HNLC zone particularly the Southern Ocean highly impact the primary productivity and phytoplankton composition. For example, in the simplified version of DARWIN biogeochemical model (Dutkiewicz et al., 2009), which consists of two species (diatoms and small phytoplankton), coupled with MITgcm model, suggested that iron supply to the surface layer of the open Southern Ocean (Antarctic Circumpolar Current region) is highly driven by eddies (Uchida et al., 2019). Eddies can improve vertical mixing, influencing the upwelling of iron from deeper waters. This transport of iron from subsurface reservoirs to the sunlit layer where phytoplankton thrive can stimulate their growth and productivity. Therefore, to improve the accuracy of iron fluxes in the Southern Ocean, it is important to consider the role of eddies in addition to the atmospheric source. In particular, the inverse energy cascade caused by meso- and sub-mesoscale baroclinic instabilities plays a crucial role in determining iron distributions (Person et al., 2019).

Additionally, a recent study by Fu (2022) which used an inverse biogeochemical model with phosphorus, carbon, and oxygen Modules (Wang et al., 2019) to investigate the effects of ocean iron fertilization (OIF) in high nutrient low chlorophyll (HNLC) regions, specifically the Southern Ocean, North Pacific, and eastern equatorial Pacific. The findings revealed that OIF led to increased productivity, causing a downward transfer of nutrients and carbon to deeper waters. This process was accompanied by improved efficiency of the soft tissue pump and an increased in the CO₂ uptake. However, these positive effects were counterbalanced by a decrease in global mean oxygen levels and an expansion of oxygen minimum zones, suggesting potential consequences for ocean oxygenation and the development of hypoxic conditions. Notably, the study used (Wang et al., 2019) biogeochemical model and exhibited strong agreement with observations (GLODAPv2), showing correlation coefficients of 0.92 for dissolved inorganic phosphorus (DIP), 0.93 for dissolved inorganic carbon (DIC), and 0.85 for oxygen (O₂). Interestingly, a comprehensive investigation conducted by (Bianchi et al., 2022) revealed the specific role of nitrogen cycling within the oxygen minimum zone (OMZ) of the HNLC Pacific Ocean, highlighting its implications for ecosystem dynamics. A new nitrogen cycling model was developed and specifically designed for integration into an ocean biogeochemical model, aiming to explore the impact of different nitrogen transformation processes on the biogeochemistry of the Eastern Tropical South Pacific OMZ. The model successfully captured typical OMZ conditions, providing valuable insights into the sensitivity of nitrogen transformations to environmental factors. However, it should be noted that while the model explicitly resolved N chemical tracers and their transformations, it did not directly incorporate the microbial communities responsible for these reactions, potentially limiting the understanding of the crucial role of microbes. Additionally, the model did not explicitly account for the conversion of dissolved CO₂ to organic matter through chemolithotrophy due to the relatively small rates compared to organic matter remineralization in the upper ocean, posing challenges in achieving a more comprehensive integration between chemolithotrophy and the carbon cycle.