Modelling the influence of submesoscale processes on phytoplankton dynamics in the northern South China Sea

Introduction

Mesoscale eddies are a ubiquitous feature and important in regulating physical and biogeochemical environments in the ocean (e.g., Falkowski et al., 1991; McGillicuddy et al., 1998; McGillicuddy et al., 2007; Chelton et al., 2011; Siegel et al., 2011; Mahadevan et al., 2012; Gaube et al., 2013; Omand et al., 2015). Regions around mesoscale eddies are known to be associated with strong current shear and strain, which can induce intense upward and downward motions at submesoscale (Lapeyre et al., 2006; Mahadevan and Tandon, 2006). These vertical motions can efficiently transport nutrients into the euphotic zone and lead to elevated biological patches (Lévy et al., 2001; Lehahn et al., 2007; Mahadevan, 2016).

In the ocean, vertical velocities are generally orders of magnitude smaller than horizontal ones (Klein and Lapeyre, 2009), and direct measurements of vertical velocities are not currently available. With the Omega equation, one can diagnose the vertical velocity by using measured density and velocity fields (e.g., Martin and Richards, 2001). Legal et al. (2007) reported an anticorrelation between vertical velocity and density gradient in a frontal region, and suggested that the large-scale strain can dynamically change small-scale filaments. The edge region of a mesoscale eddy or between eddies can be treated as a front because of the horizontal density gradient. Klein and Lapeyre (2009) thus suggested that the vertical exchange should be more efficient at the eddy edge than eddy center.

In addition to fronts, submesoscale eddies that are defined as energetic eddies with a radius smaller than the Rossby deformation radius and a localized structure, have also been observed in different parts of the ocean (McWilliams, 1985; Li et al., 2017; Zhang et al., 2022). They can be very long lived and travel far from their origins. As they can retain their core water mass during their life, they can transport waters with anomalous properties over long distances (Lukas and Santiago-Mandujano, 2001). The cumulative effect of these submesoscale eddies can potentially affect the large-scale transport and distribution of heat, nutrients, and other materials (e.g., Frenger et al., 2018; Gula et al., 2019).

The South China Sea (SCS) is an oligotrophic marginal sea adjacent to the western Pacific Ocean. Seasonally changing monsoon winds play an important role in modulating the upper ocean biology (Liu et al., 2002; Gan et al., 2006). Superimposed on the basin-scale variability, significant spatial and temporal variations of phytoplankton dynamics have been found to be related to mesoscale eddies (Ning et al., 2004; Chen et al., 2007; Huang et al., 2010; Xiu and Chai, 2011; Guo et al., 2015; He et al., 2016; Liu et al., 2017; Xiu et al., 2019). Most of the eddies show high chlorophyll concentrations in the cyclonic eddy core and low concentrations in the anticyclonic eddy core, which is likely related to the mesoscale eddy mechanism (McGillicuddy et al., 1998). A lot of studies have been conducted focusing on physical characteristics of submesoscale features in the northern SCS (Liu et al., 2010; Zhang et al., 2016; Dong and Zhong, 2018; Qiu et al., 2019; Ni et al., 2021; Zhang et al., 2022). However, less attention has been given to the biological impact from submesoscale processes at the eddy edge where enhanced shear and strain exist. Observations from Zhou et al. (2013) and Wang et al. (2018) indicated that possible submesoscale structures around an anticyclonic eddy may enhance phytoplankton production and increase carbon export; however, their samplings were too coarse in space (30-50 km) to resolve detailed submesoscale dynamics. With satellite data in the western SCS, Liu et al. (2017) showed high chlorophyll anomaly present at the northwestern periphery of anticyclonic eddies and suggested that it could be induced by the ageostrophic secondary circulation.

The submesoscale process has a typical spatial scale of 1-10 km and a temporal scale of ~1 day that are difficult to observe by coarse-resolution ship measurements or regular-frequency Argo floats. How these submesoscale processes affect phytoplankton and nutrient distributions in the northern SCS (NSCS), however, remains largely unknown. Moreover, the high energy of the submesoscale field also has significant implications for predictive modeling of oceanic pollutant pathways and concentrations. In this study, a high-resolution physical-biological model was built to investigate the influence of submesoscale features around mesoscale eddies on biological processes in the NSCS.

