Surface net heat flux estimated from drifter observations

The surface satellite-tracked drifter dataset (Lumpkin and Centurioni 2019) covering the period 1988–2013 is used to estimate terms in the mixed layer temperature equation for the global ocean. The dataset consists of latitude and longitude arrays of mixed layer temperature and current velocity components with 6-hourly temporal resolution, obtained by the procedure of optimum interpolation. During the period 1988–2013, over 17 000 drifters have been deployed globally (figure S1), and more than 24.7 million individual measurements are available (figure S2). These drifters have submerged 'holy sock' drogues centered at a depth of 15 m to reduce the large biases induced by processes at air-sea interface (e.g. winds and breaking waves) and thereby accurately measure mixed layer temperature and current velocity, even under the most intense TCs. Therefore, these globally distributed drifters allow for the possibility of accurately estimating the heat exchange across the air-sea interface from the mixed layer temperature equation.

The thickness of the mixed layer is a key factor for determination of the mixed layer temperature. The MLD used in the present study are determined by a temperature criterion from the Simple Ocean Data Assimilation version 3.3.1 monthly mean ocean reanalysis dataset, with a spatial resolution of 0.5° × 0.5° spanning from 1988 to 2013.

Satellite-derived net heat flux data with daily temporal resolution obtained from the Japanese Ocean Flux Data Sets with Use of Remote Sensing Observations (J-OFURO) research project based on observational information obtained from various earth observing satellites is used to compare with the drifter-derived heat flux. In the present study, the third-generation dataset, J-OFURO3, is adopted. The satellite data, which comes with a horizontal resolution of 0.25° × 0.25°, is interpolated onto a 0.5° × 0.5° grid to be consistent with the MLD data during 1988–2013.

Observations of TC occurrence during 1988–2013 were acquired from the best track data provided by the Joint Typhoon Warning Center for the Western Pacific Ocean, the Indian Ocean, and the Southern Hemisphere, and the National Hurricane Center and Central Pacific Hurricane Center for the Atlantic and Northeast and Central Pacific Oceans. This TC track data is produced with 6-hourly temporal resolution.

Buoy-derived net heat flux data with daily temporal resolution obtained from the Global Tropical Moored Buoy Array (GTMBA) and Ocean Climate Stations (OCSs) is used to compare with the drifter-derived and satellite-derived heat fluxes. In the present analysis, 36 buoys located in the global tropical oceans from the GTMBA are used in climatology, and 8 buoys from the GTMBA and 1 buoy from the OCS are applied under TC conditions.

In order to obtain the heat flux across the air-sea interface, we solve the mixed layer temperature equation (Dong et al 2007), which can be expressed as:

$$\begin{equation}\frac{Q}{{\rho {c_{\text{p}}}h}} + \frac{{ - \Delta T}}{h}{w_{\text{e}}} + k{\nabla ^2}T = \frac{{\partial T}}{{\partial t}} + v \cdot \nabla T,\end{equation} \tag{ 1 }$$

where T is the mixed layer temperature, v is the mixed layer current velocity including both the zonal and meridional components, ρ (1027 kg m⁻³) is the reference density of seawater, c_p (4000 J kg⁻¹ K⁻¹) is the specific heat of seawater at constant pressure, h is the MLD, and ΔT is the temperature difference across the bottom of the mixed layer (set to be 0.5 °C). Here Q is the sum of the surface net heat flux that is positive into the ocean and the downward radiative heat flux at the bottom of the mixed layer, w_e is the entrainment velocity, and k is the eddy diffusivity.

Equation (1) is a nonhomogeneous linear equation, where Q, w_e and k are considered as unknown variables, with $\frac{1}{{\rho {c_{\text{p}}}h}}$, $\frac{{ - \Delta T}}{h}$ and ${\nabla ^2}T$ being the corresponding coefficients. The three coefficients on the left-hand side and temperature tendency and horizontal advection terms on the right-hand side can be determined from the drifter measurements, MLD data and coefficient constants mentioned above. To estimate Q, w_e and k by solving equation (1), the major steps are as follows:

Step 1, estimate partial derivatives of $\frac{{\partial T}}{{\partial t}},\frac{{\partial T}}{{\partial x}},\frac{{\partial T}}{{\partial y}},\frac{{{\partial ^2}T}}{{\partial {x^2}}},\frac{{{\partial ^2}T}}{{\partial {y^2}}}$ using the central difference method based on each set of five consecutive drifter measurements along a trajectory (figure S3).

Step 2, estimate Q, w_e and k at the central location (location 4 in figure S3) by solving the following three equations derived from equation (1) with three sets of the above partial derivative estimates. Note we suppose Q, w_e and k remain unchanged during the observed period of 1.5 day

