Challenges in modelling spatiotemporally varying phytoplankton blooms in the 1 Northwestern Arabian Sea and Gulf of Oman 2 3

17 Recent years have shown an increase in harmful algal blooms in the Northwest Arabian Sea 18 and Gulf of Oman, raising the question of whether climate change will accelerate this trend. 19 This has led us to examine whether the Earth System Models used to simulate phytoplankton 20 productivity accurately capture bloom dynamics in this regionboth in terms of the annual 21 cycle and interannual variability. Satellite data (SeaWIFS ocean color) shows two 22 climatological blooms in this region, a wintertime bloom peaking in February and a 23 summertime bloom peaking in September. On a regional scale, interannual variability of the 24 wintertime bloom is dominated by cyclonic eddies which vary in location from one year to 25 another. Two coarse (1°) models with the relatively complex biogeochemistry (TOPAZ) 26 capture the annual cycle but neither eddies nor the interannual variability. An eddy-resolving 27 model (GFDL CM2.6) with a simpler biogeochemistry (miniBLING) displays larger 28 interannual variability, but overestimates the wintertime bloom and captures eddy-bloom 29

climatological blooms in this region, a wintertime bloom peaking in February and a 23 summertime bloom peaking in September. On a regional scale, interannual variability of the 24 wintertime bloom is dominated by cyclonic eddies which vary in location from one year to 25 another. Two coarse (1°) models with the relatively complex biogeochemistry (TOPAZ) 26 capture the annual cycle but neither eddies nor the interannual variability. An eddy-resolving 27 model (GFDL CM2.6) with a simpler biogeochemistry (miniBLING) displays larger 28 interannual variability, but overestimates the wintertime bloom and captures eddy-bloom 29 coupling in the south but not in the north. The models fail to capture both the magnitude of 30 the wintertime bloom and its modulation by eddies in part because of their failure to capture 31 the observed sharp thermocline/nutricline in this region. When CM2.6 is able to capture such 32 features in the Southern part of the basin, eddies modulate diffusive nutrient supply to the 33 surface (a mechanism not previously emphasized in the literature). For the model to simulate 34 the observed wintertime blooms within cyclones, it will be necessary to represent this 35 relatively unusual nutrient structure as well as the cyclonic eddies. This is a challenge in the 36 Northern Arabian Sea as it requires capturing the details of the outflow from the Persian  something that is poorly done in global models. a highly productive region (Madhupratap et al., 1996;Tang et al., 2002), with satellite 42 estimates of carbon export of 137 gC/m 2 /yr, much higher than the ~80 gC/m 2 /yr found in the 43 Subpolar North Atlantic and Pacific (Dunne et al., 2007). Peak chlorophyll-a concentrations 44 exceed 0.7 mg/m 3 in this region (Fig. 1a). 45 This region may be changing in important ways. In both the Persian Gulf and the Gulf 46 of Oman, there is evidence that harmful algal bloom (HABs) and their impacts are increasing 47 (Richlen et al., 2010). HAB occurrences have been more frequently reported in the Gulf of 48 Oman rather than Persian Gulf. A total of 66 red tide events (mostly dominated by Noctiluca 49 scintillans) have been recorded between 1976 and 2004 including 25 blooms resulting in mass 50 mortality of fish and marine organisms. Reasons for the increase in blooms include 51 aquaculture activities and industrial and sewage inputs, natural dispersal and human-aided 52 transport, long-term increases in nutrient loading and global expansion of species (Richlen et 53 al., 2010) as well as global climate change (Goes et al., 2005). The latter paper suggested that 54 increasing blooms were driven by an increase in the strength of the Asian monsoon. 55 Evaluating such a possibility and extending it into the future requires the use of Earth 56 System Models. However, such projections will only be as good as the models on which they 57 are based. In this paper we examine several models run at the Geophysical Fluid Dynamics 58 Laboratory in the Arabian Sea. We consider numerical results from four different 3D global 59 Earth system models, which we denote CORE-TOPAZ, Coupled-TOPAZ, Coupled-BLING 60 and the Geophysical Fluid Dynamics Laboratory Climate Model version 2.6 (CM2.6) 61 eddies are associated with blooms, but only during the winter. This means that interannual 126 variability in blooms will be shaped by mesoscale eddy activity and may not be predictable. 127 Results of the 3D global Earth system models are discussed in section 4. Annual cycles of 128 variation of chlorophyll-a and nutrients for all GFDL models within the whole region are 129 compared against the corresponding satellite results and field measurements. The models tend 130 to overestimate wintertime productivity, in large part due to excessive mixing. They also fail 131 to explain the bloom-SSHA relationship except in a few special cases. We argue that the 132 eddies act to modulate turbulent mixing of nutrients to the surface-a mechanism not 133 emphasized in previous literature. However, this can only occur if there is a strong and 134 relatively shallow nutricline. Since the model only simulates such a feature in the Southern 135 Arabian Sea, it does not capture the observed relationship between SSH and biology. Both the 136 overestimation of the wintertime bloom and the failure to predict its modulation by eddies can 137 thus be traced to difficulties in modeling the stratification of the Northwest Arabian Sea, most 138 likely as a result of a failure to properly simulate overflows. 139 140 2 Description of data and models 141