Data and model

To study biological responses to submesoscale structures, a coupled physical-biological model was developed for the NSCS region. The physical model was based on the Regional Ocean Modelling System (ROMS), and the biological model was based on a modified version of the Carbon Silicate Nitrate Ecosystem (CoSiNE; Chai et al., 2002) model. This modified version of CoSiNE model includes two phytoplankton groups (small phytoplankton with S1 for nitrogen based biomass and Chl1 for chlorophyll concentration, diatom with S2 for nitrogen based biomass, and Chl2 for chlorophyll concentration), two zooplankton classes (microzooplankton (Z1), mesozooplankton (Z2)), two size classes of detritus (small (SPON), large (LPON)), biogenic silica (bSi), nitrate (NO₃), ammonium (NH₄), silicate (SiOH₄), phosphate (PO₄), dissolved inorganic carbon (DIC), total alkalinity (TALK) and dissolved oxygen (DO). In the model, both small phytoplankton and diatom uptake the NO₃, NH₄, and PO₄ for growth, and diatom needs extra nutrient, SiOH₄. The microzooplankton grazes on small phytoplankton, while the mesozooplankton grazes on diatoms, microplankton, and detritus. The mortality and aggregation of phytoplankton and zooplankton are the source terms of detritus. Predation by mesozooplankton and the remineralization are the sink terms of detritus. The detailed model equations and parameters can be found in Ma et al. (2019).

The three-dimensional coupled model was set up for the NSCS (116-120° E, 18-22°N). It has a 1/108° (~ 1 km) resolution horizontally and has 30 vertical levels in terrain-following sigma-coordinates. The coupled model was initialized and one-way nested to the Hybrid Coordinate Ocean Model (HYCOM) dataset that has a spatial resolution of 1/12° and a temporal resolution of 1 day. The HYCOM model uses the Navy Coupled Ocean Data Assimilation system (NCODA) to assimilate available altimeter data, satellite, in-situ profiles from XBTs, buoys and Argo floats, which gives a better representation of the state of the ocean. In addition, eight tidal constituents (M2, S2, N2, K2, K1, O1, P1, and Q1) were used to calculate the hourly tidal elevation at the boundary. Here, M2 is the principal lunar semidiurnal constituent, S2 is the principal solar semidiurnal constituent, N2 is the larger lunar elliptic semidiurnal constituent, K2 is the luni-solar semidiurnal constituent, K1 is the luni-solar diurnal constituent, O1 is the principal lunar diurnal constituent, P1 is the principal solar diurnal constituent, and Q1 is the larger lunar elliptic constituent. For the surface forcing, the 6-hourly surface winds were obtained from the Cross-Calibrated Multi-Platform (CCMP) wind dataset with a spatial resolution of 0.25°. The surface heat and freshwater fluxes were calculated by the COARE3.0 bulk formula using both the CCMP wind and 6-hourly NCEP/NCAR reanalysis data. The initial and boundary conditions for biological variables were derived from a coarse-resolution coupled model, which covers the Pacific Ocean and runs continuously from 1993 to present (Xiu and Chai, 2011).

