Closing the sea surface mixed layer temperature budget from in situ observations alone: Operation Advection during BoBBLE

Sea surface temperature (SST) is a fundamental driver of tropical weather systems such as monsoon rainfall and tropical cyclones. However, understanding of the factors that control SST variability is lacking, especially during the monsoons when in situ observations are sparse. Here we use a ground-breaking observational approach to determine the controls on the SST variability in the southern Bay of Bengal. We achieve this through the first full closure of the ocean mixed layer energy budget derived entirely from in situ observations during the Bay of Bengal Boundary Layer Experiment (BoBBLE). Locally measured horizontal advection and entrainment contribute more significantly than expected to SST evolution and thus oceanic variability during the observation period. These processes are poorly resolved by state-of-the-art climate models, which may contribute to poor representation of monsoon rainfall variability. The novel techniques presented here provide a blueprint for future observational experiments to quantify the mixed layer heat budget on longer time scales and to evaluate these processes in models.

Here T, ρ 0 , h and c p are the temperature, mean density, the mixed layer depth (MLD) and specific heat capacity of sea water, respectively. u and v denote the horizontal components of velocity. The suffix a denotes a vertically-averaged (over the ML) quantity and the subscript −h denotes the quantity at the base of the ML (see Methods). q 0 is the net heat flux at the oceanic surface and q pen represents the penetrative loss of shortwave radiation. κ H and κ Z are the horizontal and vertical eddy diffusivities, respectively (see Methods).
To estimate this heat budget, in situ time-series observations at a single (ship or mooring) location are typically combined with satellite or model data sets, in order to obtain horizontal gradients [9][10][11][12] , but this approach introduces large uncertainties and cannot resolve important small-scale or short-term processes. Recent efforts [9][10][11] used data from moored buoys [13][14][15] but required air-sea fluxes, currents and temperature gradients from other data sources such as satellites and so could not resolve subsurface mixing, entrainment or other small scale ML processes. Additional uncertainty arises in the estimation of penetrative short wave radiation, when empirical formulas based on satellite-derived chlorophyll concentration are used 10,16 . The heat budget is often estimated over a constant depth, thus ignoring entrainment 12,17 . Further, estimates of the individual heat budget terms are usually not calculated simultaneously 10 , and a full closure of the equation is not usually feasible.

