Interactive comment on “Challenges in modelling spatiotemporally varying phytoplankton blooms in the Northwestern Arabian Sea and Gulf of Oman”

The manuscript "Challenges in modelling spatio-temporally varying phytoplankton blooms in the Northwestern Arabian Sea and Gulf of Oman" by Sedigh Marvasti et al. is a fundamentally solid model-data intercomparison. Having said that, I do find the manuscript somewhat unbalanced. Most importantly, it lacks analysis about the origin of the asymmetry between the large September bloom and the small February bloom. I think a real fundamental insight about the controls on seasonal phytoplankton blooms could be gained from such an analysis. Furthermore, the current analysis focuses too much on the possible roles of cyclonic and anti-cyclonic eddies, culminating in a long-winded speculation without a satisfactory resolution. Overall, I think that this C4281


Introduction
The region of northwestern Arabian Sea and the Gulf of Oman (15-26 • N, 56-66 • E) is a highly productive region (Madhupratap et al., 1996;Tang et al., 2002), with satellite estimates of carbon export of 137 gC m −2 yr −1 , much higher than the ∼ 80 gC m −2 yr −1 found in the Subpolar North Atlantic and Pacific (Dunne et al., 2007).Peak chlorophyll a concentrations exceed 0.7 mg m −3 in this region (Fig. 1a).
This region may be changing in important ways.In both the Persian Gulf and the Gulf of Oman, there is evidence that harmful algal bloom (HABs) and their impacts are increasing (Richlen et al., 2010).HAB occurrences have been more frequently reported in the Gulf of Oman than in the Persian Gulf.A total of 66 red tide events (mostly dominated by Noctiluca scintillans) have been recorded between 1976 and 2004 including 25 blooms resulting in mass mortality of fish and marine organisms.Reasons for the increase in blooms include aquaculture, industrial and sewage inputs, natural dispersal and human-aided transport, long-term increases in Published by Copernicus Publications on behalf of the European Geosciences Union.nutrient loading and global expansion of species (Richlen et al., 2010) as well as global climate change (Goes et al., 2005).The latter paper suggested that increasing blooms were driven by an increase in the strength of the Asian monsoon.
Evaluating such a possibility and extending it into the future requires the use of Earth System Models.However, such projections will only be as good as the models on which they are based.In this paper we examine several models run at the Geophysical Fluid Dynamics Laboratory in the Arabian Sea.We consider numerical results from five different 3-D global Earth system models, which we denote CORE-TOPAZ, Coupled-TOPAZ, Coupled-BLING, Coupled-miniBLING, and the Geophysical Fluid Dynamics Laboratory Climate Model version 2.6 (CM2.6 miniB-LING).The first two of these models use the relatively complex TOPAZ biogeochemistry, but have low resolution and do not resolve eddies, the third has a simplified biogeochemistry (BLING, Galbraith et al., 2010) which does not carry phytoplankton biomass as a separate variable while the last two models have an even simpler biogeochemistry that does not directly simulate dissolved organic matter or iron cycling (miniBLING, Galbraith et al., 2015).Only the final model resolves eddies.
The seasonal cycle is an important metric for models to be able to simulate.The Arabian Sea is influenced by a reversing monsoonal cycle (Wang and Zhao, 2008), an evaporative fresh-water flux over most of the basin, and an annual mean heat gain (Banse and McClain, 1986;Fischer et al., 2002).In summer (June-September), the southwest Monsoon (SWM) blows strongly across the northwestern Arabian Sea (Al-Azri et al., 2010).Driven by a land-sea pressure gradient, the SWM is a large-scale feature of the atmospheric circulation of the tropics, extending from a surface pressure high near 30 • S in the Southern Hemisphere northward to the surface low over Asia (Anderson and Prell, 1993).During the SWM, winds are steered by the East African highlands to form a strong low level atmospheric jet, referred to as the Findlater Jet (Bartolacci and Luther, 1999;Honjo et al., 2000), which crosses the Equator over the Indian Ocean and blows over the Arabian Sea parallel to the Omani coastline in a northeast direction (Honjo et al., 2000).The orientation of the Findlater Jet parallel to Omani coast leads to coastal upwelling along the coast and downwelling on the eastern side of the Jet in the middle of the Arabian Sea.This upwelling provides nutrients to the surface layer (Fig. 1b) (Al-Azri et al., 2013;Kawamiya and Oschlies, 2003;Madhupratap et al., 1996;Murtugudde et al., 2007;Veldhuis et al., 1997;Wang and Zhao, 2008).The SWM does not destabilize the surface layers, which are fairly stable in northern summer (Fig. 1c).
The Northeast Monsoon (NEM), which happens from December through February, is not as strong as the SWM (Dickey et al., 1998;Shalapyonok et al., 2001;Veldhuis et al., 1997).Ocean surface wind stress is lower (0.032 N m −2 in NEM compared to 0.127 N m −2 in SWM), and does not lead to upwelling like the SWM along the Omani coast.However, negative heat flux results in a destabilizing buoyancy flux, subsequent convective overturning (Barimalala et al., 2013;Kawamiya and Oschlies, 2003), and deepening and cooling to a depth of ∼ 60 m (Fig. 1c, d).This brings up nutrients and fuels a wintertime bloom.In addition, as shown in Fig. 1d in wintertime bloom the mixed layer depth (MLD) is deeper than summer.
A second metric of the bloom dynamics is the relationship between the blooms and mesoscale eddies (Al-Azri et al., 2013;Dickey et al., 1998;Hamzehei and Bidokhti, 2013;Shalapyonok et al., 2001;Goes et al., 2005).The confluence of the Persian Gulf outflow current and the East Arabian Sea Current parallel to Omani and Yemeni coastlines in Arabian Sea leads to formation of a frontal zone and formation of persistent eddies in the region.Because the size of eddies is comparable to the width of the Gulf of Oman, they can affect mixing and transport of biota on a basin scale (Fischer et al., 2002;Piontkovski et al., 2012).Piontkovski et al. (2012) suggested that the increased amplitude of the seasonal cycle of chlorophyll a might be associated with the increased variability of mesoscale eddy kinetic energy (EKE) per unit mass in the Gulf of Oman or in the western Arabian Sea.Gomes et al. (2008) noted potential anticorrelation between sea surface height and chlorophyll, but did not find a consistent relationship over time.Gaube et al. (2014) provide a global overview of how eddies influence chlorophyll blooms.They find that the effect of mesoscale eddies on the chlorophyll bloom varies both temporally and spatially.They identify four particular mechanisms that can be distinguished by linking sea surface anomalies to chlorophyll, namely eddy stirring, trapping, eddy intensification, and Ekman pumping.Although Gaube et al. (2014) find a negative correlation between chlorophyll and SSH in the Arabian Sea, they do not analyze which of these mechanisms is involved in this region, nor do they quantify the extent to which this correlation varies over the course of the season.Resplandy et al. (2011) indicated that the spatial variability associated with mesoscale eddies in the Arabian Sea produces spatial variability in the bloom and that another source of variability is found to be restratification at these structures.Advection from coastal region is identified as the mechanism providing nutrients in summer, while vertical velocities associated with mesoscale structure are found to increase the overall nutrient supply.However, this work does not make clear how the spatial distribution of the eddy nutrient supply is related to the eddies, not whether this relationship is the same in all seasons.
The structure of this paper is as follows: all data sets including ocean color data and altimeter data are explained in Sect. 2 of the paper along with the specification of five different 3-D global Earth system models.In Sect.3, the remote-sensing results are used to study the spatiotemporal variability of chlorophyll a in mesoscale structures in the study region.We find a seasonal relationship between SSHA and chlorophyll such that cyclonic eddies are associated with blooms, but only during the winter.This means that interannual variability in blooms will be shaped by mesoscale eddy activity and may not be predictable.Results of the 3-D global Earth system models are discussed in Sect. 4. Annual cycles of variation of chlorophyll a and nutrients for all GFDL models within the whole region are compared against the corresponding satellite results and field measurements.The models tend to overestimate wintertime productivity, in large part due to excessive mixing.They also fail to explain the bloom-SSHA relationship except in a few special cases.We argue that the eddies act to modulate turbulent mixing of nutrients to the surface -a mechanism not emphasized in previous literature.However, this can only occur if there is a strong and relatively shallow nutricline.Since the model only simulates such a feature in the Southern Arabian Sea, it does not capture the observed relationship between SSH and biology.Both the overestimation of the wintertime bloom and the failure to predict its modulation by eddies can thus be traced to difficulties in modeling the stratification of the Northwest Arabian Sea, most likely as a result of a failure to properly simulate overflows.

