A Comparative Study of Vertical Mixing Schemes in Modeling the Bay of Bengal Dynamics

The choice of vertical mixing scheme in ocean models plays an important role in modeling the surface and subsurface circulation and the vertical structure. This work performs a comparative study between K‐profile parametrization (KPP) and k‐ϵ mixing schemes for a regional domain in the Bay of Bengal (BoB) using the Modular Ocean Model version 5 (MOM5). It is observed that sea surface temperature (SST) and the mixed layer depth (MLD) show significant improvement with the k‐ϵ mixing scheme. Energetic analysis shows that changes in the viscous dissipation and turbulent buoyancy flux are the primary reason for improvement with k‐ϵ. The overestimation of viscous dissipation in the KPP scheme is corrected by k‐ϵ, resulting in a deeper mixed layer closer to observations. The tendency of buoyancy flux to retain stability in the water column also results in a better representation of SST in k‐ϵ. Overall, we conclude that the k‐ϵ mixing scheme works better for the BoB region.

of fresh water during summer due to greater heat absorption within a shallow mixed layer (e.g., Howden & Murtugudde, 2001). The BOBMEX experiment looked into the heat advection in BoB and its influence on MLD (Bhat, 2003;Bhat et al., 2001). Measurements taken during BoBBLE (Vijith et al., 2020;Vinayachandran et al., 2018) highlighted the SST variability and formation of the barrier layer. Salinity effects in the Indian Ocean were modeled by Murtugudde and Busalacchi (1998) and Han et al. (2001). Girishkumar et al. (2013) used RAMA buoy measurements to document temperature inversions and their influence on the mixed layer heat budget. Kumar et al. (2016) examined the long-term variability of BoB SSTs to reveal that they are influenced predominantly by seasonal forcing.
Most of the modeling studies so far for the BoB have used the K-profile parametrization (KPP) vertical mixing scheme (Behara et al., 2019;Behara & Vinayachandran, 2016;Chatterjee et al., 2013;Kurian & Vinayachandran, 2007;Mukherjee et al., 2018;Srivastava et al., 2018) to represent turbulence and the associated mixing. It has been observed from these studies that the KPP scheme shows a relatively large bias in the ocean properties, especially in the southern and central regions of the BoB. It is primarily due to the parametric nature of KPP, which does not model the regional-scale processes. Therefore, it is imperative to test other vertical mixing schemes that could improve the representation of ocean processes and associated dynamics, reflecting improvements in the ocean properties such as SST and MLD. Among the other available mixing schemes, k-ϵ turbulence closure has been known to work well for modeling ocean dynamics. The k-ϵ scheme has been adopted for the Liverpool Bay, the Mediterranean Sea, and some Pacific Ocean studies (Fernández et al., 2006;Simpson et al., 2002;Walsh et al., 2017) and it has been documented that it performs well for these domains. However, none of these studies contrasted k-ϵ with other mixing models for the BoB domain. Fernández et al. (2006) compared KPP and k-ϵ for the Mediterranean Sea and found that KPP produces slightly deeper MLDs in summer than observations, which is rectified by using k-ϵ. The study mentioned that the differences in MLDs are due to changes in values of viscosity and diffusivity at the base of the mixed layer, indicating suppressed entrainment with the KPP mixing scheme.
To the best of our knowledge, the k-ϵ mixing scheme for modeling BoB dynamics has not been reported. Given that BoB has a multitude of processes that generate turbulence and mixing, it would be prudent to test a mixing scheme that focuses on modeling the turbulent processes. Therefore, the objective of the present study is to implement k-ϵ into MOM5 and perform an intercomparison of the results with the most trusted mixing scheme in ocean modeling (KPP). Mechanisms leading to the differences observed between the two mixing schemes are diagnosed. The rest of the paper is structured as follows: Section 2 provides a model description, data set, experiments performed, and model setup information. Section 3 is focused on model output comparison for surface properties like currents, surface temperature, and mixed layer depth. Subsurface analysis with vertical temperature profile and turbulent flux quantification is discussed in Section 4. A summary and conclusion are provided in Section 5.

Data and Methodology
Modular Ocean Model (MOM) version 5 (MOM5), developed by NOAA's Geophysical Fluid Dynamics Laboratory (GFDL), is used for this study (for details on MOM5, refer to Griffies (2012)).