Bay of Bengal Boundary Layer experiment.
Here we use a groundbreaking approach (summarised schematically in Fig. 1) that combines high-resolution observations from ship-based and autonomous platforms to nearly close the surface ML heat budget and reveal the crucial importance of subsurface ocean processes for the evolution of SST in the Bay of Bengal (BoB), a region that has a prominent role in driving the Asian Summer Monsoon (ASM) 18 . Time series of oceanographic properties, including temperature, salinity, velocity, underwater radiation and subsurface mixing, along with surface fluxes of heat, were calculated from shipboard measurements onboard the RV Sindhu Sadhana in the southern BoB during the boreal summer monsoon of 2016 as part of the Bay of Bengal Boundary Layer Experiment (BoBBLE) 18 . Time-series measurements were made continuously for 11 days during 4-14 July at 89°E, 8°N (hereafter referred to as TSE; shown in Fig. 2A) using a CTD, uCTD, two ocean gliders, an Acoustic Doppler Current Profiler (ADCP), a Vertical Microstructure Profiler (VMP), an Automated Weather Station (AWS) and an underwater radiometer (see Methods). Figure 2B,C show, respectively, the locations and time of measurements of each instrument.
The horizontal advection term is a key process that is difficult to measure. To calculate this term, we adopted an Operation Advection (OA) strategy during BoBBLE, as follows. uCTD profiles were taken at intervals of approximately 1 km along westward and southward sections from TSE once per day, after dusk ( Fig. 2B,C). Additionally, two ocean gliders were deployed in an 'L' shaped configuration, approximately 15 km long, with the gliders at the ends of the 'L' and the ship at the corner. The 'L' was oriented so that the first glider was positioned upstream of the ship in the direction of the climatological (northeastward) current, with the second glider located perpendicular to this direction (Figs. 1 and 2B). Temperature profiles were obtained approximately every 2 hours from the ship CTD and every 3 hours from the gliders throughout the 11-day occupation of TSE ( Fig. 2C), allowing continuous and relatively high-frequency estimates of the horizontal temperature gradients. These were combined with current velocities from the shipboard ADCP to calculate the advection term. Independent estimates of the horizontal advection terms were obtained from the uCTD and gliders. Previous attempts with a uCTD or a CTD have involved surveying in a butterfly-shaped pattern 12,17 and took about 1.5 days to complete, during which time the ocean state would have changed significantly. The OA strategy employed during the Continental Tropical Convergence Zone (CTCZ) Programme in the BoB 19 , measured CTD profiles at four edges of a 'plus' with the time-series location as the center but lacked simultaneity and completeness in measuring all required parameters. The OA strategy that was adopted during BoBBLE took only one hour to complete each uCTD leg. This ensured independent simultaneous measurements of the horizontal temperature gradients and advection term once per day.
The 11-day time-series observations were carried out during a period of suppressed atmospheric convection with high surface solar radiation and the surface heat flux balance was mainly between the solar radiation and latent heat flux (Fig. 3A). Wind stress was moderate to weak (<0.2 Nm −2 ) with the maximum observed on 8 July (Fig. 3B). During July, the strong northeastward Southwest Monsoon Current 18,20 has been observed to advect cooler and saltier water to the TSE location (Fig. 3C,D) 21 . The near-surface currents were moderate (~0.35 m s −1 ) and flowed northeastward during the first three days of the observation and then reversed on 7 July (Fig. 3C). The zonal and meridional components of the currents were nearly uniform within the ML (Fig. 3G,H).   . This cooling and ML deepening phase did not coincide with the occurrence of a maximum in wind stress or with a net heat loss at the surface 21 . Hence, it cannot be explained directly by local surface forcing. Moreover, the cooling event occurred concurrently with the northeastward flow of near-surface currents. Hence, advection is likely to play a major role in this event. Following the cooling event, as the currents reversed and flowed southeastward, the SSS decreased, the ML shallowed and the sea surface warmed (Fig. 3C-F). The warming was superposed with a diurnal oscillation with an amplitude of ~0.1 °C 22 . By the end of the observations, the temperature, salinity and MLD were restored to their values at the beginning of the observations. closure of the heat budget. The sum of ML-averaged (see Methods) horizontal advection, vertical mixing, net surface heat flux (corrected for the loss due to penetrating short wave radiation) and entrainment reproduces the overall pattern of the tendency term (Fig. 4A,B) with a correlation of . 0 77 and root mean square difference of 0.22 ± 0.4 °C day −1 (within 95% confidence interval). This is a remarkably good closure of the ML heat budget, especially given the magnitude of small-scale variability, and is unprecedented for such short time scales.
Our analysis reveals that oceanic processes such as horizontal advection and entrainment are important in determining the net tendency of the temperature, as was evident on 5 and 9 July. On 5 July, when the negative tendency reached its peak (~−1 °C day −1 ), the horizontal advection (~−0.5 °C day −1 ) and entrainment (~−0.5 °C day −1 ) were the predominant contributors (Fig. 4A). The decomposition of these two processes reveals that the advection was in the zonal direction (Fig. 4C), towards the east, and the entrainment was caused mainly by the lateral induction (Fig. 4D) that arises from the horizontal advection across a sloping ML base. On 9 July, the horizontal advection was strongly positive, but was opposed by strong negative entrainment (largely due to lateral induction), resulting in a weak temperature tendency (Fig. 4A).
Net heat flux contributed to the warming that occurred during 8-9 July, as was expected from previous studies 23 . The diurnal variation that was observed in the tendency during the warming period (7-13 July) was driven www.nature.com/scientificreports www.nature.com/scientificreports/ by the diurnal variation in the net heat flux (Fig. 4A,E). As the ML shallowed, the penetrating short wave loss at the base of the ML increased. On 12 July, when the ML was shallower than 30 m, the loss at the base of the ML was about 20% of the radiation received at the sea surface (Fig. 4E). During these days, the chlorophyll fluorescence was low, indicating that the water was clear and allowed penetration of short-wave radiation below the ML (Fig. S1). The vertical (Fig. 4A) and horizontal mixing (Fig. 4F) were very weak throughout the observation period.
The temperature gradient was estimated using two different platforms, uCTD and gliders (Fig. 4G,H). Glider based estimates for ∂ ∂ T x / were consistently larger in magnitude, for both positive and negative values of ∂ ∂ T x / . Whereas, for ∂ ∂ T y / , there was no consistent bias. These differences represent real small-scale variability in spatial gradients. A map of high-resolution satellite SST, superimposed with the individual temperature measurements from CTD, uCTD and gliders (Fig. 5), shows that there is large small-scale spatial variability in SST gradients even at a scale of 20 km that explains the different measurements from these configurations. These differences underline the value of having two complementary observations of spatial gradients and the importance of high spatial resolution for these measurements. Besides, in situ measurements by uCTD and gliders were conducted only to the west and south of TSE (Fig. 1B)

