SENSITIVITY OF UPPER OCEAN DYNAMICS IN HIGH-RESOLUTION TROPICAL INDIAN OCEAN MODEL TO DIFFERENT FLUX PARAMETERIZATION : CASE STUDY FOR THE BAY OF BENGAL ( BOB )

Simulation experiments using a high-resolution ocean general circulation model (OGCM) of the tropical Indian Ocean (TIO) were carried out to assess the model’s sensitivity to different flux parameterization. The flux formulation proposed by Kara et al. (2000) is used in the control run (CR). One more experiment differing in the bulk fluxes formulation for the computation of momentum, freshwater and heat is carried out. In the first experiment (CR), actual wind is used for the computation of the exchange coefficient in air-sea bulk flux formulation. In the second experiment (E1), model surface current is used in the wind stress formulation to compute the turbulent air-sea fluxes for TIO region. The formulation used in E1 is the same as it is used in CR, instead of actual wind, relative wind component is used in flux formulas. Both experiments are carried out for the period 2014-2016. The OGCM is forced using the daily fields of winds, radiation and freshwater fluxes obtained from ERA-Interim Reanalysis. In this study, we examine and quantify the performance of the above-mentioned experiments with respect to observations from ARGO, satellite-based sea surface temperature (SST) and sea surface salinity (SSS) for the year 2015. We observe that the upper ocean dynamics is significantly modulated by different flux algorithms. The errors in simulated SST is reduced by ~8% to 10% in E1 compared to CR, respectively. The temperature errors in the top 20m depth are reduced by 8% in E1. It is found that this flux formulation using relative winds is effective in accurately simulating the upper ocean dynamics in strong wind regimes of the Bay of Bengal.