Satellite products
We examine the relationship of blooms and eddies using the GSM5 Maritorena et al. (2002) product based on the Sea-WIFS (Sea-viewing Wide Field-of-view Sensor) ocean color data and Sea Surface Height Anomaly (SSHA), based on altimeter data acquired from the Archiving, Validation and Interpretation of Satellite Oceanographic (AVISO) Data Center (http://www.aviso.oceanobs.com).The SSH anomaly is calculated relative to the annual cycle.
The GSM algorithm represents the normalized water leaving radiance L wN (λ) at multiple wavelengths as a nonlinear function, as following Maritorena et al. (2002), where t is the sea-air transmission factor, F 0 (λ) is the extraterrestrial solar irradiance, n w is the index of refraction of the water, seawater backscatter b bw (λ), absorption a w (λ), a * ph is the chlorophyll a (chl) specific absorption coefficient, S is the spectral decay constant for absorption by chromophoric dissolved organic materials (CDOM), η is the power-law exponent for the particulate backscattering coefficient, and λ 0 is a scaling wavelength (443 nm).The cdm absorption coefficient [a cdm (λ 0 )], and slope factor S then determine the absorption across a range of wavelengths while the particulate backscatter coefficient [b bp (λ 0 )] and coefficient η constrain the scattering.Letting λ 0 be 443 nm assuming that all terms other than chl, [a cdm (λ 0 )] and b bp (443 nm) are constant, one can then use the normalized water leaving radiance to invert for chl, a cdm , and backscatter b bp .One limitation of this approach is that if the inherent optical properties vary with time or space, this variation will introduce errors into the estimate.Following Behrenfeld et al. (2005), we convert the backscatter coefficient into units of particulate carbon biomass using the relationship p carb = 13 000 (b bp − 0.00035).

S. Sedigh Marvasti et al.: Challenges in modeling spatiotemporally varying phytoplankton blooms
Satellite-based remote sensing is the only observational method suitable for measuring physical and biological properties over large regions of the ocean.However, satellite ocean color and SST are limited to surface distributions and provide no information about the vertical structure within the ocean (McGillicuddy et al., 2001).Additionally acquiring data requires cloud-free viewing of the ocean surface, which as we will see is a problem in this region at certain times of the year.This lack of information motivates our examination of numerical models, which ideally could be used to provide estimates of the ocean state when observations are sparse as well as to extrapolate both vertically and into the future.

Numerical models
Numerical results are presented in this paper based on the output of five different 3-D global Earth system models, which we denote CORE-TOPAZ, Coupled-TOPAZ, Coupled-BLING/miniBLING and GFDL CM2.6 (miniB-LING).The first two of these models use the relatively complex TOPAZ biogeochemistry, but have low resolution and do not resolve eddies.The third and fourth use two simplified biogeochemistry codes (BLING and miniBLING) which do not carry phytoplankton biomass as a separate variable while the last model has very high resolution and uses the miniBLING simplified biogeochemistry.Below, we describe the different physical models, followed by a summary of the biogeochemical codes run within these models.

Physical model description
The ocean-ice model used in the CORE-TOPAZ model follows the corresponding components of the GFDL CM2.1 global coupled climate model (Delworth et al., 2006).The vertical resolution ranges from 10 m over the top 200 m to a maximum thickness of 250 m at 5500 m depth with 50 layers in all.The meridional resolution is 1 • , whereas the zonal resolution varies between 1 • in mid-latitudes and 1/3 • at the equator.North of 65 • , a tripolar grid is employed to avoid singularity arising from convergence of meridians at the North Pole.Up-to-date parameterizations of mixed-layer dynamics, isopycnal mixing, advection by subgridscale eddies, bottom topography, bottom flows, and lateral viscosity are included -for more detail see Griffies et al. (2005) and Gnanadesikan et al. (2006).Both the dynamics and thermodynamics sea ice are simulated with five thickness classes of sea ice being resolved.
In the CORE-TOPAZ model, surface forcing is set using the Coordinated Ocean-ice Reference Experiment (CORE) protocol (Griffies et al., 2009), where the inputs for calculating surface fluxes are taken from an atmospheric analysis data set adjusted to agree better with in situ measurements.Sensible and latent heat fluxes are then calculated using bulk formulae.Freshwater forcing is given by a combination of applied precipitation, evaporation computed us-ing bulk fluxes, and a correction diagnosed to restore surface salinities in the top 10 m to climatological monthly values over 60 d.Hence, the fluxes forcing the CORE runs could be thought of as "best guess" observationally based estimates.Such a prescription omits important feedbacks whereby the atmosphere ensures that rainfall and evaporation are consistent with each other, although the restoring correction is a crude representation of these feedbacks.We use the version of the model described in Gnanadesikan et al. (2011), which analyzed different modes of interannual variability in biological cycling across the North Pacific Ocean.
The Coupled-TOPAZ model corresponds to the control simulation of the GFDL ESM2M submitted as part of the IPCC AR5 process (Dunne et al., 2012).In this model the ocean is coupled to the atmosphere, land, and sea ice components.Gnanadesikan et al. (2014) discuss the behavior of this model in the North Atlantic, but its behavior in the Arabian Sea has not been previously analyzed.Two additional versions of this model, referred to here as Coupled-BLING/miniBLING, were run using the BLING and mini-BLING biogeochemical models described below, but with the light field given by the TOPAZ code.The differences between the 1 • models highlight differences due to biological formulation.
The ocean component of ESM2M employs the MOM4p1 code of Griffies et al. (2009) which largely mimics the CM2.1 ocean (identical horizontal and vertical resolution and parameterization of mixing).However, ESM2M ocean uses a rescaled geopotential vertical coordinate (z * ; Adcroft et al., 2004;Stacey et al., 1995) for a more robust treatment of free surface undulations.The ESM2M implementation includes updates to the K-profile parameterization (Large et al., 1994) based on Danabasoglu et al. (2006), as well as model-predicted chlorophyll modulation of short-wave radiation penetration through the water column.ESM2M also includes completely novel parameterizations relative to CM2.1, such as parameterization of submesoscale eddyinduced mixed layer restratification (Fox-Kemper et al., 2008).Instead of prescribed vertical diffusivity for interior mixing (Bryan and Lewis, 1979), ESM2M employs the Simmons et al. ( 2004) scheme along with a background diffusivity of 1.0 × 10 −5 m 2 s −1 in the tropics and 1.5 × 10 −5 m 2 s −1 poleward of 30 • latitude following a tanh curve.
The Geophysical Fluid Dynamics Laboratory Climate Model version 2.6 (CM2.6) is a high-resolution eddyresolving model.This model has the same atmosphere model and ocean Physics as CM2.5 (Delworth et al., 2012).CM2.6's ocean component has higher horizontal resolution than CM2.5, with grid spacing, which is changeable from 11 km at the equator to less than 4 km at very high latitudes.This means that the model is capable of resolving eddy features in the tropics, as we will see below.

Biogeochemical cycling codes
The TOPAZ code (Tracers of Ocean Productivity with Allometric Zooplankton code of Dunne et al., 2010), keeps track of five inorganic nutrients used by phytoplankton: nitrate and ammonia, inorganic phosphate, silicate, and dissolved iron.Additionally, the model carries three other dissolved inorganic tracers: dissolved inorganic carbon, alkalinity and dissolved oxygen.Based on the work of Dunne et al. (2007), the model also keeps track of fine lithogenic material, which plays a role in ballasting organic material and delivering it to the sediment (Armstrong et al., 2002;Klaas and Archer, 2002).The five inorganic nutrients are taken up in different ways by three classes of phytoplankton: small, large and diazotrophic.A comprehensive description of TOPAZ v2 can be found in the Supplement of Dunne et al. (2013).
TOPAZ is unusual among comprehensive Earth System Models in that it uses a highly parameterized version of grazing.Instead of grazers being explicitly simulated, grazing rates are simply taken as a function of phytoplankton biomass, with different power-law dependence for small and large phytoplankton.The grazing formulation was fit to about 40 field sites to produce a size structure that transitions realistically between being dominated by small phytoplankton and low particle export ratio at low levels of growth and large phytoplankton and high particle export ratio in nutrient and light-replete conditions.At equilibrium, the resulting parameterization produces biomass that is a function of growth rate (linear for small plankton, cubic for large).A similar scaling in particle size spectrum was seen across ecosystems by Kostadinov et al. (2009).In contrast to models that explicitly simulate zooplankton, TOPAZ does not depend on poorly known zooplankton behavioral parameters (such as handling efficiency or grazing half-saturation) or on the details of how different trophic levels interact.
Even though it does not simulate zooplankton explicitly, TOPAZ still carries over two dozen tracers, making it extremely expensive to run in high-resolution simulations.For this reason Galbraith et al. (2010) developed the Biogeochemistry with Light Iron Nutrients and Gasses (BLING) model, which parameterizes the entire ecosystem.The original version of BLING has only five explicit tracers: dissolved inorganic phosphorus (PO 4 ), dissolved organic phosphorus (DOP), dissolved Iron (Fe), dissolved inorganic carbon (DIC), and oxygen (O 2 ).It includes the impacts of macronutrient and micronutrient limitation and light limitation on phytoplankton by using these to calculate a growth rate.Using the same machinery as TOPAZ, it then uses this growth rate and implicit treatment of community structure to estimate phytoplankton biomass, and uses this biomass to calculate the rate at which nutrient is taken up by plankton and cycled through the ecosystem.
The miniBLING code (Galbraith et al., 2015) represents a further simplification.In this model the iron field is taken from a lower-resolution version of the model (an approxi-mation which has limited impact in the Arabian Sea, where phytoplankton is generally not iron-limited) and so Fe is not treated prognostically.Additionally the DOP pool is eliminated.Simulations using the ESM2M physical model show that control simulations of oxygen and surface nutrients produced by the miniBLING and BLING models are very similar to those produced in the same model with TOPAZ (Galbraith et al., 2015).Galbraith et al. (2015) also show that BLING and miniBLING simulate very similar patterns of oxygen change and anthropogenic uptake in a simulation where CO 2 is increased by 1 % per year until it is twice the preindustrial concentration.
It should be noted that simplified BLING and miniBLING codes neglect some processes that may be important.Only nonliving components are advected and mixed by the ocean circulation, which could result in inaccurate distribution of biology in frontal regions at high resolution.Additionally, as will be discussed below, the lack of a biomass variable may lead to overestimating how rapidly plankton inventories can grow.Also, the rich behavior of the nitrogen cycle with its interaction with iron, phosphorus and oxygen cannot be simulated with one macronutrient tracer (Behrenfeld, 2010).Specifying iron limitation, as done in miniBLING, may also have some impacts in our region.As extensively discussed by Naqvi et al. (2010) there is a possibility of iron limitation over the southern parts of the Omani shelf and in the offshore region during the latter part of the Southwest Mansoon, which can result in high nitrate-low chlorophyll conditions.The western equatorial and southern tropical region of the Indian Ocean are iron-limited and the Arabian Sea (southern parts) may become iron-limited under strong upwelling conditions (Wiggert et al., 2006).

Annual cycle and interannual variability
We begin by using the GSM5 satellite data to examine the annual cycle and interannual variability in two different regions, the whole NW Arabian Sea (56-66 • E, 15-26 • N) and a smaller region including the Gulf of Oman, (60-62 • E, 22-26 • N).As shown in Fig. 2a to c for whole region, clear annual cycles of chlorophyll a, backscattering and CDOM are observed.Even larger annual cycles of variation of chlorophyll a, backscattering and CDOM are seen in the smaller region, as shown in Fig. 2d to f.More pronounced interannual variability is observed in the smaller region as opposed to the larger region.
The annual variations of all parameters are broadly consistent with each other.The maximum values associated with the summer bloom are generally seen in September, with values of 1.0 mg m −3 , 50 mgC m −3 , and 0.1 m −1 for chlorophyll, particulate carbon and CDOM, respectively, within the whole region.Within the smaller region, the values are 1.25 mg m −3 , 65 mgC m −3 , and 0.125 m −1 for chlorophyll, particulate carbon, and CDOM, respectively.For 2 years of 2001 and 2002, the particulate carbon values (∼ 90 mgC m −3 ) are much higher than the average of the other months over both regions, but the chlorophyll does not show pronounced peaks.A winter bloom is also pronounced in February as a second maximum in a yearly cycle, where the magnitudes are about 0.07 mg m −3 , 40 mgC m −3 , and 0.07 m −1 for chlorophyll, particulate carbon and CDOM, respectively, within the whole region, and about 0.09-1.5 mg m −3 , 55-80 mgC m −3 , and 0.11-0.14m −1 for chlorophyll, particulate carbon and CDOM, respectively, within the smaller region.That the summer bloom in the both regions is stronger than the winter bloom has been discussed by Al-Azri et al. (2010) and Levy et al. (2007).

Variability of chlorophyll a in mesoscale structures
Mesoscale structures can be seen in the Northwest Arabian Sea in both the SeaWiFS chlorophyll a distribution and AVISO sea surface height anomaly.Over the course of 2001 (Fig. 3), both a summer bloom (which most likely starts in August and ends in ∼ October) and a winter bloom (which starts in January and goes away in April) can be seen in chlorophyll a.In March, the last month of the winter bloom, chlorophyll a concentrations are high over the entire region in both the anticyclones (warm eddies with positive SSHA) and the cyclones (cold eddies with negative SSHA).The observed bloom in March terminates abruptly in April, although the observations show that eddies are still active in the region.In June, July and August, the satellite ocean color data are not available due to excessive cloudiness.In September, the last month of the summer bloom, most of the region including cyclones and anticyclones and coastal regions had high chlorophyll a concentration.However in the following months the bloom persists only within the cold eddies and disappears over the warm eddies (a phenomenon also seen in Sargasso Sea by McGillicuddy et al., 2001).The relationship between sea surface chlorophyll a and eddies for the other years between 1998 and 2005 during the month of November, is shown in Fig. 4. The relationship between blooms and SSHA is clear and striking.Note particularly the difference between 1998 and 2001, when the location of high and low chlorophyll regions switches relative to the Ras al Hadd.This difference in bloom location is perfectly reflected in the different locations of the eddies.

Chlorophyll-Sea Surface Height Anomaly (SSHA) cross-correlation
The seasonal relationship between chlorophyll and SSHA can be seen in the monthly variation of the spatial crosscorrelation between the two variables over the entire northwest Arabian Sea.Chlorophyll-SSHA cross-correlations between 1998 and 2005 in the satellite data are shown in Fig. 5a.To check that the chlorophyll results are not an artifact of the remote-sensing inversion, two other related parameters, the backscattering coefficient (BBP) and chromophoric dissolved organic matter (CDOM) are also crosscorrelated with SSHA, as depicted in Fig. 5b and c After the winter bloom (typically April and May), the cross-correlation is positive or very small, which suggests no relation between the mesoscale eddies and the blooms.As discussed in Kumar et al. (2001), low primary production is observed after termination of winter cooling during Spring Inter-Mansoon (SIM) (see also Gomes et al., 2008).This result would be also consistent with SIM producing weak atmospheric forcing in the region.
In contrast, after the summer bloom (typically October-December) as the average values of chlorophyll a decrease, chlorophyll and SSHA become relatively highly anticorrelated.The reason for the anti-correlation is the persistence of chlorophyll at the regions with negative SSHA that are typically considered to be cyclonic (cold) eddies and disappearance of chlorophyll a in positive SSHA that are assumed to be anti-cyclonic (warm) eddies.Particle backscatter also provides almost same cross-correlation suggesting that the chlorophyll a signal does not result purely from photoadaptation.Moreover, the CDOM-SSHA cross-correlation shows the same annual cycle although with smaller peak values.
The spatial relationship between blooms and eddies seen in the Northern Arabian Sea can be compared with the patterns noted by Gaube et al. (2014).Their eddy stirring mechanism involves advection of high and low chlorophyll signals around an eddy, resulting in a low which is offset from the center of an anticyclone and a high which is offset from the center of a cyclone.Ekman pumping would be expected to produce negative anomalies in cyclones with a positive "halo" and positive anomalies in anticyclones with a negative "halo" Gaube et al. (2014), Fig. 2).Trapping of chlorophyll involves eddies retaining the properties that they had when shed from a boundary current, which would generally imply low values in anticyclones and high values in cyclones.Eddy intensification would be expected to produce the same picture, as cyclones would see rising nutriclines in the center but anticyclones would see deepening nutriclines.The basic picture seen in the Arabian Sea is inconsistent with the first two mechanisms but is potentially consistent with the second two.However, without in situ data it is impossible to validate either of these mechanisms.