Discussion
In summary, we have demonstrated that a novel combination of ship-based and autonomous platforms can effectively close the surface ML heat budget at unprecedented temporal resolution from in situ measurements alone. The time-series observations were carried out during a convectively suppressed phase of the boreal summer intraseasonal oscillation 18 , a cloud-free period, during which the sea surface is likely to warm due to intense solar radiation. However, our estimates revealed that SST in the BoB is significantly modulated by lateral advection and entrainment. During the initial stage of the observation, these oceanic processes inhibited the surface warming and thus delayed the transition towards a convectively active phase over the central BoB. Hence, within the BoB, the oceanic processes are key mechanisms that set the basin-scale SST distribution and influence convection over the region 24,25 , and our results underline that these processes cannot be neglected at any scale. Notably, the synoptic-scale convective systems that originate in the BoB during the active phase of the ASM propagate westward or northwestward and cause rainfall over the land 3,26 .
However, such small scale processes are difficult to resolve in coupled climate models and may thus contribute to biases in simulation of the monsoon 27 . Forecasting monsoon rainfall at intraseasonal to seasonal time scales requires accurate simulation of SST evolution 28 , and inadequate representation of these processes may also hinder simulations of the influence of climate change on monsoon rainfall 27 . The proof-of-concept study presented here (summarised in Fig. 1) serves as a blueprint for future observational campaigns that aim at determining ML energy budget over longer time scales. Such studies will help to identify sources of model biases in simulating SST and its variability, and thus aid improvements in model simulations of coupled ocean-atmosphere processes. This approach may be used for studying monsoon processes as here, and also for other phenomena such as: the ocean component of the Madden-Julian Oscillation 29 and interannual modes of variability such as El Niño-Southern Oscillation 30 and the Indian Ocean Dipole 31 .