Model Domain Specification
The model domain considered for the present study is a part of the Indian Ocean from 10°S to 25°N and 75° to 95°E. While the model run considers this entire domain, the analysis and discussions in the forthcoming sections are only shown for 0° to 25°N and 75° to 95°E. The Arakawa B-grid is used for computing tracers and velocity. A uniform horizontal grid resolution of 0.25° is employed in zonal and meridional directions. The nonuniform vertical grid with a 5 m uniform grid till 60 m depth is utilized. The Indian Ocean ETOPO2V2 (Sindhu et al., 2007) is used for ocean topography.
For open ocean interactions, the radiation boundary condition is used with the Orlanski (1976) Srivastava et al., 2018;Thompson et al., 2006) relaxed tracers at the boundaries to prevent systematic drift. This study also has relaxed temperature and salinity with a 30-day timescale. This flux correction method considers both models simulated and prescribed data as explained by Griffies (2003).
The model setup uses the β-plane approximation to represent the Coriolis effect. As base case, the model uses the KPP scheme for simulating vertical mixing with a background vertical viscosity of 10 −4 m 2 s −1 ; background diffusivity is set to 10 −5 m 2 s −1 ; constant for pure convection is 1.8 (Large et al., 1994); critical Richardson number is chosen as 0.25 (Kundu & Cohen, 1990). For comparative purposes, another vertical mixing scheme from the General Ocean Turbulence Model (GOTM)  is used, namely, the k-ϵ turbulence model. Technical documents with source code of GOTM are available in public domain: https://gotm.net and the coupling between GOTM and MOM has been described in Villarreal et al. (2005). Following the work by Rodi (1987), the constants used in k-ϵ model are chosen as σ k = 1, and σ ϵ = 1.3; empirical constants, c ϵ1 = 1.44, c ϵ2 = 1.92, c ϵ3 = 0.

Data Set Used in Study
The model uses annual averaged temperature and salinity profiles taken from World Ocean Atlas (WOA18) Zweng et al., 2019) to initialize the model simulation. For shortwave penetration into the upper ocean, a chlorophyll-based scheme (Morel & Antoine, 1994) has been used. Wind stresses are prescribed from Quick Scatterometer (QuikSCAT), short wave, long wave, sensible, and latent heat fluxes from WHOI OAFlux, precipitation from Tropical Rainfall Measuring Mission (TRMM), and chlorophyll from MODIS data set. All input parameters are available at http://apdrc.soest.hawaii.edu/data/data.php. Simple Ocean Data Assimilation (SODA) reanalysis data (Carton et al., 2018) and North Indian Ocean Atlas (NIOA) monthly Climatology (Chatterjee et al., 2012) are used to validate model output. The differences between the model and validation data are analyzed to depict model simulation accuracy and ocean physics. The model simulations are spun up by converting the model forcing data from 2003 to 2009 to climatological data for performing the control run for 10 years. To ensure the stability of total kinetic energy from the control run, experimental runs with interannual data for 2003-2008 are performed, and the output is analyzed for these 6 years.

Ocean Properties
We denote the seasons as DJF (December-January-February), MAM (March-April-May), JJAS (June-July-August-September), and ON (October-November) based on the direction of wind flow. The response in resultant ocean currents was analyzed (figure not shown). Depth averaged surface currents agree with observed basin-scale seasonal features like EICC and the strong equatorial currents such as the spring and fall Wyrtki Jets. There are no significant differences between the k-ϵ mixing scheme and KPP mixing scheme in ocean currents (Supporting Information S1). Also, sea surface salinity (SSS) (Supporting Information S1) has no significant difference between the two mixing schemes compared with the NIOA data.