INTRODUCTION
Momentum and heat flux exchange at the air-sea interface constitute key connector between the ocean and atmosphere and play a significant role in the circulation and the energy exchange between these two systems (Gill, 1982;Fairall et al., 1994Fairall et al., , 1996aFairall et al., , 1996bFairall et al., , 2003;;Kara et al., 2000;Siedler et al., 2013).In ocean general circulation models (OGCMs), surface wind stress and the heat fluxes are the major driving forces that modify the upper ocean dynamics and thermodynamics (Yuen et al., 1992;Chen et al., 1994;Shriver and Hurlburt, 1997;Swenson and Hansen, 1999;Kara et al., 2000).Therefore, both momentum and heat fluxes are essential for turbulent mixing in the upper ocean (Barnier, 1998).
In OGCM simulations, turbulent air-sea fluxes (TASFs) are generally implemented either one of the three ways: i) using a prescribed fields that are generated from a climatological data set (da Silva et al., 1994), atmospheric model (Kallberg, 1998), and remotely sensed and field (like satellite and buoy) observations (Cumming et al., 1997;Weller et al., 1998); ii) using bulk parametrizations that depend on the surface characteristic at the air-sea interface (e.g., heat fluxes (Hogan and Rosmond, 1991); and iii) using turbulence-based air-sea bulk flux formulations (Fairall et al., 1996a(Fairall et al., , 1996b(Fairall et al., , 2003)).The first method is exclusively used for OGCM simulations, where a gridded field of surface wind stress and heat fluxes are interpolated to the model grid at each time step.Other two * Corresponding author methods are generally used in OGCMs, where the essential atmospheric properties are supplied through a prescribed field or some other sources such as available meteorological surface state parameters, and a bulk transfer coefficient (Brodeau et al., 2017).
The benefit of using prescribed fields is the assurance of accurate forcing and computationally it is less expensive for the OGCM simulations.The major drawback of using prescribed fields prevent the cyclic feedback mechanisms between the coupled ocean-atmosphere system.This limitation allows the model and bulk formulas to determine the TASFs at each time integrating step of the OGCM simulation (Kara et al., 2000).McCreary et al. (1993) suggested that using model SST in the computation of surface fluxes of an OGCM simulation automatically provides a physically realistic tendency toward the correct SST simulation, that avoid model SST to drift.Therefore, the accurate notation of bulk formulation for the computation of TASFs is more important when performing large-scale OGCM simulations for the tropical oceans (Atlas, 1987;Atlas et al., 1987;Price et al., 1987) and the good knowledge of TASFs boosts the performance and the longevity of OGCM.
In general practice, the magnitude of ocean surface currents is at least one order smaller with respect to 10m wind and assumed to have a negligible effect on TASFs.A number of recent experimental studies have examined with this approximation and found that ocean surface currents in the wind stress formulation can minimize a positive bias in the estimate of input forcing to the ocean and drive the upper ocean dynamics (Luo et al., 2005;Duhaut and Straub, 2006;Dawe and Thompson 2006;Zhai and Greatbatch, 2007;Hughes and Wilson 2008;Scott and Xu, 2009;Zhai et al., 2012;Munday and Zhai, 2015) and the effect is so-called relative wind-stress effect).Pacanowski (1987) and Luo et al. (2005) have shown that the relative wind-stress effect in a revised wind-stress formulation reduces the strength of equatorial current (by ~30%) to a more realistic level, which has direct impact on sea surface temperature (SST) and the upper ocean dynamics, thereby improving the model simulations of the tropical ocean.Zhai and Greatbatch (2007) and Rath et al. (2013) reported a similar mechanism of reduction in the magnitude when they included ocean surface currents in the wind-stress formulation.The relative windstress effect brings relatively a small change in the wind-stress magnitude and heat flux by about ~2% when averaged over the North Pacific, it systematically reduces the magnitude of ocean surface currents, particularly in high oceanic wind zone (Zhai and Greatbatch, 2007;Eden and Dietze, 2009;Munday and Zhai, 2015;Xu et al., 2016).Duhaut and Straub (2006), Zhai and Greatbatch (2007), and Xu and Scott (2008) reported a significant reduction on the wind energy input by revised wind-stress formulation and also a damping effect on the eddy kinetic energy (EKE) by about ~10% by eddy/wind interaction in their northwestern Atlantic simulation.Eden and Dietze (2009) and Rath et al. (2013) quantify a similar amount of reduction in EKE in the sub-polar region.In the tropical region, the reduction is much higher (i.e., ~50%) or even more.
In comparison, there have been numerous modelling studies using bulk aerodynamic algorithm for the calculation of TASFs, accounting for the realistic simulation of the dynamics and thermodynamics of couple GCMs with better accuracy and ease of computation (Zeng et al., 1998;Kara et al., 2000;Shinoda et al., 2013;Chen et al., 2015;Jensen et al., 2016Jensen et al., , 2018)).Therefore, ocean surface currents and SST have significant impact on the parametrization of bulk formulas (Dawe and Thompson, 2006;Munday and Zhai, 2015;Brodeau et al., 2017).Since direct measurements of TASFs are difficult, unpredictable and highly sensitive to small variations in the physical condition of the medium, they cannot be used to build global or regional climatologies of TASFs (Brunke et al., 2003(Brunke et al., , 2011;;Brodeau et al., 2017).Rather, this information gives reasonable facts before the development and improvement of bulk air-sea flux parameterizations (Fairall et al., 1996b(Fairall et al., , 2003(Fairall et al., , 2011;;Edson et al., 2013).In a 0.20 0 North Pacific regional model, Dawe and Thompson (2006) found that including model surface currents in the bulk formulas actually reduces heat fluxes by about ~10% in the Kuroshio region, although the basin averaged heat flux reduction is quite small (i.e., 1%-2%).Interestingly, the reverse effect was observed in the tropical Pacific region where the heat fluxes increase because of surface warming that brings changes in the surface dynamics.Dawe and Thompson (2006) further suggested that changes in surface turbulent heat fluxes due to ocean surface currents in the bulk formulation may result from not only the direct effect of including ocean surface currents in the heat flux computation but also the indirect effect of including ocean surface currents in the wind stress computation, that modulate SST and hence surface heat fluxes via ocean dynamics.However, the relative importance of the direct and indirect effect are still uncommon and unpredictable.
Differences in the state-of-the-art bulk algorithms have been repeatedly shown to lead to a relatively large spread in the computed TASFs (Blanc,, 1985(Blanc,, , 1986(Blanc,, , 1987;;Zeng et al., 1998;Eymard et al., 1999;Brunke et al., 2002Brunke et al., , 2003)).In spite of this, we ignore large uncertainties in meteorological surface state variables from different sources (Josey et al., 1999(Josey et al., , 2014;;Berry and Kent, 2005), which has now become a common practice to neglect the accuracy of various constants and parameter approximations used in bulk formulas.On the other hand, the uncertainties associated with these conditions and assumptions are not properly diagnosed to constitute actionable information for concerned users.
So far, the focus of previous studies on this topic has been primarily to quantify the performance of coupled GCM, and to some extent, on Tropical Ocean using air-sea bulk flux algorithm.There have been few studies reporting the sensitivity of the boundary layer dynamics of the coupled system to TASFs flux algorithm.Furthermore, the previous studies often rely on the sensitivity of North Pacific, North Atlantic and the Southern Ocean model simulations [e.g., Dawe and Thompson, 2006;Zhai and Greatbatch, 2007;Eden and Dietze, 2009;Munday and Zhai, 2015] to bulk flux algorithm.Here we investigate for the first time the sensitivity of upper ocean dynamics to different bulk flux algorithm in the Bay of Bengal.With this concern in mind, we examine and quantify the performance of the dynamics of upper ocean simulation to different bulk flux algorithm with relative wind effect using high-resolution (i.e., 0.10 0 ) stand-alone Tropical Indian Ocean model.provides a brief model description and experiment design.Section 4 discusses the results and important findings.Section 5 concludes with a brief summary of our results.

Satellite Data
A daily product of Kd_490nm prepared from the merged (MERIS, MODIS Aqua, SeaWiFS and VIIRS) fully normalized remote sensing reflectance products at 443, 490, 555 and 670 nm using the semi-analytical method of Lee and Arnone (2005).Kd_490nm is available from GlobColour through the HERMES web interface.In this study, we have used Level-3 (L3) products of regular spatial resolution of 4km.The merged data product is available at 8-day intervals (ftp://hermes.acri.fr/?class=archive).The merged data had gaps which were filled by making a daily climatological dataset for the global domain.
A daily product of the Group for High-Resolution Sea Surface Temperature (GHRSST) data incorporate SST observations from various available sources like microwave and infra-red sensors.The Level-4 (L4) product is generated using various objective analysis techniques to produce gap-free SST maps over the global oceans.In this study, we have used L4 GHRSST product with a regular spatial resolution of 1 km (Donlon et al., 2014).The data was downloaded from the Physical Oceanography Distributed Active Archive Center (PODAAC) (https://podaac.jpl.nasa.gov/GHRSST).
Monthly climatology of river runoff for Ganga-Brahmaputra derived from altimeter observations (Papa et al., 2010) has been used in the model.

In-situ Data
The Coriolis Ocean Dataset for Re-Analysis (CORA) product is a global dataset of in-situ temperature and salinity measurements.The CORA observations come from different sources collected by the Coriolis Data Centre in collaboration with the In-Situ Thematic Centre of the Copernicus Marine Service (CMEMS).This dataset is revised every year by the R&D Coriolis team, and data are extracted in the NetCDF Argo format at a given date.The dataset contains data from different types of instruments: mainly Argo floats, Expendable Bathythermographs (XBT), Conductivity-Temperature-Depth (CTD) and eXpendable-Conductivity-Temperature-Depth (XCTD) profilers, and moorings (Cabanes et al., 2013).The data set is available on a monthly basis.In order to re-initialize the OGCM a monthly climatology of temperature and salinity was generated using the CORA dataset from 2006 -2012.
In order to validate model subsurface parameters, in-situ measurements of temperature and salinity profiles were used.These datasets were obtained from the National Oceanographic Data Center (NODC) under the Global Temperature and Salinity Profile Program (GTSPP) that provides temperature and salinity profile data.These data are compiled from several in-situ platforms including XBTs, CTDs, drifters, ARGO floats (Sun et al., 2010).The data are available from Ocean Archive System (OAS) (ftp://ftp.nodc.noaa.gov/pub/gtspp/indian/).

Model Description
The model used in this study is based on the Modular Ocean Model (MOM, version 3.0) (Pacanowski and Griffies, 1999).This is a primitive equation, a hydrostatic model with a free surface.For this study, the model is configured for the Indian Ocean (27 0 E -142 0 E, 35 0 S -26 0 N) with a variable horizontal resolution that it is uniform 0.1 0 x 0.1 0 in the region 55 0 E -100 0 E and 10 0 S -23 0 N and reaches a minimum up to 1.0 0 in the rest of the domain.It has 50 z levels in the vertical grid (vertical resolution ranging from 2.5 m at the surface to 5320 m at the bottom, with top 14 levels in the upper 36 m, shown in Table -3).The model uses the Arakawa B-grid and leapfrog scheme for numerical simulations which makes it easy to solve the momentum equations and the well mixed, weakly stratified layer.The vertical mixing based on nonlocal K-profile parameterization (Large et al., 1994) known as KPP mixing scheme has been used.Bathymetry data are taken from ETOPO-2.(Data Announcement 88MGG-02, Digital relief of the Surface of the Earth, NOAA, National Geophysical Data Center, Boulder, Colorado, 1988).
Using climatological temperature and salinity from Levitus (1982) and CORA dataset (as discussed in section 2.2.2), the model was re-initialized (Balmaseda et al., 2007), and spun up from rest and forced with ERA-Interim daily analyzed wind fields, air temperature, specific humidity, net solar and thermal radiation to obtain the initial condition for 2014.The daily climatology of diffused attenuation coefficient (Kd_490; as discussed in section 2.2.1) was used to determine the penetrative depth of blue-green component of solar radiation in the model.The penetrative depth of shortwave radiation is exponentially computed using shortwave parameterization scheme.The detail description of shortwave parameterization is mentioned in section 3.2 (Paulson and Simpson, 1977).Kd_490 varies with space and time, depending upon the nutrient-rich active zone and biogenic sediment deposits in different regions.The other attenuation parameter like R is then determined from Jerlov (1968) look-up table, including the use of 0. 35 m for ζ1 in prescribing the absorption of the red light.Monthly climatology of altimeter derived river discharge data from Ganga, Brahmaputra (Papa et al., 2010) was used in the model.Runoff discharge from other prominent rivers that flow into BoB (Irrawaddy, Damodar and Godavari) has been taken from Vorosmarty et al. (1998).The river discharge is distributed horizontally around the exact river mouth locations in the first model layer.The model was then further integrated from 2014 until December 2016.This base experiment is referred to as a control run (CR).

Numerical Experiments
To investigate the sensitivity of upper ocean dynamics using TIO model to different flux parameterization, we performed a sensitive experiments E1 with reference to control run (CR).The details description of these experiment set up is mentioned in Table 2.In the first experiment (CR), the daily climatology of Kd_490nm (as discussed in section 2.2.1) and 10m atmospheric wind are considered for the simulation.In CR, the actual wind component is used for the computation of the exchange coefficient in the bulk flux formulation.In the second experiment (E1), the exchange coefficients are computed based on the relative wind effect on air-sea bulk flux formulation proposed by Dawe and Thompson (2006) (Eq. 1).The bulk flux formulation used in CR and E1 is well wellknown Kara et al. (2000) algorithm.Both CR and E1 are initialized with the same the initial condition and ERA-Interim forcing and simulated for 3 years from 2014 to 2016.And the analysis is done for the year 2015.),   is the specific humidity of 2m air,   is the saturated humidity at SST, sensible heat   , and latent heat   , respectively,  10 −   is the relative wind (wind velocity at reference (i.e., 10m) height with respect to the first depth level of model surface velocity),   −   and   −   are the differences in potential temperature and humidity between the atmosphere, and the ocean surface (first depth level of OGCM), respectively, and heat flux () is positive for the warmer ocean.Equation ( 1) states that the surface momentum and turbulent air-sea heat fluxes (TASFs) depends on the relative wind component.Different bulk formulae's, proposed by climate modeler's (Large and Pond, 1982;Zhang and McFarlane, 1995;Zeng et al., 1998;Fairall et al., 1994Fairall et al., , 1996aFairall et al., , 1996bFairall et al., , 2003)), differ mainly in the methods used to compute the transfer coefficients (like,   ,  ℎ and   ) that vary with the magnitude and physical condition of the atmospheric boundary layer (Dawe and Thompson, 2006).

Shortwave Parameterization Scheme
The subsurface heating parameterization used here was proposed by Paulson and Simpson (1977).It is represented as: where Iz is the downward irradiance Io is the incident irradiance, ζ1 and ζ2 are attenuation lengths for blue-green and red light respectively and R is an empirical constant.Diffused attenuation coefficient usually kept vertically constant but spatially vary with water types Jerlov (1976).Depending on the wave-length and biogenic components of the water, Jerlov classifies the oceanic water types into five types.Generally, the water types are classified as Jerlov I, IA, IB, II and III.Type 1(Jerlov I) represents clear water in open ocean type V (Jerlov III) is high turbid water.
The value of R is taken to be 0.58 and that of ζ1 and ζ2 is 0.35m and 23m for type I (Jerlov I) waters respectively.The divergence of downward irradiance may be depicted as: this is added to the source term as where   is the ocean surface net shortwave radiation, this is then finally added to the tracer equation where U,   ,   and   are the ocean currents, source term, horizontal and vertical eddy diffusivity, respectively.The lefthand side of ( 5) is the local temperature gradient and on the right-hand side of (5), the first, second, third and fourth terms are characterized as the advection, turbulent, diffusion and the source terms.

Momentum Flux Differences between the CR and E1
Including ocean model surface current in the bulk formulas leads to a slight but wide-spread weakening of the mean surface wind stress (Figure 1).This weakening in surface wind stress is maximum in those regions where ocean surface currents are relatively strong.For example, the strength of the mean wind stress averaged over the South-West BoB near the Southern tip of Srilanka region (5 0 -8 0 N, 80 0 E-86 0 E) is about ~8% to 12% weaker in E1 than in CR, in agreement with the 5%-10% difference reported by Dawe and Thompson (2006).Similar percentage decreases are also found in the East India Coastal region is about ~6% to 8%.When averaged over the BoB (5 0 -25 0 N, 78 0 E-100 0 E), the mean wind stress in E1 is about almost 4% weaker than those in CR.Interestingly, figure 2(b) makes it clear that the differences in surface heat flux between CR and E1 are not a direct effect of including ocean surface currents in the turbulent heat flux calculations, but mainly an indirect effect of including ocean surface currents in the wind stress calculation via surface current differences between the two experiments Dawe and Thompson (2006).Differences in surface heat flux between CR and E1 are closely linked to their SST differences (Figure 3).The pronounced heat flux differences in the entire BoB region and the western boundary current regions are in opposite phase to the SST differences, that shows changes in SST is one of the primary contricontributor which modulate surface heat fluxes in the basin.This further infers the SST differences between CR and E1, especially in the coastal region, due to wind stress differences (Figure 1) induced by relative wind in the flux formulation.
Figure 3. Time-mean sea surface temperature (SST, 0C) of (a) CR and the time-mean differences between (b) CR and E1 for the year 2015 The results above demonstrate that the feedback of model surface currents in the wind stress formulation dominates the differences in surface heat fluxes and the SST over the direct effect incorporating model surface currents in the heat flux calculation.As the differences are very small in magnitude, it is clear that the sudden response of surface heat flux is almost entirely governed by the relative wind component in the bulk formulation.

Comparison with In-Situ Observations
The vertical profiles of the domain-averaged (  -  N,   E-  E), temperature root mean square error (RMSE) between the experiments (i.e., CR and E1) and ARGO profiles are illustrated in Figure 5 (a).The two runs (E1 and CR) have similar, small cold biases near the sea surface.In the subsurface, in between 20-100m, error in temperature profile is increased by 0.5 0 C, this may be due to decrease in wind stress magnitude that hinder the upper ocean vertical mixing and cause anomalous warming in the subsurface.At 100m depth, the temperature error in E1 (~ 1.7 0 C) is slightly less compared to CR (~ 2.0 0 C).One can also see that at this depth the natural variability of temperature is more than the errors in the E1.The greatest differences occur at the depths ranging from 20m to 100m.Correlation has significantly improved in the upper ocean in the E1 up to a depth of 10m.A different approach such as the one proposed by Fairall et al. 1994Fairall et al. , 1996aFairall et al. , 1996bFairall et al. , 2003. to compute the exchange coefficients is required to be tested in the OGCM.

Sl
The paper is organized as follows.Section 2 provides a brief description of input forcing and validation data sets.Section 3

Figure 1 .
Figure 1.Magnitude of time-mean surface wind stress (dyn/cm2) of (a) CR and differences in the magnitude of the time-mean between (b) CR and E1 for the year 2015 4.2 Heat Flux Differences between the CR and E1 Figure 2 shows the time-mean (or annual average) net surface heat flux averaged over the year 2015 in CR and the differences between the E1.Although the overall patterns of the mean surface heat fluxes in E1 experiments are very similar, including ocean surface currents in air-sea flux calculations lead to a significant reduction in surface heat loss in the eastern and western basin of BoB region.Whereas in the central BoB it is increased by ~10 -14%.When averaged over the BoB (5 0 -25 0 N, 78 0 E-100 0 E), the mean heat flux in E1 is about almost ~6% higher than those in CR.

Figure 2 .
Figure 2. Time-mean Heat flux (Watt/m2) of (a) CR and the time-mean differences between (b) CR and E1 for the year 2015

Figure 4
Figure4shows the spatial map of RMSE between modelsimulated SST and GHRSST for the year 2015.It is clearly observed that the relative wind effect has actually minimized the SST error by around ~ to 12% in the south-west BoB region.In the central BoB, the errors have reduced by about ~6%-8%.When averaged over the entire BoB (5 0 -25 0 N, 78 0 E-100 0 E), the mean SST error in E1 is about almost 2% lesser than those in CR.However, errors are quite large in both CR and E1 near the coastal regions for reasons that are not yet explored in this study.Both CR and E1 have higher correlation values in the head as well as in central bay (i.e., > 0.90, figure5).

Figure 4 .
Figure 4. RMSE map of SST for CR and E1 with GHRSST for the year 2015 Figure 5 (b) represents the vertical profiles of root mean square error (RMSE) and correlation of simulated salinity from two different runs.The surface salinity error and correlation value in E1 is increased by 4% at this location as compared to CR. Below the surface, there is not much difference in salinity from E1 and CR at this location.

Figure 5 .
Figure 5. Vertical RMSE profile of simulated model a) Temperature and b) Salinity for the year 2015 for Bay of Bengal (BoB) 5. CONCLUSION In this study, we have investigated the impact of including model surface currents in the bulk formulas on air-sea fluxes and the ocean general circulation model using a high resolution standalone TIO model.By comparing the model simulation that include and exclude model surface currents in air-sea flux formulation, we find changes in wind stress, heat flux and sea surface temperature are small when integrated over basin, but in some confined regions it alter by ~10%-15% of the mean or even more.Wind stress is reduced over high wind zone region over the Bay of Bengal region, resulting in ~10% reduction in the magnitude of coastal current.Heat flux to the atmosphere in the central Bay of Bengal reduced by ~10%-15%.We find in most of the Bay of Bengal region, SST increase when the relative wind stress effect is included, although the warming in our model is increased by 6%.The model's increase in both Bay of Bengal SST and outgoing heat flux is interesting to note that increased heat flux to the atmosphere would tend to lower SST.It appears that including relative wind component in the heat flux formulas initially reduces heat flux out of the ocean, causing SST to warmer by ~9 %.The warmer SST further opposes the relative wind stress effect, hence increasing the outgoing heat flux.The difference in wind stress on air-sea fluxes formulation in Bay of Bengal can significantly modify the upper ocean circulation pattern that affects the vertical mixing and the mixed layer depth (MLD).At the same time, though, that heating recharges the submerged warmer water that, when exposed by storms, may add additional thermal energy to storm growth.It is seen that although at several places the errors in SST have reduced by including the effect of ocean currents in the bulk formulae, there are still large errors in the model simulations in other regions.A different approach such as the one proposed byFairall et al. 1994Fairall et al. , 1996aFairall et al. , 1996bFairall et al. , 2003. to compute the exchange coefficients is required to be tested in the OGCM.

Table 3 :
Vertical discretization of the water column in our model.transfer or exchange coefficients for wind stress ,  10 is the 10m atmospheric wind,   is the model surface current,   is 2m air temperature,   is model surface temperature (or model SST