instruments.
Temperature and salinity in the upper 500 m of the water column at the time-series location, TSE (8°N, 89°E), were repeatedly measured at approximately 3 hourly intervals using a CTD (SeaBird Electronics, SBE 9/11+) deployed from RV Sindhu Sadhana. We deployed two gliders, equipped with CT sensors, about 16 km southwest (7.9°N, 88.9°E) and southeast of TSE (7.9°N, 89.1°E), respectively. These gliders measured profiles from the surface to about 1000 m at a vertical resolution of 0.5-1 m. Temperature and salinity data from the gliders and CTD were optimally interpolated to a two-dimensional time-depth grid, with grid spacings of 1 hr and 1 m. Further details are in Matthews et al. 22 . However, the effective temporal resolution of these gridded data sets is approximately 2-3 hr, corresponding to the approximate time intervals between glider and ship successive CTD  profiles. Two vertical sections of temperature and salinity were measured, every day at around 1400 hours (UTC) and 1900 hours (UTC), respectively running southwards and westward of TSE, using a uCTD (Ocean Sciences-Teledyne underway CTD) towed by the ship. Each of these vertical sections was approximately 10 km long and took approximately one hour to complete. We measured temperature and salinity profiles using the uCTD while the ship was sailing at a steady speed of 6 knots. The uCTD probe was allowed to drop freely for 2 minutes at a vertical rate of 1.5-2.5 m s −1 from the surface to a depth of 180-300 m, which encompassed the ML. Ocean current velocity was measured using an ADCP with an operational frequency of 150 kHz. The data were processed 32 to derive current velocity in the top 11-180 m depth at a vertical resolution of 2 m and temporal resolution of 2 minutes. Data were corrected for mis-alignment angle and data collected during sudden accelerations were discarded. The data were then temporally averaged over a one-hour interval to be consistent with the CTD and glider data. To estimate vertical eddy diffusivity (K Z ) of temperature, a loosely tethered VMP (Rockland Scientific, VMP-250), equipped with two airfoil shear probes and a CT sensor, was used. Vertical profiles were measured to a maximum of 250 m, five times a day (0330, 0730, 1200, 1800 and 2330 hours UTC), using the VMP. An AWS, onboard the vessel, installed at an approximate height of 15 m above the sea surface, gave measurements of wind velocity (RM Young), atmospheric temperature (YSI), pressure (Honeywell), relative humidity (Rotronic), and all the components of the surface heat fluxes (LI-COR infra-red gas analyzer with 3D sonic anemometer-based eddy covariance system). To quantify the penetrative shortwave radiation, a hyperspectral underwater radiometer (Satlantic HyperProII) was used. The instrument assembly was equipped with sensors measuring downwelling irradiance (Ed), upwelling radiance (Lu), chlorophyll fluorescence, coloured dissolved organic matter fluorescence, back-scattering, conductivity and temperature. Three profiles were measured using the radiometer in the upper 150 m every day between 0600 and 0700 hours (UTC) at TSE.

Heat budget equation. The heat budget equation averaged in the space-time dependent oceanic ML, applying incommpressibility and Boussinesq approximations, is written as follows (See the Appendix of Moisan and
Niiler, 1998 33 for detailed derivation).
The terms in the Eq. 2, from left to right, are tendency (rate of change of temperature averaged over the MLD), horizontal advection, horizontal turbulent mixing, entrainment and net heat flux. The net heat flux comprises q 0 , the net heat flux at the ocean surface, and − q h , which represents the sum of penetrative loss of shortwave radiation (q pen ) and vertical turbulent mixing at the base of the ML. Moisan and Niiler, 1998 33 , had integrated the temperature equation from a constant isotherm located below the thermocline to the surface and hence, ignored the penetrative shortwave radiation and entrainment. However, these terms are not negligible, if the lower limit of integration is the base of the ML. We adopt the most widely used parameterization based on the second derivative of temperature for the vertical and horizontal turbulent mixing. The vertical turbulent mixing term, averaged in the ML, is given by, Here, κ Z is the vertical eddy diffusivity, estimated from the vertical microstructure measurements. We estimate the second derivative of the temperature using the uCTD data. Assuming a spatially and temporally invariant horizontal eddy diffusivity (κ H ), the horizontal turbulent mixing term, averaged in the ML is given by, Applying the Leibniz rule, www.nature.com/scientificreports www.nature.com/scientificreports/ A scaling of the three terms on the RHS of Eq. 6 suggests that the first term is higher than the other terms by an order of magnitude, and, therefore, we neglect the last two terms. Similarly, Substituting Eqs. 7 and 8 in 5, the horizontal mixing term averaged in the ML becomes, The entrainment term has three components, as given below.