Temporal variability
Time series of chlorophyll a, phosphate and nitrate for all GFDL models are shown in Fig. 7a to c within the whole region and compared against the corresponding GSM5 satellite results or WOA09.Note that the 8 years of the model output, selected as the last 8 years of the run, would not be expected to correspond to the 8 actual years in the satellite data.The annual cycles of chlorophyll a and biomass are quite similar to each other in all GFDL models, insofar as they show two distinct blooms in yearly cycle.The maximum values that can be considered as a winter bloom in the whole region are mostly seen around February (Piontkovski et al., 2011), with values of 0.32-0.38,0.48-0.62,1.0-2.0,1.5-2.2,0.8-1.6, and 0.6-0.75mg m −3 for chlorophyll in CORE-TOPAZ, Coupled-TOPAZ, CM2.6 (miniBLING), miniBLING (Low resolution), BLING and satellite data, respectively.A summer bloom is also pronounced in September as a second maximum in the yearly cycle over the whole region, with peak magnitudes of about 0.25-0.52,0.65-0.7,0.65-1.15,0.8-1.15,0.5-0.75, and 0.75-1.3mg m −3 for chlorophyll across the different data sets.
Notice the results from the BLING model run in the coarser resolution ESM2M code (purple lines).The differences between BLING and miniBLING (light blue lines) in this code are just due to having fixed iron and no dissolved phosphorus in miniBLING.The light field in these ESM2M runs is computed from using TOPAZ-derived chlorophyll, so that all three models see identical physical conditions.Both BLING and miniBLING in ESM2M produce an asymmetry in chlorophyll between February and September that is similar to that produced in CM2.6 miniBLING.This asymmetry is not seen in TOPAZ.Analysis of what drives this asymmetry shows that it is not straightforward.All of the model runs show an asymmetry in the nutrient concentrations that is in the opposite direction as the observations, with higher nutrients in February than in September, as shown in Fig. 7b.As we will show later in the manuscript, this is probably associated with the models mixing to excessive depth during the wintertime.However, in TOPAZ this does not produce an asymmetry in chlorophyll, while in BLING and miniBLING it does.There are two possible reasons for this: 1.The equilibrium assumption, which means that biomass in both BLING and miniBLING is not directly simulated.In TOPAZ, the growth of plankton during the spring is limited by the biomass of phytoplankton, whereas in the fall TOPAZ continues to have higher heterotrophic biomass (diagnosed from growth rates over previous months) that then grazes the plankton.
In BLING and miniBLING, by contrast, the biomass responds almost instantaneously to changes in growth conditions.This means that if the growth rate increases from 0.05 to 0.1 day −1 over the course of a month, the biomass associated with large phytoplankton will increase eightfold, even though the additional growth should only be enough to give an increase of a factor of 30 days × 0.05 day −1 = 1.5.Possibilities for addressing this effect include replacing the DOP tracer with a biomass tracer, which could then be partitioned between the different phytoplankton boxes based on the temporally smoothed growth rate, or increasing the timescale over which the growth rate is smoothed when biomass is calculated.light, whereas in BLING it is calculated using the mixed layer average light.Preliminary results with a very coarse resolution model using BLING show that this reduces the summer-winter asymmetry slightly, but is not sufficient to make the February bloom smaller than the September bloom.This effect will also be addressed in future research.
It is likely that all three of these factors -too deep winter mixed layers leading to too high nutrients, too little light limitation and instantaneous response to changes in growth conditions, are all responsible for the overly strong blooms in boreal winter in the Arabian Sea.
To get a better sense of the mechanisms driving the blooms in the model, the biomass (mol P kg −1 ) of the miniBLING CM2.6 model is compared with the light intensity in the mixed layer and the light-saturated photosynthesis rate (carbon specific) (s −1 ) in Fig. 8a and b for January of year 195.The two terms in Fig. 8 are the two terms in the model that affect growth rate.Because biomass in the miniBLING model is a function of growth rate only, understanding the variation in two terms is sufficient to understand what drives the variation of biomass in the model.The biomass production and mixed layer light intensity (Fig. 8a) are not meaningfully correlated parameters.On the other hand, the biomass and the light-saturated carbon specific growth rate (Fig. 8b; indicating the degree of nutrient limitation) are positively correlated.From this, it can be concluded that the blooms in this region are more driven by nutrient rather than light, consistent with, for example, Gomes et al. (2008).This suggests in turn that it is likely biases in nutrient supply that drive biases in productivity.
We can get more insight into nutrient biases by examining the individual tendency terms associated with advection, vertical diffusion and subgridscale eddy fluxes and time rate of change of nutrients.For simplicity, in this paper we combine the vertical diffusive flux associated with small-scale mixing from the background diffusion with that due to the mixed layer parameterization.Figure 9 shows PO 4 advection, diffusion and time tendency flux terms for the whole region (56-66 • E, 15-26 • N) over a typical year.We calculate these by integrating the time tendency terms for phosphate over the top 50 m.The results show that the dominant source in whole region during the winter bloom is diffusion, suggesting the model predicts excessively strong mixing during the wintertime.By contrast, the advection dominates diffusion during the summer bloom, supplying the majority of nutrients during the months of July and August.The fact that the summertime bloom is close to observations suggests that the model correctly simulates this wind-driven upwelling.In addition to having annual cycles that are different from observations, the models also differ from data in terms of interannual variability.As shown in Fig. 10, low-resolution models (CORE-and coupled-TOPAZ) provide an almost uniform seasonal coefficient of variation (mean C.o.Vs are 0.15 and 0.18, respectively), while both data and eddy resolving CM2.6 models show higher interannual variability and seasonal changes (mean C.o.Vs are 0.35 and 0.5, respectively).The C.o.Vs are particularly higher during the winter and summer blooms in the observations, while the low-resolution models do not see these signals.In other words, the low-resolution models fail to get enough variability, while the high-resolution models produce too much interannual variability.Together with Fig. 4, this statistical analysis suggests that eddies are necessary to explain the variability in the data as opposed to the low-resolution models, but that the high-resolution model does not properly capture this variability.Below, we examine the relationship of eddies and blooms in the high-resolution models.4.2 Blooms and sea surface height in CM2.6