Sea Surface Temperature
Sea surface temperature is an important property determined by the surface and subsurface ocean processes, especially in the BoB, where the high mean SSTs drive a strong response in the atmosphere. The BoB is known to exhibit a variety of mesoscale processes that also vary meridionally, which control the SST in this region. The model and observations show cooler SST in DJF and warmer SST in MAM. This is due to the variation in the amount of radiation received and the latent heat losses in different seasons (Murtugudde & Busalacchi, 1999). The north-south temperature gradient in DJF and warmer temperature near the equator in all seasons are also well represented in our model simulations. Figure 1a shows the comparison of observational SST from NIOA data and KPP model output. We observe that the KPP model tends to predict cooler SST (positive bias) except in the North Bay for all seasons. KPP shows warmer SST (negative bias) near the northern Bay than observational data. This negative bias could partly be due to the absence of river runoff forcing in the model setup. However, it is prudent to note that salinity restoration has been implemented in our study that acts as a proxy for river runoff. For similar forcing conditions, Figure 1b shows the SST difference between KPP and k-ϵ mixing schemes. In all the seasons, bias is reduced in the k-ϵ mixing scheme compared to KPP. The improvement is especially evident in central and southern BoB, where k-ϵ shows warmer SST (negative bias) that is closer to observations.

Mixed Layer Depth
Next, we consider MLD variability in the two vertical mixing schemes. The temperature criteria (Δ = 0.2 • C) (Courtois et al., 2017;de Boyer Montégut et al., 2004) has been used for MLD computation. Both the mixing schemes well replicate the basic pattern of MLD, but the differences are subtle and important. During DJF, MLD is higher in the northern part of the Bay, whereas, during JJAS, MLD is higher in the southern part of BoB. During MAM, the model-generated MLD is relatively on the shallow side. These MLD seasonal patterns show a positive correlation with seasonal winds. The more the wind stresses, deeper the MLD. Figure 1c shows the difference between SODA reanalysis and KPP model simulation. It clearly shows that the KPP mixing scheme renders shallower (positive bias) MLD than reanalysis data. This bias is larger in JJAS and the near-equatorial region when winds are stronger. In some regions, the bias is as high as 30 m, which is significantly large. Figure 1d shows comparison of MLD from KPP and k-ϵ mixing schemes. There is a clear reduction in bias, by as much as 10 m, with implementation of the k-ϵ scheme. We see that the k-ϵ scheme has helped deepen the MLD by up to ∼10 m over most of the domain for all the seasons, thereby leading to a closer match with the reanalysis data. The impact on the MLD heat and momentum budgets and the circulation are significant (not shown).

Characterization of Vertical Structure and Turbulent Fluxes
In order to understand and quantify the impact of the k-ϵ mixing scheme in improving the SST and MLD, we look at the vertical structure and turbulent fluxes. Figure 2 shows 5° to 10°N averaged longitude-depth seasonal plots for JJAS and ON seasons, where we see the most significant improvement with the k-ϵ mixing scheme. Figure 2 shows that near 85°E, the observed thermocline is deeper, while the KPP model leads to a shallow thermocline. The k-ϵ scheme deepens the thermocline, resulting in a better vertical structure closer to observations. In general, we note that the model output from k-ϵ better represents the vertical structure of temperature compared to KPP. This is one possible reason for improving SST and MLD while using the k-ϵ mixing scheme. The N 2 vertical structure figure is attached (Supporting Information S1). We see that at the base of the mixed layer a negative bias is seen, which is a representation of more stable water column produced by the k-ϵ mixing scheme. The stable water column shows a better representation of the vertical temperature profile, thereby resulting in improvements in both SST and MLD.