Satellite products 142
We examine the relationship of blooms and eddies using the GSM5 Maritorena et al. (2002) 143 product based on the SeaWIFS (Sea-viewing Wide Field-of-view Sensor) ocean color data 144 and Sea Surface Height Anomaly (SSHA), based on altimeter data acquired from the 145 Archiving, Validation and Interpretation of Satellite Oceanographic (AVISO) Data Center 146 (http://www.aviso.oceanobs.com). The SSH anomaly is calculated relative to the annual 147 cycle. 148 The GSM algorithm represents the normalized water leaving radiance L wN (λ) at multiple 149 wavelengths as a nonlinear function, as following (Maritorena et al., 2002), 150 (1) 151 where t is the sea-air transmission factor, F 0 (λ) is the extraterrestrial solar irradiance, n w is the 152 index of refraction of the water, seawater backscatter b bw (λ), absorption a w (λ), a ph * is the 153 chlorophyll-a (chl) specific absorption coefficient, S is the spectral decay constant for 154 absorption by chromophoric dissolved organic materials (CDOM), η is the power-law 155 exponent for the particulate backscattering coefficient, and λ 0 is a scaling wavelength (443 156 nm). The cdm absorption coefficient [a cdm (λ 0 )], and slope factor S then determine the 157 absorption across a range of wavelengths while the particulate backscatter coefficient [b bp (λ 0 )] 158 and coefficient η constrain the scattering. Letting λ 0 be 443 nm assuming that all terms other 159 than chl, [a cdm (λ 0 )] and b bp (443nm) are constant, one can then use the water leaving radiance 160 to invert for chl, a cdm , and backscatter. One limitation of this approach is that if the inherent 161 optical properties vary with time, this variation will introduce errors into the estimate. 162 Following Behrenfeld et al. (2005), we convert the backscatter coefficient into units of 163 particulate carbon biomass using the relationship p carb =13000(b bp -0.00035). 164 Satellite-based remote sensing is the only observational method suitable for measuring 165 physical and biological properties over large regions of the ocean. However, satellite ocean 166 color and SST are limited to surface distributions and provide no information about the 167 vertical structure within the ocean (McGillicuddy et al., 2001). Additionally acquiring data 168 requires cloud-free viewing of the ocean surface, which as we will see is a problem in this 169 region at certain times of the year. This lack of information motivates our examination of 170 numerical models, which ideally could be used to provide estimates of the ocean state when 171 observations are sparse as well as to extrapolate both vertically and into the future. 172 173

Numerical models 174
Numerical results are presented in this paper based on the output of four different 3D 175 global Earth system models, which we denote CORE-TOPAZ, Coupled-TOPAZ, Coupled-176 BLING/miniBLING and GFDL CM2.6 (miniBLING). The first two of these models use the 177 relatively complex TOPAZ biogeochemistry, but have low resolution and do not resolve 178 eddies. The third uses two simplified biogeochemistry codes (BLING and miniBLING) which 179 do not carry phytoplankton biomass as a separate variable while the last model has very high 180 resolution and uses the miniBLING simplified biogeochemistry. Below, we describe the 181 different physical models, followed by a summary of the biogeochemical codes run within 182 these models. The ocean-ice model used in the CORE-TOPAZ model follows the corresponding 187 components of the GFDL CM2.1 global coupled climate model (Delworth et al., 2006). The 188 vertical resolution ranges from 10 m over the top 200 m to a maximum thickness of 250 m at 189 5500 m depth with 50 layers in all. The meridional resolution is 1 o , whereas the zonal 190 resolution varies between 1 o in mid-latitudes and 1 / 3 o at the equator. North of 65 o , a tripolar 191 grid is employed to avoid singularity arising from convergence of meridians at the North 192 Pole. Up-to-date parameterizations of mixed-layer dynamics, isopycnal mixing, advection by 193 subgridscale eddies, bottom topography, bottom flows, and lateral viscosity are included-for 194 more detail see Griffies et al. (2005) and Gnanadesikan et al. (2006). Both the dynamics and 195 thermodynamics sea ice are simulated with of five thickness classes of sea ice being resolved. 196 In the CORE-TOPAZ model, surface forcing is set using the Coordinated Ocean-ice 197 Reference Experiment (CORE) protocol (Griffies et al., 2009) forcing the CORE runs could be thought of as "best guess" observationally based estimates. 204 Such a prescription omits important feedbacks whereby the atmosphere ensures that rainfall 205 and evaporation are consistent with each other, although the restoring correction is a crude 206 representation of these feedbacks. We use the version of the model described in Gnanadesikan 207 et al. (2011), which analyzed different modes of interannual variability in biological cycling 208 across the Pacific Ocean. 209 The Coupled-TOPAZ model corresponds to the control simulation of the GFDL 210 ESM2M submitted as part of the IPCC AR5 process (Dunne et al., 2012). In this model the 211 ocean is coupled to the atmosphere, land, and sea ice components. Gnanadesikan et al. (2014) 212 discuss the behavior of this model in the North Atlantic, but its behavior in the Arabian Sea 213 has not been previously analyzed. Two additional version of this model, referred to here as 214 Coupled-BLING, was run using the BLING and mini-BLING biogeochemical models 215 described below, but with the light field given by the TOPAZ code. The difference between 216 the 1-degree models emphasizes differences due to biological formulation. 217 The ocean component of ESM2M employs the MOM4p1 code of Griffies et al. (2009) 218 which largely mimics the CM2.1 ocean (identical horizontal and vertical resolution and 219 parameterization of mixing). However, ESM2M ocean uses a rescaled geopotential vertical 220 coordinate (z*; Adcroft et al., 2004;Stacey et al., 1995) for a more robust treatment of free 221 surface undulations. The ESM2M implementation includes updates to the K-profile 222 parameterization (Large et al., 1994) based on Danabasoglu et al. (2006), as well as model-223 predicted chlorophyll modulation of short-wave radiation penetration through the water 224 column. ESM2M also includes completely novel parameterizations relative to CM2.1, such as 225 parameterization of submesoscale eddy-induced mixed layer restratification (Fox-Kemper et 226 al., 2008). Instead of prescribed vertical diffusivity for interior mixing (Bryan and Lewis, 227 1979), ESM2M employs the Simmons et al. (2004) scheme along with a background 228 diffusivity of 1.0 x 10 -5 m 2 s -1 in the tropics and 1.5 x 10 -5 m 2 s -1 poleward of 30° latitude 229 following a tanh curve. 230 The Geophysical Fluid Dynamics Laboratory Climate Model version 2.6 (CM2.6) is a 231 high-resolution eddy-resolving model. This model has the same atmosphere model and ocean 232 Physics as CM2.5 (Delworth et al., 2012). CM2.6's ocean component has higher horizontal 233 resolution than CM2.5, with grid spacing, which is changeable from 11 km at the equator to 234 less than 4 km at very high latitudes. This means that the model is capable of resolving eddy 235 features in the tropics, as we will see below. 236 237 2.2.2 Biogeochemical Cycling codes 238 The TOPAZ code (Tracers of Ocean Productivity with Allometric Zooplankton code 239 of Dunne et al., 2010), keeps track of five inorganic nutrients used by phytoplankton: nitrate 240 and ammonia, inorganic phosphate, silicate, and dissolved iron. Additionally, the model 241 carries three other dissolved inorganic tracers: dissolved inorganic carbon, alkalinity and 242 dissolved oxygen. Based on the work of Dunne et al. (2007), the model also keeps track of 243 fine lithogenic material, which plays a role in ballasting organic material and delivering it to 244 the sediment (Armstrong et al., 2002;Klaas and Archer, 2002). The five inorganic nutrients 245 are taken up in different ways by three classes of phytoplankton: small, large and 246 diazotrophic. A comprehensive description of TOPAZ v2 can be found in the supplemental 247 material of Dunne et al. (2013). 248 TOPAZ is unusual among comprehensive Earth System Models in that it uses a highly 249 parameterized version of grazing. Instead of grazers being explicitly simulated, grazing rates 250 are simply taken as a function of phytoplankton biomass, with different power-law 251 dependence for small and large phytoplankton. The grazing formulation was fit to about 40 252 field sites to produce a size structure that transitions realistically between being dominated by 253 small phytoplankton and low particle export ratio at low levels of growth and large 254 phytoplankton and high particle export ratio in nutrient and light-replete conditions. At 255 equilibrium, the resulting parameterization produces biomass that is a function of growth rate 256 (linear for small plankton, cubic for large). A similar scaling in particle size spectrum was 257 seen across ecosystems by Kostadinov et al. (2009). In contrast to models that explicitly 258 simulate zooplankton, it does not depend poorly known zooplankton behavioral parameters 259 (such as handling efficiency or grazing half-saturation) on to the details of how different 260 trophic levels interact. 261 Even though it does not simulate zooplankton explicitly, TOPAZ still carries over two 262 dozen tracers, making it extremely expensive to run in high-resolution simulations. For this 263 reason Galbraith et al. (2010)

developed the Biogeochemistry with Light Iron Nutrients and 264
Gasses (BLING) model (Galbraith et al., 2010), which parameterizes the entire ecosystem. 265 The original version of BLING has only five explicit tracers: dissolved inorganic phosphorus 266 (PO4), dissolved organic phosphorus (DOP), dissolved Iron (Fe), DIC (dissolved inorganic 267 carbon), and oxygen (O 2 ). It includes the impacts of macronutrient and micronutrient 268 limitation and light limitation on phytoplankton by using these to calculate a growth rate. 269 Using the same machinery as TOPAZ, it then uses this growth rate and implicit treatment of 270 community structure to estimate phytoplankton biomass, and uses this biomass to calculate 271 the rate at which nutrient is taken up by plankton and cycled through the ecosystem. 272 The miniBLING code (Galbraith et al., 2015) represents a further simplification. In 273 this model the iron field is taken from a lower-resolution version of the model (an 274 approximation which has limited impact in the Arabian Sea, where phytoplankton are 275 generally not iron-limited) and so Fe is not treated prognostically. Additionally the DOP pool 276 is eliminated. Simulations using the ESM2M physical model show that control simulations of 277 oxygen and surface nutrients produced by miniBLING and BLING models are very similar to 278 those produced in the same model with TOPAZ (Galbraith et al., 2015). This manuscript also 279 shows that BLING and miniBLING simulate very similar patterns of oxygen change and 280 anthropogenic uptake in a simulation where CO2 is increased by 1% per year until it doubles. 281 It should be noted that simplified BLING and miniBLING codes neglect some 282 processes that may be important. Only nonliving components are advected and mixed by the 283 ocean circulation, which could result in inaccurate distribution of biology in frontal regions at 284 high resolution. Also, the rich behavior of the nitrogen cycle with its interaction with iron, 285 phosphorus and oxygen cannot be simulated with one macronutrient tracer Behrenfeld (2010). 286 Specifying iron limitation, as done in miniBLING, may also have some impacts in our region. winter bloom is also pronounced in February as a second maximum in a yearly cycle, where 311 the magnitudes are about 0.07 mg/m 3 , 40 mgC/m 3 , and 0.07 mgC/m 3 for chlorophyll, 312 particulate carbon and CDOM, respectively, within the whole region, and about 0.09~1.5 313 mg/m 3 , 55~80 mgC/m 3 , and 0.11~0.14 mgC/m 3 for chlorophyll, particulate carbon and 314 CDOM, respectively, within the smaller region. That the summer bloom in the both regions is 315 stronger than the winter bloom has been discussed by Al-Azri et al. (2010), and Levy et al. 316 (2007). 317 Hadd. This difference in bloom location is perfectly reflected in the different patterns of the 336 eddies. 337

Chlorophyll-Sea Surface Height Anomaly (SSHA) cross-correlation 338
The seasonal relationship between chlorophyll and SSHA can be seen in monthly variation of 339 the spatial cross-correlation between the two variables over the entire northwest Arabian Sea. positive or very small, which suggests no relation between the mesoscale eddies, and the 357 blooms. As discussed in Kumar et al. (2001), low primary production is observed after 358 termination of winter cooling during Spring Inter-Mansoon (SIM) (Gomes et al., 2008). This 359 result would be also consistent with SIM producing weak atmospheric forcing in the region. 360 In contrast, after the summer bloom (typically October-December) as the average values of 361 chlorophyll-a decrease, chlorophyll and SSHA become relatively highly anti-correlated. The 362 reason for the anti-correlation is the persistence of chlorophyll at the regions with negative 363 SSHA that typically considered to be cyclonic (cold) eddies and disappearance of 364 chlorophyll-a in positive SSHA that assumed to be anti-cyclonic (warm) eddies. Particle 365 backscatter also provides almost same cross-correlation and average value results suggesting 366 that the chlorophyll-a signal does not result purely from photo-adaptation. Moreover, the 367 CDOM-SSHA cross-correlation shows the same patterns for bloom times but typically shows 368 lower correlation. 369 The spatial relationship between blooms and eddies seen in the Northern Arabian Sea can be 370 compared with the patterns noted by Gaube et al. (2014). Their eddy stirring mechanism 371 involves advection of high and low chlorophyll signals around an eddy, resulting in a low 372 which is offset from the center of an anticyclone and a high which is offset from the center of 373 a cyclone. Ekman pumping would be expected to produce negative anomalies in cyclones 374 with a positive "halo" and positive anomalies in anticyclones with a negative "halo" (Gaube 375 et al. (2014) , Fig. 2). Trapping of chlorophyll involves eddies retaining the properties that 376 they had when shed from a boundary current, which would generally imply low values in 377 anticyclones and high values in cyclones. Eddy intensification would be expected to produce 378 the same picture, as cyclones would see rising nutriclines in the center but anticyclones would 379 see deepening nutriclines. The basic picture seen in the Arabian Sea is inconsistent with the 380 first two mechanisms but is potentially consistent with the second two. However, without in-381 situ data it is impossible to validate either of these mechanisms. Notice the results from the BLING model run in the coarser resolution ESM2M code (purple 399 lines). The differences between BLING and miniBLING (light blue lines) in this code are just 400 due to having fixed iron in miniBLING. The light field in these ESM2M runs is computed 401 from using TOPAZ-derived chlorophyll, so that all three models see identical physical 402 conditions. Both BLING and miniBLING in ESM2M produce an asymmetry in chlorophyll 403 between February and September that is similar to that produced in CM2.6 miniBLING. This 404 asymmetry is not seen in TOPAZ. Analysis of what drives this asymmetry shows that it is not 405 straightforward. All of the model runs show an asymmetry in the nutrient concentrations that 406 is in the opposite direction as the observations, with higher nutrients in February than in 407 September, as shown in Fig 7b. As we will show later in the manuscript, this is probably 408 associated with the models mixing to excessive depth during the wintertime. However, in 409 TOPAZ this does not produce an asymmetry in chlorophyll, while in BLING and miniBLING 410 it does. There are two possible reasons for this: 411 1) The equilibrium assumption, which means that biomass in both BLING and miniBLING is 412 not directly simulated. In TOPAZ, the growth of plankton during the spring is limited by the 413 biomass of phytoplankton, whereas in the fall TOPAZ continues to have higher heterotrophic 414 biomass (diagnosed from growth rates over previous months) that then grazes the plankton. In 415 BLING and miniBLING, by contrast, the biomass responds almost instantaneously to changes 416 in growth conditions. 417 2) Different handling of light limitation. In TOPAZ light limitation is calculated using the 418 instantaneous local light, whereas in BLING it is calculated using the mixed layer average 419 light. Preliminary results with a very coarse resolution model using BLING show that this 420 reduces the summer-winter asymmetry slightly, but is not sufficient to make the February 421 bloom smaller than the September bloom. 422 It is likely that all three of these factors-too deep winter mixed layers leading to too high 423 nutrients, too little light limitation and instantaneous response to changes in growth 424 conditions, are all responsible for the overly strong blooms in boreal winter in the Arabian 425 Sea. 426 To get a better sense of the mechanisms driving the blooms in the model, the biomass (mol P 427 kg -1 ) is compared with the light intensity in the mixed layer and the light-saturated 428 photosynthesis rate (carbon specific) (s -1 ) in Figs. 8a and b for January of year 195. The two 429 terms in Fig. 8 are the two terms in the model that affect growth rate. Because biomass in the 430 miniBLING model is a function of growth rate only, it is in fact sufficient to understand what 431 drives the growth in the model. The biomass production and mixed layer light intensity (Fig.  432 8a) are not meaningfully correlated parameters. On the other hand, the biomass and the light-433 saturated carbon specific growth rate ( Fig. 8b; indicating the degree of nutrient limitation) are 434 positively correlated. From this, it can be concluded that the blooms in this region are more 435 driven by nutrient rather than light, consistent with, for example, Gomes et al. (2008). This 436 suggests in turn that it is likely biases in nutrient supply that drive biases in productivity. 437 We can get more insight into nutrient biases by examining the individual tendency terms 438 associated with advection, vertical diffusion and subgridscale eddy fluxes and time rate of 439 change of nutrients. For simplicity, in this paper we combine the vertical diffusive flux 440 associated with small-scale mixing with that due to the mixed layer parameterization. In addition to having annual cycles that are different from observations, the models also differ 450 from data in terms of interannual variability. As shown in Fig. 10