$$\begin{align*}{\frac{Q}{{\rho {c_{\text{p}}}h \left. \right|}}_3} + \frac{{ - \Delta T}}{{{{\left. h \right|}_3}}}{w_{\text{e}}} + k{\left. {{\nabla ^2}T} \right|_3} &= {\left. {\frac{{\partial T}}{{\partial t}}} \right|_3} + {\left. v \right|_3} \cdot {\left. {\nabla T} \right|_3}{\text{ ,}}\nonumber\\ {\frac{Q}{{\rho {c_{\text{p}}}h \left. \right|}}_4} + \frac{{ - \Delta T}}{{{{\left. h \right|}_4}}}{w_{\text{e}}} + k{\left. {{\nabla ^2}T} \right|_4} &= {\left. {\frac{{\partial T}}{{\partial t}}} \right|_4} + {\left. v \right|_4} \cdot {\left. {\nabla T} \right|_{4\,}},\nonumber\\ {\frac{Q}{{\rho {c_{\text{p}}}h \left. \right|}}_5} + \frac{{ - \Delta T}}{{{{\left. h \right|}_5}}}{w_{\text{e}}} + k{\left. {{\nabla ^2}T} \right|_5} &= {\left. {\frac{{\partial T}}{{\partial t}}} \right|_5} + {\left. v \right|_5} \cdot {\left. {\nabla T} \right|_5}{ } .\end{align*}$$

Step 3, repeat the above processes along each drifter trajectory. Estimated Q, w_e and k from all drifters are further grouped into bins of 0.5° × 0.5°globally. Here the 0.5° × 0.5° bins are adopted to be consistent with the spatial resolution of the MLD.

Step 4, the same as described in step 3 but under TC conditions. The estimated Q, w_e and k under TCs are grouped into bins of 5° × 2.5° globally. Because the average size of a TC is around 500 km in diameter, 5° in zonal direction is used. Half of this distance is applied in the meridional direction due to less coherence of ocean variability. Specifically, we also grouped these estimated Q, w_e and k in the along-track and cross-track coordinate system under TCs.

Global distributions of the climatological entrainment velocity, eddy diffusivity and net heat flux during the study period (1988–2013) are obtained. To reduce the noise from mesoscale features and focus on its larger-scale pattern, these climatologic values are smoothed in space using a mean filter with the smoothing lengths of 1° in latitude and 2° in longitude.

Large values of the entrainment velocity appear in regions of strong currents such as near the Kuroshio Current, the Gulf Stream, and the Antarctic Circumpolar Current (ACC) (figure 1(a)). High rates of mixed layer deepening may contribute to the strong entrainment in the Western Boundary Current (WBC) regions, such as the Kuroshio Current (Qu 2003), while the strong prevailing westerlies and the associated intense upwelling may contribute to the strong entrainment in the ACC regions. On average, the long-term mean entrainment velocity is about 6.2 × 10⁻⁴m s⁻¹ on a global scale, which is on the same order of magnitude as that found by Nagai et al (2005). High values of the eddy diffusivity appear in the western boundaries around the global ocean, the ACC region, as well as along the global equatorial belts (figure 1(b)). These large values are likely to correspond with the generation of mesoscale eddies, baroclinic instability, and tropical instability waves, respectively (Zhurbas and Oh 2004, Zhurbas et al 2014). In general, the spatial pattern of the eddy diffusivity bears much resemblance to that obtained by Zhurbas et al (2014). The global distributions of both the entrainment velocity and eddy diffusivity verify that the mixed layer temperature equation has the potential to estimate the surface net heat flux accurately from the drifter dataset.

Zoom In Zoom Out Reset image size

The drifter-derived net heat flux in climatology is shown in figure 2(a). The ocean gains heat primarily in the tropics, especially in the equatorial cold tongue regions, and releases heat primarily at high latitudes especially in the significant WBC systems such as the Kuroshio Current and Gulf Stream in the Northern Hemisphere and the Agulhas Current in the Southern Hemisphere. The high-latitude North Atlantic Ocean appears to be another region of significant heat transfer into the atmosphere from the underlying ocean, which is likely related to the deep-water formation there where heat is released to make surface water dense enough to sink. In addition, heat absorption into the ocean also occurs in the four major eastern boundary upwelling regions; namely the California and Peru-Chile upwelling regions in the Pacific Ocean, the Canary and Benguela upwelling regions in the Atlantic Ocean, and the Somali current upwelling region on the western boundary of the Indian Ocean. The upwelling in these regions results in low sea surface temperatures (SSTs) that lead to stronger heat absorption from the atmosphere.

Zoom In Zoom Out Reset image size

To see whether the spatial pattern of the drifter-derived heat flux is similar with that of the satellite-derived heat flux, figure 2 shows the globally distributed Q_net derived from drifters (figure 2(a)) and from satellites using the J-OFURO3 flux product (figure 2(b)). The spatial structure of the drifter-derived heat flux bears considerable resemblance to the satellite-derived heat flux and has a high spatial correlation coefficient of 0.83, which is statistically significant at the 99% confidence level. However, the drifter-derived heat flux is much smaller in magnitude than the satellite-derived heat flux over nearly the entire global oceans, especially for the EC and WBC regions (figure 2(c)). Quantitatively, the drifter-derived heat fluxes have a global mean positive value of 16.8 W m⁻² and negative value of −14.4 W m⁻², while the associated satellite-derived values are 36.3 W m⁻² and −28.7 W m⁻², respectively. Consistently, the spatially averaged drifter-derived heat flux in the regions of Kuroshio Current, Gulf Stream, Agulhas Current, Equatorial Current and ACC are −13.8, −25.3, −7.0, 18.2, and 5.7 W m⁻², respectively, all of which are weaker than the corresponding satellite observations of −23.7, −43.3, −18.3, 48.3 and 16.8 W m⁻². In addition, both the two heat flux estimates show similar seasonal variations with heat released to the atmosphere in each hemisphere's winter and transferred into the ocean in each hemisphere's summer at all latitudes, but with different magnitudes.