Quantification of Turbulent Fluxes
Turbulent fluxes play an important role in governing the mixing in the ocean and hence the controlling the evolution of SST and MLD. The turbulent kinetic energy, also referred to as the eddy kinetic energy (EKE), formulated below in Equation 1, provides a mechanistic approach to understanding the role of turbulent fluxes on ocean energetics and stability (Pope, 2000). 10.1029/2022EA002327 6 of 10 of pressure, and density, ρ 0 is the reference density (taken at the surface), ν is the kinematic viscosity, and e ij is the strain tensor described as follows: In this study, we look into the seasonal variations of K, P, B, and ϵ to understand the mechanisms of turbulence production and dissipation. The following methodology is employed to compute the seasonal mean of EKE budget terms. First, the velocity and density fields are time-averaged over a particular season (highlight JJAS) to obtain a single mean-field. Any mean-field, say , takes the form ( ) = 1 ∫ , where T is the total duration, and dt is the time step at which each field is available. The fluctuation field is obtained as First, we consider the eddy kinetic energy, K, which denotes the turbulent energy available for mixing. Figure 3a shows the difference in K between KPP and k-ϵ mixing schemes. It is noticed that there is no basin-wide change in K since the difference is mostly zero. We see some differences very close to the equator, which could be attributed to strong equatorial dynamics. This indicates that the total energy available for mixing is the same in both the mixing schemes. However, the contributions of the individual terms, namely, P, B, and ϵ, are important in understanding the improvements brought about by the k-ϵ scheme.
The turbulent buoyancy flux, B, represents the potential energy in the system and denotes the water column's stability. Figure 3b shows the difference in buoyancy flux between the two mixing schemes. The basin-wide difference is primarily negative. This means that the k-ϵ mixing scheme retains stability in the water column compared to KPP. The relatively stable water column better represents the vertical temperature profile, thereby resulting in improvements in both SST and MLD. The improved stability in k-ϵ is a direct consequence of the increase in the potential energy that helps maintain the thermocline.
The dissipation, ϵ, is a sink term in the EKE budget equation. It is a process by which eddy kinetic energy is lost to viscosity. Figure 3c shows the difference in ϵ. It is seen that the KPP model shows higher viscous dissipation compared to k-ϵ. This fundamentally means that energy is getting dissipated at a higher rate within the water column, thereby restricting the energy available for mixing to a shallower depth. This is the reason for the underestimation of MLD in the KPP model. On the other hand, k-ϵ allows the energy to penetrate deeper depths, resulting in a better representation of the vertical structure and the MLD closer to observations. The turbulent production flux, P, does not show any difference between KPP and k-ϵ, so we do not discuss it further. The negligible difference in P means that the shear production mechanism in both KPP and k-ϵ is similar, which is understandable since no significant change in ocean currents is seen from one scheme to the other.
The quantification of turbulent fluxes provides another avenue for understanding the superior performance of k-ϵ over KPP for estimating SST and MLD in the BoB. Along with the improvement in the vertical structure (see Section 4.1), the better representation of turbulent buoyancy flux, B and viscous dissipation, ϵ are the main reasons for k-ϵ vertical mixing scheme showing better results.

Summary and Conclusion
A comparative study of two different vertical mixing schemes, namely KPP and k-ϵ, was carried out for the Bay of Bengal (BoB) domain using MOM5. It was noted that both the mixing schemes do an excellent job of capturing the basin-wide patterns of surface and subsurface features like the currents, SST, SSS, and MLD. Quantification of SST and MLD revealed that the k-ϵ mixing scheme shows a closer match to observations than KPP for all seasons. In some regions of BoB, a significant change in SST and MLD was brought about by implementing k-ϵ, which showed the superiority of this scheme compared to KPP. The reasons for better performance of k-ϵ were probed. It was documented that the improvements in the vertical temperature structure and better representation Daily (4 km) Chlorophyll concentration, oci algorithm Table 1 Details of Forced Variable of turbulent buoyancy flux and viscous dissipation are the prime reasons for the k-ϵ mixing scheme showing a closer match to observations. This study concludes that the k-ϵ mixing scheme is better for modeling BoB dynamics. Further contrasts are needed at subseasonal and interannual timescales to establish the superiority of k-ϵ under varying background states, especially since the BoB SSTs in monsoon variability at these timescales is of great interest.

Data Availability Statement
The details of the forcing data set used in this study is given in Table 1. This includes the type of variable, version, frequency and data name. In order to access this data set, which is freely available use the below link: http:// apdrc.soest.hawaii.edu/data/data.php. For example, to obtain 'Initial temperature' data go to. http://apdrc.soest. hawaii.edu/data/data.php select data set from Table 1 that is, "World Ocean Atlas (Levitus; WOA)" and click on "LAS" − > then select version "NODC World Ocean Atlas 2018" − > select resolution "0.25°" − > then select data frequency "Annual" − > select year range "2005-2017" − > select data name "Objectively analyzed mean sea water temperature" and select "Next." OR click on variable name in Table 1; select the same variable from list and click "Next." Finally select file format, specify latitude and longitude, select exact year and submit query to download the data. The NIOA data used for validation can be accessed from the link below: https://publication-data.nio.org/s/q7g6t84j4YTkiGz. All the information about the Modular Ocean Model V5 employed for this study is available at below link: https://github.com/mom-ocean/MOM5.