Large-scale correlation 464
The relationship between SSHA and chlorophyll is quite different in the model as compared 465 to the satellite. Monthly variation in the cross-correlation of chlorophyll and SSHA for eight 466 consequent years in CM2.6 is shown in Fig. 11. As in the remote sensing, the model shows  Fig. 12c shows that E1 moved westward 494 during this period of time, and that the chlorophyll concentration was kept high within the 495 central part of the eddy. E1 appears to be created by the passage of a cyclone, similar to the 496 eddy observed by Wang and Zhao (2008)  However, at other latitudes, the largest blooms offshore are found along gradients in SSH 500 rather than being associated with maxima or minima. This suggests a different mechanism for 501 producing blooms in the model. Following Gaube et al. (2014). it appears that the eddy 502 stirring mechanism is dominant. Satellite data (i.e. see Fig. 3 for the month of May) provide 503 some hints for the existence of advective chlorophyll plumes being advected away from 504 coastal regions. As shown in Figs. 12a and b, high velocities in the marginal region between 505 adjacent cyclonic and anticyclonic eddies can cause such plumes in the GFDL models as well. 506 Why is the model only able to simulate the relationship between SSH and chlorophyll in the 507 Southern part of the domain? We hypothesize this is due to differences in stratification 508 between the two regions. The average water temperature (colors) and the macronutrient ( Fig. 13. In 512 the northern part of the region (see Fig 13a and b), the GFDL model provides a reasonably 513 good estimation of the mean temperature field near the surface, but subsurface temperatures 514 are not as consistent as there is far too little stratification. This is also associated with a very 515 weak nutricline in CM2.6. Variations in isopycnal depth will therefore not lead to big 516 differences in nutrient supply. Figs. 13c and d show the same results for the southern part of 517 the region. Unlike the northern part of the domain, the temperature gradient over these depths 518 is well estimated by CM2.6. While the nutricline is still too weak there is some gradient in 519 nutrients between 80 and 120m. 520 As seen in (Fig. 1d) both the ARGO and WOA09 wintertime mixed layer depth is 521 considerably deeper than that of summer, reaching a maximum of 65m. However, in the 522 northern regions the MLD seems to be too deep in winter, reaching values of 130-150 m. This 523 suggests that the overly deep mixed layer in the northern part of the region may explain both 524 the tendency towards an overly strong winter bloom and the failure of mesoscale eddies in 525 modulating chlorophyll blooms. If we look during the time period where we have eddies E1 526 and E2 (Nov- Dec. year 197, Fig 6c,d) we see shallower mixed layers associated with both 527 eddies. 528 Both the temperature and mixed layer biases in the northern part of the Arabian Sea may 529 result from having too much water from Persian Gulf in this region. This can be seen in the 530 yearly averaged subsurface salinity-density distribution over the region, shown in Figs 13e 531 and f for both WOA09 data and CM2.6 (model year 197), respectively. Fig. 13e shows two 532 separate tongues of salty water, one near the surface and one at the depth of ~300m. These 533 salty water signals are consistent with the seasonal cycle of Persian Gulf outflow as discussed 534 in Ezam et al. (2010). On the other hand, CM2.6 shows one subsurface salty water signal 535 from the northern part, which is deep and strong enough to result in weak stratification in the 536 north to a depth of 250m, as shown in Fig. 13f. These results suggest that a sharp thermocline 537 and nutricline is necessary for eddy activity to modulate the mixing of nutrients to the surface. 538 We test the idea that a sharper thermocline could modulate mixing of nutrients to the surface 539 by looking at the sources of nutrient in the southern part of Arabian Sea where eddy-bloom 540 relationships are seen. Accordingly, the region containing eddy E1 in Fig. 12 is analyzed to 541 determine the physical mechanisms by which nutrient is transported into the surface layer. 542 The bloom associated with eddies E1 and E2 do not fit with any of the mechanisms 550 highlighted in Gaube et al. (2014). We first consider the mechanism of trapping. Eddy E1 is 551 generated in the ocean interior, not as a result of coastal upwelling. As shown in Fig. 15, the 552 nutrient supply rate ranges between 5 and 8 mmol/m 2 /month in the eddy. The concentrations 553 in this eddy are only 0.01 µM (5 mmol/m 2 ) over the top 50 m. It cannot be the case that the 554 nutrients in the eddy can last for several months as a result of "trapping", there must be a 555 continuous supply. Moreover although eddy E2 shows a horizontal advection signal in 556 November (with a positive ring around the edge in Fig. 12a), the signal in December has the 557 opposite sign. Eddy intensification is also an unlikely mechanism for explaining the blooms, 558 as dSSH/dt is relatively small (particularly if we track the minimum SSH associated with E1 559 in Fig. 12c or E2 in Fig. 12d). Finally, Ekman pumping signatures in Gaube et al. (2014) have 560 the opposite sign as what is seen in E1 and E2. 561 Our results also contrast with those in Resplandy et al. (2011). The focus in Resplandy et al. 562 (2011) is on the productivity driven by horizontal and vertical advection in summer and 563 mostly vertical advection in winter. This contradicts our finding of a primary diffusive source 564 of nutrient in winter although it is consistent the finding of advective source of nutrients in 565 summer. We point out that in our model, the only two eddies that actually look like what we 566 see in the satellite observations involve enhanced mixing from below. This is a different result 567 from Levy et al. (2014) and Resplandy et al. (2011). Moreover it is not clear whether these 568 papers get the seasonal correlation with SSH or not. Resplandy et al. (2011) do not focus on 569 structures at the eddy scale as they are more concerned with the net impact of eddies. 570

571
To summarize, we hypothesize that 572 1. The reason that blooms are found in cyclones in the Arabian Sea during the NEM is that 573 the dominant source of nutrients to the surface, i.e. mixing (Barimalala et al., 2013;574 Kawamiya and Oschlies, 2003) is concentrated there. 575 2. Interannual variability in wintertime blooms in the Northwest Arabian Sea is controlled by 576 the combined presence of these eddies and strength of wintertime cooling. 577 3. Excessive mixing (resulting in too weak a thermocline) prevents mixing from being 578 modulated by eddies in the model except occasionally in the southern part of our region. In 579 the real world the modulation of mixing seen in Fig. 14

extends into the Northwest Arabian 580
Sea and the Gulf of Oman. 581

Conclusions 582
Our analysis of bloom variability in the northwestern Arabian Sea and Gulf of Oman has 583 illustrated a set of both similar and dissimilar descriptive features between satellites and 584 among a suite of models and explored the various mechanisms involved. Satellite analyses 585 verified the existence of two blooms, the stronger one associated with the Southwest monsoon 586 and the weaker one associated with the Northeast Monsoon as also shown by Madhupratap et 587 al. (1996), Kawamiya and Oschlies (2003), Murtugudde et al. (2007), and Al-Azri et al. 588 (2010). We demonstrate a pronounced anti-correlation between SSHA and chlorophyll 589 blooms in during certain times in northern winter but a much weaker relationship in other 590 months (typically northern summer) with the relationship disappearing as the blooms vanish 591 in the months of April and May (northern spring). While the depth of thermocline and 592 nutricline and also stratification are affecting the convection during the Northeast Monsoon 593 (Dickey et al., 1998;Kumar et al., 2001;Wiggert et al., 2002), we show that a thin 594 nutricline/thermocline and a strong stratification are also required to enable cold eddies to 595 bring nutrients to euophotic zone and develop phytoplankton blooms. During the wintertime 596 monsoon, while both cooling in the winter and eddies control the blooms, variability in bloom 597 location will arise from variability in the location of eddies, and so may not be predictable. In 598 contrast, during the Southwest Monsoon the dominant upwelling associated with the intense 599 environmental forcing supersedes the effect of eddies and the activity of the cold eddies is not 600 pronounced. 601 Understanding of this phenomenon has been sought using three different 3D ocean-602 atmosphere models, including a CORE-forced ocean with the TOPAZ biogeochemistry, a 603 coupled model with the TOPAZ biogeochemisty and CM2.6. Because the coarse models with 604 TOPAZ are not able to capture eddies and the interannual variability, CM2.6 (miniBLING), a 605 eddy-resolving high resolution model, was also considered for simulating the spatial and 606 temporal changes of the bloom in the region. This model simulates the two blooms seen in the 607 data and shows that the nutrients driving the northern summer bloom are supplied by 608 advection while those driving the wintertime bloom are supplied by vertical diffusion. 609 However, this model is unable to simulate the seasonal relationship observed in the satellite 610 products between blooms and sea surface height. Although there is some anti-correlation, it 611 tends to be associated with larger spatial scales and not really related to eddies. Instead, 612 eddies in the model usually wrap the chlorophyll around themselves, producing high 613 chlorophyll concentrations around their edges and not at their centers. Comparing the model 614 results to field measurements (WOA09) showed that the model does not account for the 615 strong thermocline and nutricline in the northern part of the region. In the wintertime, this 616 leads to excessive convective supply of nutrients and too strong of a bloom. However, for a 617 few cases eddies with a bloom at the center are tracked in the southern part of the domain. In 618 this region consistency is observed between the model results and the field data. Analysis of 619 the term balances in mixed layer show that eddies in this region modulate the diffusive supply 620 of nutrients. We suggest that what happens in the model in the Southern Arabian Sea actually 621 describes the Arabian Sea as a whole according to the observations and the field data. The 622 model misses the eddy signal in the north because it lacks a thin nutricline, motions of which 623 will lead to differences in nutrient supply. In the real world, eddies modulate the diffusive 624 supply of nutrients during the wintertime and there is more mixing in the eddy centers along 625 with the diffusive supply provided by the cooling in the wintertime. Accordingly, there is a 626 potential to improve the numerical models by better simulating the Persian Gulf Outflow to 627 produce a sharper thermocline, allowing more realistic nutrient supply. 628 It is worth noting that regional models, (such as Resplandy et al. (2011)) do have the potential 629 to better simulate the hydrography of the Northern Arabian Sea. Because such models are 630 very tightly constrained through "sponges" that restore hydrography at the boundaries, they 631 may not have the problems that global models do at representing the effects of overflows that 632 they do not properly simulate. However, such models cannot by themselves simulate the 633 effects of changing climate, which in turn changes the boundary conditions. For this reason, 634 global models must still be used for projection, making it important to identify the reasons 635 that they are not going to work.