ML tendency Lateral induction Vertical advection
The first term arises from the tendency (rate of change) of MLD. The second term results from horizontal advection across a sloping ML base and is often referred to as "lateral induction" 34 . The third term is the vertical advection.
Substituting Eqs. 4, 9 and 10 in Eq. 2 and re-arranging, we get,  21 , with the same MLD criteria, using the VMP data collected during BoBBLE shows that the total kinetic energy (TKE) dissipation rate was greater than 10 −7 W kg −1 in the ML, about two orders of magnitude higher than that was estimated just below the MLD. Further, they showed that a maxima in vertical shear estimated using the ADCP data coincided with the MLD. These results justify the MLD criteria chosen in this paper.
Tendency. The rate of change of temperature (tendency) at TSE was estimated from the CTD data. The data were first averaged in the ML and then the time derivative was estimated at hourly resolution using a centered differencing scheme. The evolution of SST was coherent with the ML temperature, but was marginally warmer by 0.06 °C on average, and exhibited stronger diurnal cycling (Fig. S2A). Note that we have deliberately chosen a definition of the ML that remains deeper than the diurnal warm layers, as these are not the focus of this study, so these differences are expected. The tendency of the SST (Fig. S2B) correlated well (0.90) with the tendency of the ML temperature (hereafter referred to as tendency).
Horizontal gradient from uCTD profiles. Any field such as temperature, T(x), varying in one dimension, can be expressed using the Taylor series expansion as follows, Applying a first order linear approximation, we can rewrite Eq. 12 as, x We made several estimates of ΔT utilizing temperature profiles measured by uCTD spread at different locations along a linear transect (see the schematic diagram in Fig. S3). Assuming that there is an offset of δ and a Scientific RepoRtS | (2020) 10:7062 | https://doi.org/10.1038/s41598-020-63320-0 www.nature.com/scientificreports www.nature.com/scientificreports/ random variability of ε x ( ), the observed temperature difference between the two points can be written as follows.
where the offset δ is due to sampling variability and is typically very small (Figs. S4 and S5). The ε x ( ) accounts for random variability and can be identified with the higher order terms O[Δx 2 ] from Eq. 13. Here, suffix i refers to the i th uCTD profile along the linear transect and suffix 1 denotes the first profile, closest to TSE. A similar expression can be written for the meridional gradient of temperature in the y-direction. The observed ΔT obs were then plotted as a function of Δx using the westward uCTD section (Fig. S4) and Δy using the southward section (Fig. S5) for each day. Straight lines were then fitted assuming that the scatter satisfies the first order approximation given in Eq. 13. The slope of the line (Eq. 14) gives the temperature gradient and the intercept gives the offset (δ). Any scatter about the straight line will be due to random variability, including the higher order terms. This method of linear fitting assumes that the horizontal temperature gradient has a spatial scale of about 10 km and therefore any variability whose length scale is smaller than 10 km is an error. The temperature gradient along the meridional direction and the errors associated were also estimated in a similar way. We find that the data in Figs. S4 and S5 do lie approximately along straight lines, and conclude that the linear approximation made in Eq. 13 is valid. Here, we could get the estimates of the uncertainties involved in the gradient calculation, which is not possible with glider based gradient estimates. The offsets on 5 July in the estimate of dT dx / and dT dy / were −0.04 and −0.004 °C and the random variability were 0.03 and 0.01 °C. The offsets were removed while presenting the final estimate of the gradients.
Horizontal advection from glider profiles. We positioned the ship and two gliders in an L-shaped configuration (see Fig. S6). As the gliders are not stationary, their positions were different each time a profile was measured. Therefore, the temperature (averaged in ML) difference between a glider and CTD is given by  If the observations were aligned on the x and y axes, Eq. 16 simplifies to the familiar finite difference estimates. The averaged velocity components were then multiplied to the respective temperature gradients to obtain the components of horizontal temperature advection. Before averaging the ADCP-derived velocity components in the ML, it was assumed that the top 11 m (depth of the first bin of ADCP measurement) of the water column has a velocity equal to the velocity at 11 m. Alternative extrapolation methods were tested but found to be unrealistic.
Horizontal mixing. We utilize the uCTD data to estimate the second horizontal derivative of ML-averaged temperature. This method is similar to the method followed to compute the horizontal gradient from uCTD profiles, with difference in temperature replaced by a first derivative estimated using a forward difference scheme. We made several estimates of ΔT/Δx utilizing temperature profiles measured by uCTD spread at different locations along the x-axis. Then ΔT/Δy was plotted as a function of Δx (similar to Eq. 14).
obs uCTD uCTD The second derivative of temperature (∂ ∂ T x / 2 2 ) is the slope of the line that can be fitted to the observed scatter of ΔT/Δx versus Δx. δ, the intercept of the fit is the offset and ε x ( ), measure of the scatter, is the random variability. The observed ΔT/Δx and ΔT/Δy were then plotted as functions of Δx and Δy, respectively (Figs. S7 and S8) leading to an estimate of ∂ ∂ T x / 2 2 of order 10 −9 °C m −2 . κ H is a function of length scale 37 . For a length scale of 10 km, the order of magnitude of κ H is about 10 m 2 s −1 . The horizontal mixing term is then of order 10 −3 °C day −1 and has a negligible contribution compared to horizontal advection (Fig. 4F). Earlier studies also suggest that this term is very small [38][39][40] .
Vertical mixing. To estimate the vertical diffusion coefficient using the VMP data, we followed the methodology described in George et al. 21 . Turbulent kinetic energy dissipation rate (ε) was estimated from the velocity shears measured by the VMP following Roget et al. 41 assuming isotopic turbulence. The vertical diffusion coefficient was then calculated using the relationship 42 , K Z = Γε N / 2 . The mixing efficiency (Γ) was taken to be a constant (0.2) and the Brunt-Väisälä frequency (N 2 ) was estimated using the CTD data. The finite difference decomposition of the vertical mixing term (Eq. 4) was carried out using a forward difference scheme applied at the base of the ML.
Vertical advection was estimated as the tendency of 27 °C isotherm located in the thermocline, assuming that the base of the ML oscillates at the same rate as the thermocline 43 . The assumption here is that entrainment occurs when there is a difference in the vertical movement of the ML and the vertical movement of the thermocline. By adding the ML tendency and the vertical advection term (which have opposite signs), we compute the difference in vertical motion of these two surfaces. We assume that this difference in vertical motion reflects the entrainment of water into the ML, since adiabatic waves will move both surfaces in concert. Internal waves that stretch or compress both the thermocline and the ML do not lead to a transport of heat through the base of the ML and thus do not directly influence ML temperature. In short, the assumption is that the ocean is adiabatic in the interior and the water parcels do not cross the pycnocline.
Net surface heat flux. Estimation of sensible and turbulent latent heat fluxes were done using the eddy covariance method [44][45][46] .
Penetrating short wave radiation. q pen was estimated based on Lotliker et al. 47 . The diffuse attenuation coefficient of photosynthetically available radiation (PAR) at wavelengths between 400-700 nm (k PAR ) was estimated using the following equation. is the spectral downwelling irradiance at depth z and wavelength λ measured using the hyperspectral profiling radiometer. λ + E ( 0 , ) d represents the spectral downwelling irradiance just above the sea surface and was calculated using the Fresnel reflection albedo (α = . 0 043) for irradiance from sun and sky as given below.
represents the spectral downwelling irradiance just below the sea surface. k PAR varied between 0.0635 to 0.0759 m −1 during the time-series observation at TSE. The daily values of k PAR were then used to fit exponential curves and were compared with downwelling irradiance averaged over the visible spectrum during 4-13 July 2016 at TSE (see Fig. S9). Finally, the penetrative short wave loss at the base of the ML was estimated as a function of the net shortwave radiation received at the sea surface. It was found that the estimated q pen inversely correlated (−0.81 with 99% significance) with the variation of chlorophyll in the ML (Fig. S1).