The coupled model was integrated from 1 January 2015 for one year and daily averaged model outputs were used for analysis. Model outputs were further used to drive an offline particle tracking code (TRACMASS; https://www.tracmass.org/) to examine vertical motions. The TRACMASS code was developed by Döös (1995), which computes numerically the trajectory through each grid cell by solving a differential equation that depends on the velocities on the grid box wall. The chosen water trajectories from any location can be followed along the path both forward and backward in time.

Daily sea-level anomaly (SLA) field with a grid of 1/4° by 1/4° was obtained from the Archiving, Validation and Interpretation of Satellite Data in Oceanography (AVISO). The finite-size Lyapunov exponent (FSLE) product was also obtained from AVISO, which is commonly used to diagnose regions of large stretching and straining by advection (d’Ovidio et al., 2004). The FSLE is defined as the inverse time of separation of two particles from their initial distance to final distance. The particles were advected by altimetry velocities and their trajectories are computed by backward-time integrating the altimetry velocities. Thus, regions with large FSLE generally correspond to the regions where surface current divergence are strong. Chlorophyll concentration data was obtained from a multi-sensor (SeaWiFS, MODIS and MERIS) merged product with a spatial resolution of 4 km and a temporal resolution of 1 day, developed by the Ocean Color Climate Change Initiative (OC_CCI) (Lavender et al., 2015). For analysis, the SLA and FSLE were interpolated onto the chlorophyll data grid.

Results

Submesoscale features between eddies

On 28 April 2015, two cyclonic eddies were present to the west of Luzon Strait, with one located on the north (CE1) and the other one on the south (CE2). There was another anticyclonic eddy (AC1) present to the west of the two cyclonic eddies (Figure 1). The water property in AC1 was similar to that of Kuroshio water, suggesting that the AC1 was generated by Kuroshio intrusion (Shu et al., 2016). When tracing back to February 2015, we can see that the two CEs were generated locally in the NSCS (Shu et al., 2016). The locations of the three eddies were relatively stable between late April and early May. The strength of the two CEs gradually reduced after early May, while the AC maintained its strength and propagated to the west.

Figure 1 Time evolutions of SLA (units: m) overlaid by geostrophic current anomalies on 28 April (A), 04 May (B), 10 May (C) and 16 May 2015 (D).

Satellite-derived surface chlorophyll concentration showed clearly localized high patches with magnitude higher than 0.15 mg m^-3 in the edge regions between CE1 and CE2, and between AC1 and the CEs (Figures 2A, B). Due to cloud contamination, we can only show four-day-averaged chlorophyll concentration, which is not fine enough to resolve detailed submesoscale dynamics. For the daily chlorophyll data, the mean cloud coverage in the study region in April and May was about 76% and it reduced to about 31% when using the four-day-averaged data. Unlike physical variables, it usually takes days for phytoplankton to show biological changes in response to dynamical forcings. Four-day-averaged chlorophyll is probably able to reflect biological changes to submesoscale dynamics, but locations might not be consistent with physical variables (Mahadevan, 2016; Lévy et al., 2018). The spatial decoupling between upwelling, phytoplankton new production, and export production across a submesoscale front has been reported (Estapa et al., 2015). Large horizontal stretching and straining represented by large FSLE calculated from altimeter data were found in the edge area between these eddies (Figures 2C, D). These regions are very dynamic where unstable and stable manifolds cross each other and are prone to frontogenesis that is accompanied with large vertical velocities (>10 m d^-1; Lehahn et al., 2007; d’Ovidio et al., 2009).

Figure 2 Surface chlorophyll concentration (units: mg m^-3) overlaid by SLA (units: m; black contours for positive and blue ones for negative values) in (A, B), and overlaid by FSLE contours larger than 0.1 d^-1 in (C, D) during 5-12 May. The red line in (B) denotes the position of the transect for further analysis.

A transect through the edge region between CE1 and CE2 also showed that chlorophyll concentration in the CE was generally higher than that in the AC. The highest chlorophyll was present in the edge region with small negative SLAs and large FSLE values (Figure 3). High strain between eddies may create submesoscale upwelling/downwelling motion that facilitates nutrient transport vertically. Because the surface NSCS in May is generally in a nutrient-limited condition, submesoscale nutrient transport is likely to stimulate phytoplankton growth and generate localized high chlorophyll patches.

Figure 3 Spatial variability of SLA (units: m), surface chlorophyll concentration (units: mg m^-3), and FSLE (units: day^-1) along the transect shown in Figure 2B.

To investigate the existence of upwelling, movements of passive particles were tracked in the model. Although the model missed the SLA magnitudes slightly, it reproduced CE1, CE2 and AC1 reasonably well (Figure 4). On 4 May, over 50,000 passive particles were released in the model at 50 m in the area between CE1 and CE2, covering both the core and edge regions (Figure 4). The movements of these particles were tracked over time in the model. These particles were first stretched along the hyperbolic region and then advected away mostly at the eddy edge. Only a relatively small number of particles stayed in the eddy core as it propagated. We tracked those particles that entered the mixed layer after seven days from their release, and found that they were generally the ones released at the eddy edge regions (Figure 5). Consistent with this pattern, the probability of particles staying in the mixed layer during the one-month period is particularly higher at the eddy edge regions (Figure 6).

Figure 4 Comparisons of spatial distributions of SLA (units: m) from altimeter data (A, B) and the model outputs (C, D) on 28 April and 8 May. The red rectangles in (A, B) depict the domain of the model.

Figure 5 Lagrangian particle releasing results on 04 May (A), 08 May (B), 12 May (C) and 16 May 2015 (D). The particles were released at the surface of a depth of 50 m. The black dots show the particle locations on different days. The shadings are SLA (units: m). The magenta circles in (A) are the particles being brought into the mixed layer 7 days after being released.

Figure 6 The spatial distribution of particle probability of staying in the mixed layer in May 2015 (A) and the spatial distribution of modeled mixed layer depth (m) on 4 May 2015 (B). The particles were released at 50-m depth on 4 May. The contours are the SLA isolines with solid and dashed lines denoting positive and negative values, respectively.

Driving mechanisms

The submesoscale nutrient transport was further examined in the model. On May 8, the nitrate concentration at 50 m showed localized high patches (Figure 7). These patches were at a submesoscale length scale and present around mesoscale eddies. To examine the physical processes driving these nitrate patches, two cases were chosen. We used two moving boxes (50×50 km) to encompass two nitrate patches and follow their movements (white box for case1, red box for case2).

Figure 7 Spatial distributions of nitrate concentration (mmol m^-3) at 50 m from the model outputs. The white and red rectangles mark the locations of two high nitrate patches for analysis.

For case1, the high nitrate patch was associated with a coherent submesoscale eddy that can be captured from surface SLA (Figure 7). This eddy has a radius of ~20 km. It was formed near the Dongsha Island, probably related to the current-topography interactions (Zhang et al., 2022). From 5 May, it was advected to the southwest by the current at the edge of mesoscale eddies. This submesoscale eddy associated with high nitrate concentration eventually merged into the mesoscale cyclonic eddy, CE2, providing a significant contribution of nitrate input flux to CE2. While the submesoscale eddy propagated, the distribution of vertical velocity (w) at 80 m displayed a diapole pattern with positive w at the leading edge and negative w at its trailing edge (Figure 8). The similar distribution pattern of current divergence (δ=u_x+v_y ; where u, v are the velocity components in x, y, i.e., east and north, directions) with positive and negative values corresponding to upwelling and downwelling, respectively, suggests that current divergence is the possible mechanism leading to the upwelling and downwelling processes (Figure 8).

Figure 8 Distributions of vertical velocity (A units: m d^-1) at 80 m and normalized current divergence (A) at 50 m following the propagation the submesoscale eddy of case1. The arrows in (B) illustrates the general moving directions of the box.

Following the propagation of the submesoscale eddy, both positive and negative w decreased with time in a similar manner, probably suggesting a connected circulation cell (Figure 9A). Nevertheless, the box-averaged w stayed positive during the eddy propagation, which also drives the positive nitrate flux vertically Figure. The vorticity (ζ=v_x−u_y ) normalized by the Coriolis frequency (f) indicated that the submesoscale eddy experienced both the developing stage before 29 April and the decay stage afterwards (Figure 9B). During its developing stage, nitrate level at 50 m was relatively stable and it started to increase during its decay stage while propagating (Figure 9C). During its propagation, the net nitrate flux was generally positive due to the relatively smaller negative horizontal flux, which resulted in the accumulation of nitrate concentration in the upper layer (Figure 9D). Following the movement of the submesoscale eddy, modeled mean chlorophyll concentration in the upper 50 m increased from 0.1 mg m^-3 to 0.24 mg m^-3, which was in a similar magnitude as the satellite data.

Figure 9 Temporal variations of box-averaged vertical velocities (A; units: m d^-1), normalized vorticity by f (B), nitrate concentration (C; mmol m^-3), and nitrate fluxes (D; units: mmol s^-1) following the movement of case1 patch. Here, the horizontal nitrate flux was calculated from surface to 80 m and the vertical flux was calculated at 80 m.

For case2, the high nitrate patch that can reach higher than 2.0 mmol m^-3 was not associated with any coherent circulation structures (red box in Figure 7). It showed up on 5 May and was eventually stretched away by the horizontal current between mesoscale eddies. The box averaged nitrate concentration showed a similar temporal pattern (Figure 10A). The increase of nitrate from 5 May was induced by the positive net nitrate flux (horizontal plus vertical flux) into the box that was largely driven by the increase of positive w (Figures 10B, C). On 4 May, there was an increase of north-south density gradient (Figure 10E), which can consequently lead to enhanced horizontal strain ($St=\sqrt{{\left(u_x-v_y\right)}^2+{\left(v_x+u_y\right)}^2}$) and further induce secondary circulation with strong upward water motions (Figures 10B, D). This mechanism is consistent with frontal dynamics. We further composited the w at 80 m from 6 May to 12 May when the case2 patch moved generally in a zonal direction along the front between a mesoscale anticyclonic eddy in the north and a mesoscale cyclonic eddy in the south. The spatial distribution pattern of composited w demonstrated positive values in the north and negative values in the south in correspondence to the low density in the north and high density in the south, respectively (Figure 11). Therefore, the elevated nitrate concentration in case2 was mainly induced by the frontal dynamics. The modeled chlorophyll change in the case2 patch generally followed the change of nitrate concentration over time.

Figure 10 Temporal variations of box-averaged nitrate concentration (A; unit: mmol m^-3), vertical velocities (B; units: m d^-1), net nitrate advection flux (C; unit: mmol s^-1), normalized horizontal strain by f (D), and north-south density gradient (E; unit: kg m^-4) of case2 patch.

Figure 11 Composted mean of modeled vertical velocity (unit: m d^-1) at 80 m for case2 during 6 May~12 May, 2015.

Discussion and conclusions

Mesoscale eddies are known to induce perturbations in biogeochemistry in the eddy core through mesoscale dynamics. In addition to mesoscale responses, high-resolution chlorophyll images often show patchy distribution patterns that vary over a distance of a few kilometers at a temporal scale of days, which are suggested to link to submesoscale dynamics. At submesoscale, the spatial and temporal variations of tracers such as phytoplankton chlorophyll is often related (Mahadevan, 2016). Thus, to resolve submesoscale structures, both spatial and temporal resolutions need to be fine enough.

Two mechanisms have been previously suggested to generate submesoscale structures. One is the mesoscale-driven frontogenesis (Lapeyre et al., 2006; Mahadevan and Tandon, 2006), and the other is the mixed-layer instability (MLI; Boccaletti et al., 2007). The vertical scale of MLI is the mixed layer depth. The mixed layer depth during the study period was generally shallower than 40 m, except in the regions where anticyclonic eddies were present (Figure 6B). From the model results, we found strong fluctuations of isopycnals and upwelling/downwelling motions below the mixed layer. These fluctuations of isopycnals were conspicuously observed in the edge region between mesoscale eddies, where enhanced shear and strain can sharpen existing horizontal density gradient and give rise to submesoscale upwelling/downwelling motions. It is thus likely that the submesoscale processes were induced by frontogenesis and baroclinic instabilities (Ramachandran et al., 2014; Zhang et al., 2021). Other studies also indicated that the frontogenesis plays a key role in generating submesoscale processes in the SCS (Zhong et al., 2017).

By simplifying the omega equation in a frontal region, we can obtain its two-dimensional version, that is,

This omega equation is derived by assuming that the x variation (along the front) in the data is much smaller than its y variation (across the front). It states that mesoscale flow can drive the growth of submesoscale density gradient associated with vertical velocity. With further simplification of this equation, Legal et al. (2007) derived an expression indicating the robust correlation between vertical velocity and density anomaly. It is consistent with the findings of this study, in which intense vertical velocity was illustrated by particle trajectories between eddies. The modeling results for submesoscale case2 have shown that such vertical motion is often linked with the vertical transport of tracers, such as nutrients, consistent with previous studies (Lévy et al., 2001; Lapeyre and Klein, 2006).

In addition, the modeling results showed that submesoscale eddy could be another mechanism inducing strong vertical nutrient transport at submesoscale. Although case1 was not generated at the eddy edge, it was carried away by the flow around the mesoscale eddy after formation. During its movements, divergence induced upwelling at the leading edge kept injecting nutrient from subsurface to the upper layer. Previous satellite study in the western SCS has observed such small-scale eddies associated with high surface chlorophyll concentrations (Liu et al., 2015). The submesoscale eddy could also be a potential nutrient source to mesoscale eddy once they merge together.

Biological responses to submesoscale processes are more complicated than physical tracers. Nutrients upwelled from subsurface can be taken up by different phytoplankton species that are further modulated by ecosystem dynamics and other biogeochemical processes, which are all subject to different temporal scales (Lévy et al., 2018). The short time scales associated with strong submesoscale vertical nutrient injections and the biological response time favour dominance by the large phytoplankton size class (diatoms) of the model (Figure 12; Lévy et al., 2012; Guo et al., 2022). Different physical dynamics between mesoscale and submesoscale processes can thus lead to spatial and temporal heterogeneities in community structure and ecosystem dynamics (d’Ovidio et al., 2010; Lévy et al., 2012; Clayton et al., 2014; Mousing et al., 2016; Clayton et al., 2017). Moreover, the intensity of nutrient injection induced by submesoscale processes is also determined by the background nutricline depths that are set by mesoscale or large-scale processes. Consequently, regions with high physical straining and stretching are not always associated with enhanced total chlorophyll concentration (e.g., Guo et al., 2022).

Figure 12 Modeled diatom ratio (%) in total phytoplankton averaged in the upper 50 m on 9 May 2015. The red and white rectangles depict two submesoscale nitrate patches. The contours are the SLA isolines with solid and dashed lines denoting positive and negative values, respectively.

In conclusion, with a high-resolution coupled model, we examined submesoscale processes between mesoscale eddies in the NSCS. By calculating the FSLE from altimeter data and tracking Lagrangian particles in the three-dimensional model, the submesoscale upwelling motion was revealed between eddies. We found that both frontal dynamics and submesoscale coherent eddies can inject subsurface nutrients into the upper layer, stimulate phytoplankton growth, and alter community structure. As the coupled model didn’t include data assimilation, model discrepancies were found in simulating the strength and locations of mesoscale and submesoscale features. This study can serve as a process-orientated study to investigate potential dynamics and impacts of submesoscale features on biology. This study shows the strong physical-biological links at submesoscale in the NSCS, highlighting the importance of submesoscale processes on marine ecosystem and the need for high-resolution observations.

Data availability statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author contributions

Conceptualization: PX, WM; Methodology: LG, PX. Writing: PX. Review and editing: LG, WM. Funding acquisition: PX. All authors contributed to the article and approved the submitted version.

Funding

This research was supported by the Key Special Project for Introduced Talents Team of Southern Marine Science and Engineering Guangdong Laboratory (Guangzhou) (GML2019ZD0305) and the National Natural Science Foundation of China (41890805, 41730536).

Acknowledgments

The Ocean Color Climate Change Initiative (OC-CCI, Version 3.1) data were obtained from http://www.esa-oceancolour-cci.org/. The Cross-Calibrated Multi-Platform (CCMP) Version-2.0 data were obtained from www.remss.com. Six-hourly NCEP/NCAR Reanalysis surface forcing data were obtained from https://www.esrl.noaa.gov/psd/. The SLA and FSLE data were obtained from https://www.aviso.altimetry.fr/.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References