To further compare the global heat budget balance (Isemer et al 1989, Josey et al 1999, Schuckmann et al 2016, Liu et al 2017, Valdivieso et al 2017) between the drifter-derived and satellite-derived heat fluxes, the global mean Q_net is calculated based on the following area-weighted average equation:

$$\begin{equation}\overline {{Q_{{\text{net}}}}} = \frac{{\sum\limits_{i,j} {Q_{{\text{net}}}}\left( {i,j} \right) \cdot \cos \left[ {{\text{lat}}(j)} \right]}}{{\sum\limits_{i,j} \cos \left[ {{\text{lat}}(j)} \right]}}.\end{equation} \tag{ 2 }$$

Results show that the globally averaged long-term mean of Q_net for the drifter-derived and satellite-derived heat fluxes (averaged in areas where there are no missing data for either of the two fluxes) are 3.9 W m⁻² and 23.5 W m⁻², respectively. Obviously, the drifter-derived heat flux is closer to the physical constraint on the global mean Q_net of 2–3 W m⁻² (Serreze et al 2007, Bengtsson et al 2013), which shows its better performance in closing the heat budget at the ocean surface.

Because in situ observations are difficult and rare during the passage of a TC, the composite of Q_net is computed by using drifter observations within 7R_max of TC tracks in the TC-coordinate system. The left-hand panels of figures 3(a)–(c) show the composite Q_net estimates computed for all TCs, tropical storms (TSs) and category 1–5 TCs, respectively. These TC related air-sea heat fluxes all have an obvious strong left-to-right asymmetry with the maximum values located on the right side of TC tracks. Such left-to-right asymmetry primarily arises from the asymmetric distributions of TC wind structure, with the strongest wind stress occurring on the right side in the Northern Hemisphere and left side in the Southern Hemisphere of TC tracks (Price 1981), respectively. It should be noted that in this analysis, the left-to-right asymmetric spatial pattern in the Southern Hemisphere has been adjusted to be consistent with that in the Northern Hemisphere. Specifically, the maximum heat release increases from 80 W m⁻² to 175 W m⁻² as TCs intensify from TSs to category 1–5 TCs. Detailed comparisons of the drifter-derived net heat fluxes under various intensity TCs are listed in table 1.

Zoom In Zoom Out Reset image size

Further, we extend the extreme air-sea heat fluxes occurring under TC conditions to the main TC basins. Figure 4(a) shows the spatial distribution of the drifter-derived heat fluxes resulting from the passage of TCs. In general, heat is released from the underlying ocean surface with the magnitudes decreasing poleward, which may be correlated with the spatial distribution of the SSTs that almost decrease linearly with the increasing latitude. The satellite-derived heat flux shows a similar spatial pattern but with smaller magnitudes (figure 4(b)), and there are regions in midlatitudes of both the two hemispheres where heat is even released from the atmosphere to the ocean. Specifically, the spatially averaged mean Q_net are −186, −124, −112, and −85 W m⁻² at 5° N, 10° N, 20° N, and 30° N, respectively, all of which are larger in magnitude than the corresponding satellite-derived results of −89, −63, −29, and −21 W m⁻². This suggests that the drifter-derived heat flux estimates as a result of the passage of TCs is 2–4 times those derived from satellite data. Because satellite measurements are strongly biased by drop-outs from cloud coverage, the accuracy of the satellite-derived heat flux under TCs may be seriously degraded. However, the drifter measurements are of high accuracy, with the use of drogues centered at a depth of 15 m to reduce the large bias induced by wave processes at the air-sea interface and the Global Position System in localization, even under TCs. Therefore, the drifter-derived heat fluxes should be more realistic than the satellite-derived heat fluxes under TC conditions.

Zoom In Zoom Out Reset image size

Furthermore, taking the heat flux obtained from 8 buoys of the GTMBA and 1 buoy of the OCS as a subset for evaluation, we further compare the drifter-derived and satellite-derived heat fluxes under TCs. As shown in figure S4, the errors of the drifter-derived heat flux are much smaller compared with those of the satellite-derived heat flux, with the mean errors being 12.3 and 128.9 W m⁻², respectively. Additionally, as listed in table S1, the mean absolute errors, the root mean-square errors and the percent biases of these two heat fluxes further evidence that the drifter-derived heat flux performs better than the satellite-derived heat flux under TCs. These can also be seen in the normalized Taylor diagram presenting a comparison of the buoy-derived heat flux with the corresponding drifter-derived and satellite-derived heat fluxes under TC conditions (figure S5).