Large-scale correlation
The relationship between SSHA and chlorophyll is quite different in the model as compared to the satellite.Monthly variation in the cross-correlation of chlorophyll and SSHA for 8 consequent years in CM2.6 is shown in Fig. 11.As in the remote sensing, the model shows annual cycles of variation in the cross-correlation, suggesting a repeatable yearly phenomenon in the region.However the structure of this annual cycle is not consistent with the satellite data.The model predicts several months (i.e.March-August) with anti-correlation for most of the years, but with values less than 0.5, smaller than the peak anti-correlation values in satellite results.The model also predicts that several other months (i.e.October-February) should have no or even positive correlation, while the satellite shows strong negative correlations during these months.

Blooms in mesoscale structures
Why does the GFDL CM2.6 model not produce the same relationship between SSHA and chlorophyll as the satellite?We can gain some insight by examining snapshots of the two fields.In Fig. 12a and b, sea surface chlorophyll a concentration and sea surface height anomaly (SSHA) are shown at two snapshots of time, 9 November and 28 December for model year 195.Comparing the figures with the corresponding satellite results in Fig. 3 for the months of November and December, we see that the southern part of the GFDL model is more similar to the satellite data, with high concentrations of chlorophyll a tending to be located at the center of cyclones.In contrast, in the northern part of the region, the GFDL model predicts high chlorophyll at the edges of the cyclones as well as in the center of anticyclones.The eddy structures have smaller diameters in GFDL results than the field observations, though it is not clear whether this rep-resents smoothing in the AVISO product or some physical weakness of the model.
We now focus on the few examples in our model output where chlorophyll blooms are found in the center of cyclonic eddies.These are denoted as E1 and E2 in Fig. 12a and b.To track the movement of the selected eddies, E1 and E2, over the time from 9 November to 28 December, modeled chlorophyll and SSHA are shown in Fig. 12c and d along two different latitudes, 16 • N (for E1) and 19 • N (for E2). Figure 12c shows that E1 moved westward during this period of time, and that the chlorophyll concentration was kept high within the central part of the eddy.E1 appears to be created by the passage of a cyclone, similar to the eddy observed by Wang and Zhao (2008) in the aftermath of Cyclone Gonu.Similarly, as shown in Fig. 12d, E2 was a persistent eddy with both central and edge blooms during the month of November that started to move towards the west during December along 19 • N.However, at other latitudes, the largest blooms offshore are found along gradients in SSH rather than being associated with maxima or minima.This suggests a different mechanism for producing blooms in the model.Following Gaube et al. (2014) it appears that the eddy stirring mechanism is dominant.Satellite data (i.e.see Fig. 3 for the month of May) provide some hints of high-chlorophyll plumes being advected away from coastal regions.As shown in Fig. 12a and b, high velocities in the marginal region between adjacent cyclonic and anticyclonic eddies can cause such plumes in the GFDL models as well.
Why is the model only able to simulate the relationship between SSH and chlorophyll in the southern part of the domain?We hypothesize this is due to differences in stratification between the two regions.The average water temperature (colors) and the macronutrient (  region are shown in Fig. 13.In the northern part of the region (see Fig. 13a and b), the GFDL model provides a reasonably good estimation of the mean temperature field near the surface, but subsurface temperatures are not as consistent as there is far too little stratification.This is also associated with a very weak nutricline in CM2.6.Variations in isopycnal depth will therefore not lead to big differences in nutrient supply.Figure 13c and d show the same fields for the southern part of the region.Unlike the northern part of the domain, the temperature gradient over these depths is well estimated by CM2.6.While the nutricline is still too weak there is some gradient in nutrients between 80 and 120 m.
As seen in Fig. 1d both the ARGO and WOA09 wintertime mixed layer depth is considerably deeper than the summertime mixed layer depth, reaching a maximum of 65 m.However, in the northern regions of the model the MLD seems to be too deep in winter, reaching values of 130-150 m.This suggests that the overly deep mixed layer in the northern part of the region may explain both the tendency towards an overly strong winter bloom and the failure of mesoscale eddies in modulating chlorophyll blooms.If we look during the time period where we have eddies E1 and E2 (November-December year 197, Fig. 6c, d) we see shallower mixed layers associated with both eddies.Both the temperature and mixed layer biases in the northern part of the Arabian Sea may result from having too much water from the Persian Gulf in this region.This can be seen in the yearly averaged subsurface salinity-density distribution over the region, shown in Fig. 13e and f for both WOA09 data and CM2.6 (model year 197), respectively.Figure 13e shows two separate tongues of salty water, one near the surface and one at the depth of ∼ 300 m.These salty water signals are consistent with the seasonal cycle of Persian Gulf outflow as discussed in Ezam et al. (2010).On the other hand, CM2.6 shows one subsurface salty water signal from the northern part, which is deep and strong enough to result in weak stratification in the north to a depth of 250 m, as shown in Fig. 13f.These results suggest that a sharp thermocline and nutricline is necessary for eddy activity to modulate the mixing of nutrients to the surface.We test the idea that a sharper thermocline could modulate mixing of nutrients to the surface by looking at the sources of nutrient in the southern part of Arabian Sea where eddybloom relationships similar to observations are occasionally seen.Accordingly, the region containing eddy E1 in Fig. 12 is analyzed to determine the physical mechanisms by which nu-trient is transported into the surface layer.Figure 14  in 197 it is larger and positive (∼ 10 mol m −2 month −1 ) in the eddy while the advective flux is actually negative in this region.By contrast in Year 198, there is no cyclonic eddy and the diffusive fluxes are much smaller.
The bloom associated with eddies E1 and E2 do not fit with any of the mechanisms highlighted in Gaube et al. (2014).We first consider the mechanism of trapping.Eddy E1 is generated in the ocean interior, not as a result of coastal upwelling.As shown in Fig. 15, the nutrient supply rate ranges between 5 and 8 mmol m −2 month −1 in the eddy.The concentrations in this eddy are only 0.01 µM (5 mmol m −2 ) over the top 50 m.It cannot be the case that the nutrients in the eddy can last for several months as a result of "trapping", there must be a continuous supply.Moreover although eddy E2 shows a horizontal advection signal in November (with a positive ring around the edge in Fig. 12a), the signal in December has the opposite sign.Eddy intensification is also an unlikely mechanism for explaining the blooms, as dSSH/dt is relatively small (particularly if we track the minimum SSH associated with E1 in Fig. 12c or E2 in Fig. 12d).Finally, Ekman pumping signatures in Gaube et al. (2014) have the opposite sign as what is seen in E1 and E2.
Our results also contrast with those in Resplandy et al. (2011).The focus in Resplandy et al. (2011) is on the productivity driven by horizontal and vertical advection in summer and mostly vertical advection in winter.This contradicts our finding of a primary diffusive source of nutrient in winter although it is consistent with the finding of advective source of nutrients in summer.We point out that in our model, the only two eddies that actually look like what we see in the satellite observations involve enhanced mixing from below.This is a different result from Levy et al. (2014) and Resplandy et al. (2011).Moreover it is not clear whether these papers get the seasonal correlation with SSH or not.Resplandy et al. (2011) do not focus on structures at the eddy scale as they are more concerned with the net impact of eddies.
To summarize, we hypothesize that 1.The reason that blooms are found in cyclones in the Arabian Sea during the NEM is that the dominant source of nutrients to the surface, i.e. mixing (Barimalala et al., 2013;Kawamiya and Oschlies, 2003) is concentrated there.
2. Interannual variability in wintertime blooms in the Northwest Arabian Sea is controlled by the combined presence of these eddies and strength of wintertime cooling.
3. Excessive mixing (resulting in too weak a thermocline) prevents mixing from being modulated by eddies in the model except occasionally in the southern part of our region.In the real world the modulation of mixing seen in Fig. 14 extends into the Northwest Arabian Sea and the Gulf of Oman.

Conclusions
Our analysis of bloom variability in the northwestern Arabian Sea and Gulf of Oman has illustrated both similar and dissimilar descriptive features between satellites and a suite of models and explored the various mechanisms involved.Satellite analyses demonstrate the existence of two blooms, the stronger one associated with the Southwest Monsoon and the weaker one associated with the Northeast Monsoon as also shown by Madhupratap et al. (1996); Kawamiya and Oschlies (2003); Murtugudde et al. (2007);and Al-Azri et al. (2010).We demonstrate a pronounced anti-correlation between SSHA and chlorophyll blooms during certain times in northern winter but a much weaker relationship in other months (typically northern summer) with the relationship disappearing as the blooms vanish in the months of April and May (northern spring).While the depth of thermocline and nutricline and also the stratification affect the convection during the Northeast Monsoon (Dickey et al., 1998;Kumar et al., 2001;Wiggert et al., 2002), we show that a thin nutricline and/or thermocline and a strong stratification are also

S. Sedigh Marvasti et al.: Challenges in modeling spatiotemporally varying phytoplankton blooms
required to enable cold eddies to bring nutrients to euphotic zone and develop phytoplankton blooms.During the wintertime monsoon, while both cooling in the winter and eddies control the blooms, variability in bloom location will arise from variability in the location of eddies, and so may not be predictable.In contrast, during the Southwest Monsoon the dominant upwelling associated with the intense environmental forcing supersedes the effect of eddies and the activity of the cold eddies is not pronounced.Understanding of this phenomenon has been sought using five different 3-D ocean-atmosphere models, including a CORE-forced ocean with the TOPAZ biogeochemistry, a coupled model with the TOPAZ biogeochemistry and CM2.6.Because the coarse models with TOPAZ are not able to capture eddies and the interannual variability, CM2.6 (miniBLING), an eddy-resolving high-resolution model, was also considered for simulating the spatial and temporal changes of the bloom in the region.This model simulates the two blooms seen in the data and shows that the nutrients driving the northern summer bloom are supplied by advection while those driving the wintertime bloom are supplied by vertical diffusion.However, this model is unable to simulate the seasonal relationship observed in the satellite products between blooms and sea surface height.Although there is some anti-correlation, it tends to be associated with larger spatial scales and not really related to eddies.Instead, eddies in the model usually wrap the chlorophyll around themselves, producing high chlorophyll concentrations around their edges and not at their centers.Comparing the model results to field measurements (WOA09) showed that the model does not account for the strong thermocline and nutricline in the northern part of the region.In the wintertime, this leads to excessive convective supply of nutrients and too strong of a bloom.However, for a few cases, eddies with blooms at the center are tracked in the southern part of the domain.In this region, consistency is observed between the model results and the field data.Analysis of the term balances in mixed layer show that eddies in this region modulate the diffusive supply of nutrients.We suggest that what happens in the model in the Southern Arabian Sea actually describes the Arabian Sea as a whole according to the observations and the field data.The model misses the eddy signal in the north because it lacks a thin nutricline, motions of which will lead to differences in nutrient supply.In the real world, eddies modulate the diffusive supply of nutrients during the wintertime and there is more mixing in the eddy centers along with the diffusive supply provided by the cooling in the wintertime.
Accordingly, there is a potential to improve the numerical models by better simulating the Persian Gulf Outflow to produce a sharper thermocline, allowing more realistic nutrient supply.Overflows are difficult to simulate in level-coordinate models because they are prone to excessive entrainment of the dense plume (Winton et al., 1998).While significant effort has gone into simulating the Denmark Straits overflow at coarse resolution (Legg et al., 2009), our results show that smaller overflows such as the Persian Gulf may be regionally significant.This may provide further impetus for using isopycnal models in high resolution simulations, as such models can potentially simulate such overflows with greater fidelity.
It is worth noting that regional models, (such as Resplandy et al., 2011) do have the potential to better simulate the hydrography of the Northern Arabian Sea.Because such models are very tightly constrained through "sponges" that restore hydrography at the boundaries, they may not have the problems that global models do at representing the effects of overflows that they do not properly simulate.However, such models cannot by themselves simulate the effects of changing climate, which in turn changes the boundary conditions.For this reason, global models must still be used for projection, making it important to identify the reasons that they are not going to work.

Figure 3 .
Figure 3. Satellite chlorophyll a in mg m −3 (colors) and sea-surface height anomaly (SSHA, contours) in cm (contour interval = 5 cm) in the Gulf of Oman and Northwest Arabian Sea over the course of 2001.

Figure 4 .
Figure 4. Chlorophyll a in mg m −3 (colors) and sea surface height anomaly (SSHA, contours) in cm (contour interval = 5 cm) in the Gulf of Oman and Northwest Arabian Sea during November in different years.

Figure 8 .
Figure 8. Modeled biomass in CM2.6 in P units (mol P kg −1 ) vs. (a) Mixed layer irradiance (Wm −2 ); (b) Light-Saturated photosynthesis rate (carbon specific) (s −1 ) 56-66 • E, 15-26 • N for January of year 195.In the model, biomass is a function of the growth rate smoothed over several days, and the light-saturated photosynthesis rate indicates the extent to which this growth rate is controlled by nutrient limitation.

Figure 12 .
Figure 12.CM2.6 (miniBLING) Surface chlorophyll a concentration and sea surface height anomaly (SSHA) November and December during a year where the observed eddy-bloom interaction is seen in the southern part of the Arabian Sea.

Figure 13 .
Figure 13.(a-d) Seawater temperature (colors, • C) and phosphate (PO 4 ) concentration (contours, µM) for the northern (top row) and southern (middle row) parts of the central Arabian Sea.(e-f) yearly averaged subsurface distribution of salinity (black contours) and potential density (red contours).Left-hand column shows observations with the Persian Gulf plume centered at 300 m, right-hand column results from CM2.6 model with a much broader mixing of salinity over the top 200 m.

Figure 14 .
Figure 14.Surface chlorophyll in mg m −3 .(a-b) Advective flux of phosphate to top 50 m in mol m −2 (c-d, colors), and diffusive flux of phosphate in mol m −2 (e-f, colors) with sea surface height (contours in m, contour interval 0.02 m) for eddy E1 (63-66 • E, 15-18 • N) for the month of December during the two CM2.6 model years 197 and 198.