Accelerated estimation of sea-spray-mediated heat flux using Gaussian quadrature: case studies with a coupled CFSv2.0-WW3 system

Sea-spray-mediated heat flux plays an important role in air–sea heat transfer. Heat flux integrated over the droplet size spectrum can simulate well the total heat flux induced by sea spray droplets. Previously, a fast algorithm of spray flux assuming single-radius droplets (A15) was widely used, as the full-size spectrum integral is computationally expensive. Based on the Gaussian quadrature (GQ) method, a new fast algorithm (SPRAY-GQ) of sea-spray-mediated heat flux is derived. The performance of SPRAY-GQ is evaluated by comparing heat fluxes with those estimated from the widely used A15. The new algorithm shows a better agreement with the original spectrum integral. To further evaluate the numerical errors of A15 and SPRAY-GQ, the two algorithms are implemented into the coupled Climate Forecast System model version 2.0 (CFSv2.0) and WAVEWATCH III (WW3) system, and a series of 56 d simulations in summer and winter are conducted and compared. The comparisons with satellite measurements and reanalysis data show that the SPRAY-GQ algorithm could lead to more reasonable simulation than the A15 algorithm by modifying air–sea heat flux. For experiments based on SPRAY-GQ, the sea surface temperature at middle to high latitudes of both hemispheres, particularly in summer, is significantly improved compared with the experiments based on A15. The simulation of 10 m wind speed and significant wave height at middle to low latitudes of the Northern Hemisphere after the first 2 weeks is improved as well. These improvements are due to the reduced numerical errors. The computational time of SPRAY-GQ is about the same as that of A15. Therefore, the newly developed SPRAY-GQ algorithm has potential to be used for the calculation of spray-mediated heat flux in coupled models.

1 Introduction

Sea spray droplets, ejected from oceans, include film drops, jet drops and spume drops (Veron, 2015). The first two types of droplets are generated from bubble bursting caused by ocean surface wave breaking, with radius ranging from 0.5 to 50 µm (Resch and Afeti, 1991; Thorpe, 1992; Melville, 1996; Spiel, 1997; Andreas, 1998; Lhuissier and Villermaux, 2012). Spume drops are generated by strong winds (> 7–11 m s⁻¹) which directly tear the wave crests, with larger radius ranging from tens to hundreds of micrometers (Koga, 1981; Andreas et al., 1995; Andreas, 1998). Sea spray droplets play an important role in weather and climate processes (Fox-Kemper et al., 2022). On one hand, sea spray droplets contribute to local marine aerosols and subsequently modify the local radiation balance (Fairall et al., 1983; Burk, 1984; Fairall and Larsen, 1984). On the other hand, sea spray droplets affect the fluxes of heat, momentum, salt, and freshwater between the atmosphere and ocean (Andreas, 1992; Andreas et al., 2008; Andreas, 2010; Andreas et al., 2015; Ling and Kao, 1976; Fairall et al., 1994; Andreas and Decosmo, 2002).

The sea-spray-mediated heat transfer mainly occurs within the droplet evaporation layer (DEL) near the sea surface (Andreas and Decosmo, 1999, 2002; Fairall et al., 1994). Sea spray droplets with the same temperature as the ocean surface can lead to sensible heat flux in DEL, while water evaporated from these droplets can further release latent heat to the atmosphere (Andreas, 1992; Borisenkov, 1974; Bortkovskii, 1973; Wu, 1974; Monahan and Van Patten, 1988; Ling and Kao, 1976). Part of the sea-spray-mediated sensible heat is absorbed by droplet evaporation, which further increases the air–sea temperature difference and thus increases the sea-spray-mediated sensible heat flux (Fairall et al., 1994; Andreas and Decosmo, 2002). Since strong winds produce more sea spray droplets with a larger radius, sea-spray-mediated heat fluxes increase with wind speed (Fairall et al., 1994) and contribute to more than 10 % of the total surface heat flux after reaching the threshold speed (> 11 m s⁻¹ for sensible heat flux and > 13 m s⁻¹ for latent heat flux) (Andreas et al., 2008). In addition, when a droplet is released into the air, it is accelerated due to surface winds (Edson and Andreas, 1997; Fairall et al., 1994; Van Eijk et al., 2011; Wu et al., 2017). If the droplet could fall back into the ocean, additional momentum would be injected into the ocean from the atmosphere (Andreas, 1992, 2004).

The usual bulk parameterizations in numerical models for surface fluxes only include the interfacial (turbulent) fluxes (e.g., Fairall et al., 1996) while neglecting the significant contributions of sea spray droplets in DEL (Andreas et al., 2008; Fairall et al., 1994; Smith, 1997; Emanuel, 1995). Andreas and Emanuel (2001) implemented sea-spray-mediated heat flux and momentum flux parameterizations into a simple tropical cyclone model and found that the sea-spray-mediated heat flux can significantly enhance tropical cyclone intensity. The similar enhancement of tropical cyclone intensity was also noticed in recent regional coupling systems by including sea-spray-mediated heat flux (Xu et al., 2021a; Liu et al., 2012; Garg et al., 2018; Zhao et al., 2017). In the First Institute of Oceanography Earth System Model, Bao et al. (2020) first incorporated the sea-spray-mediated heat flux in global climate simulation. Following Bao et al. (2020), Song et al. (2022) found that the sea-spray-mediated heat flux can lead to cooling at the air–sea interface and westerlies strengthening in the Southern Ocean, and it thus improves estimates of sea surface temperature (SST).

Since the parameterization of sea-spray-mediated heat flux derived from observations requires full-size spectral integral and thus is computationally expensive for large-scale models (Table 1, details in Sect. 4.2; Andreas, 1989, 1990, 1992; Andreas et al., 2015), a simplified algorithm based on a single radius of sea spray droplets (Andreas et al., 2015, 2008) is widely used in atmosphere–ocean coupling systems (Xu et al., 2021a; Liu et al., 2012; Garg et al., 2018; Zhao et al., 2017; Song et al., 2022; Bao et al., 2020) and is apt to produce numerical errors. To reduce these numerical errors induced by the single radius of sea spray droplets, we develop a new fast algorithm of sea-spray-mediated heat flux based on the Gaussian quadrature (GQ) method, a fast and accurate way to calculate spectral integral. The GQ method has been successfully used for the estimation of domain-averaged radiative flux profiles (Li and Barker, 2018). The performance of the GQ-based fast algorithm of the sea-spray-mediated heat flux is evaluated and compared with the simplified algorithm for single radius of Andreas et al. (2015), referred to as A15 hereafter. The results are first compared with the original parameterization using full-size spectral integral (A92, hereafter). Then, the parameterizations with different algorithms are implemented in a global coupled atmosphere–ocean–wave system (Shi et al., 2022), and the results are compared with global satellite measurements and reanalysis data.

Table 1 The run time of the coupled Climate Forecast System model version 2.0 (CFSv2.0) and WAVEWATCH III (WW3) system's global experiments for 7 d forecast with different parameterizations.

The rest of the paper is structured as follows: observation and reanalysis data for comparisons are introduced in Sect. 2, the derivation of the GQ-based fast algorithm and the global coupling system are described in Sect. 3, the performance of the new fast algorithm is evaluated in Sect. 4, and finally, a summary and discussion are given in Sect. 5.

2 Data

The fifth-generation European Centre for Medium-range Weather Forecasts (ECMWF) Reanalysis (ERA5; Hersbach et al., 2020) 10 m wind speed (WSP10), 2 m air temperature (T02), 2 m dew point temperature, surface pressure, and significant wave height (SWH) with a spatial resolution of 0.5^∘ are used. Additionally, WSP10, T02 and 2 m specific humidity (SPH) data from the Objectively Analyzed air–sea Fluxes (OAFlux) products (Yu et al., 2008) are also applied for comparison with 1^∘ × 1^∘ resolution. The daily average satellite Optimum Interpolation SST (OISST) data are obtained from the National Oceanic and Atmospheric Administration (NOAA) with a spatial resolution of 0.25^∘ (Reynolds et al., 2007). The global monthly mean salinity observations from the European Space Agency (ESA; https://climate.esa.int/sites/default/files/SSS_cci-D1.1-URD-v1r4_signed-accepted.pdf, last access: 18 March 2023) are applied. Besides these, we also use the monthly global ocean Remote Sensing Systems (RSS) satellite data products for WSP10 (https://data.remss.com/wind/monthly_1deg/, last access: 18 March 2023) and the reprocessed L4 satellite measurements for SWH (https://doi.org/10.48670/moi-00177, Collecte Localisation Satellites, 2021), to validate the simulation results and ERA5 data.

3 Methods

3.1 Development of a fast algorithm based on GQ

The effects of sea spray droplets on sensible and latent heat fluxes (H_S,SP, H_L,SP) contribute to the total turbulent sensible and latent heat fluxes (H_S,T, H_L,T) at the air–sea interface. That is,

$$\begin{array}{}\text{(1)}& {\displaystyle }{H}_{\mathrm{S},\mathrm{T}}=H_{\mathrm{S}}+H_{\mathrm{S},\mathrm{SP}}\text{(2)}& {\displaystyle }{H}_{\mathrm{L},\mathrm{T}}=H_{\mathrm{L}}+H_{\mathrm{L},\mathrm{SP}},\end{array}$$

where H_S and H_L are the sensible and latent heat fluxes at the air–sea interface due to the air–sea differences of temperature and humidity. Based on observations of total turbulent heat fluxes and the Coupled Ocean–Atmosphere Response Experiment (COARE) algorithm (Andreas et al., 2015; Fairall et al., 1996), A92 integrates the sea-spray-mediated sensible and latent heat flux spectrums over initial droplet radius (Q_S(r₀) and Q_L(r₀)) to estimate H_S,SP and H_L,SP (details in Appendix A; Andreas, 1989, 1990, 1992; Andreas and Decosmo, 2002). The distributions of Q_S(r₀) and Q_L(r₀) spectrums as functions of initial droplet radius r₀ under various atmosphere and ocean states are shown in Fig. 1, indicating that Q_S and Q_L spectrums are more sensitive to the change of WSP10 and less sensitive to other variables, including T02, 2 m relative humidity, SST, surface air pressure and sea surface salinity.

Figure 1 The radius-specific sea-spray-mediated sensible (Q_S: black) and latent (Q_L: red) heat fluxes as functions of initial radius r₀: U₁₀, Ta, RH, Tw, P and S are 10 m wind speed, 2 m air temperature, 2 m relative humidity, sea surface temperature, surface air pressure and surface salinity, respectively.

Since the calculation of H_S,SP and H_L,SP in A92 is computationally expensive due to full-size spectral integral (Eqs. A5–A6 of Appendix A), it is difficult to apply A92 directly in coupled modeling systems. A15 (Andreas et al., 2015) developed a fast algorithm by using a single representative droplet radius (details in Appendix B), which was widely adopted in recent regional and global coupling systems (Xu et al., 2021a; Liu et al., 2012; Garg et al., 2018; Zhao et al., 2017; Song et al., 2022; Bao et al., 2020). In this study, we apply a three-node GQ method (details in Appendix C) to develop a new fast algorithm to approximate the full-size spectral integral of A92. Notably, GQ can converge exponentially to the actual integral only for a smooth function, which is a prerequisite for GQ (McClarren, 2018). As functions of r₀, Q_S(r₀) and Q_L(r₀) are not smooth (Fig. 1), so data sorting from largest to smallest is required. After sorting, local Q_S(r₀) and Q_L(r₀) become Q_{S_sort}(m) and Q_{L_sort}(m), and then GQ can be used to estimate the integral of Q_{S_sort}(m) and Q_{L_sort}(m). Note that the independent variable m is not equivalent to the original r₀ but only indicates the position. In this way, according to Appendix C, m₁=443, m₂=251 and m₃=58 are three GQ nodes of Q_{S_sort}(m) and Q_{L_sort}(m), and we can get the corresponding r₀ for local Q_S(Q_L), denoted as r_S1(r_L1), r_S2(r_L2) and r_S3(r_L3). However, the sorting leads to high complexity of GQ comparable to A92, and the values of r_S1(r_L1), r_S2(r_L2) and r_S3(r_L3) vary under various atmosphere and ocean states in the globe. Therefore, it is necessary to find the general approximate values of r_S1(r_L1), r_S2(r_L2) and r_S3(r_L3) via global statistical analyses, to avoid the sorting during application.

To derive the general approximate values of r_S1(r_L1), r_S2(r_L2) and r_S3(r_L3), we calculate the distribution of the sea-spray-mediated heat flux spectral following A92, based on the global daily WSP10, T02, 2 m dew point temperature, surface pressure and SWH of ERA5 and OISST from 1 to 31 August 2018. Since the sea-spray-mediated heat flux is not sensitive to salinity (Fig. 1e and f) and only monthly observational data are available, the ESA monthly salinity is applied. From the global spectrums, we sort Q_S and Q_L from largest to smallest to obtain local r_S1, r_S2 and r_S3 (r_L1, r_L2 and r_L3) for every grid point, whose global distribution of occurrence frequency in percentage is shown in Fig. 2. It is noted that except for r_L3, all other five nodes have frequency roughly concentrated at a constant (peak frequency > 65 % in Fig. 2a, b, d–f; Eqs. 3 and 4), while for r_L3, there is a 92.53 % concentration between 55 and 90 µm (Fig. 2c). Then we found that r_L3 (55–90 µm) is related to WSP10 (Fig. S1 in the Supplement): thereby, we set the approximate values as

$$\begin{array}{}\text{(3)}& {\displaystyle }{r}_{\mathrm{S}\mathrm{1}}=\mathrm{459.056},r_{\mathrm{S}\mathrm{2}}=\mathrm{294.185},r_{\mathrm{S}\mathrm{3}}=\mathrm{166.771},\text{(4)}& {\displaystyle }{r}_{\mathrm{L}\mathrm{1}}=\mathrm{443.914},r_{\mathrm{L}\mathrm{2}}=\mathrm{251.0498},\text{(5)}& {\displaystyle }{r}_{\mathrm{L}\mathrm{3}}=\left\{\begin{array}{l}\mathrm{60.310}\mathrm{WSP}{\mathrm{10}}^{\mathrm{0.1161}},\mathrm{WSP}\mathrm{10}\ge \mathrm{2}{\mathrm{m}\mathrm{s}}^{-\mathrm{1}}\\ \mathrm{58.086},\text{WSP10}<\mathrm{2}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}\end{array}\right.,\end{array}$$

where the unit of the radius is in micrometers. Afterwards, we directly use Eqs. (3)–(5) to approximate the full-size spectral integral of A92 without sorting as

$$\begin{array}{}\text{(6)}& {\displaystyle }\underset{a}{\overset{b}{\int }}{Q}_{\mathrm{S}}\left(r_{\mathrm{0}}\right)\mathrm{d}{r}_{\mathrm{0}}\approx {\displaystyle \frac{b-a}{\mathrm{2}}}\sum _{i=\mathrm{1}}^{\mathrm{3}}{\mathit{\omega }}_i{Q}_{\mathrm{S}}\left(r_{\mathrm{S}i}\right)\text{(7)}& {\displaystyle }\underset{a}{\overset{b}{\int }}{Q}_{\mathrm{L}}\left(r_{\mathrm{0}}\right)\mathrm{d}{r}_{\mathrm{0}}\approx {\displaystyle \frac{b-a}{\mathrm{2}}}\sum _{i=\mathrm{1}}^{\mathrm{3}}{\mathit{\omega }}_i{Q}_{\mathrm{L}}\left(r_{\mathrm{L}i}\right).\end{array}$$

Here, a and b are the lower and upper limits of r₀, which are set to 2 and 500 µm based on Andreas (1990), and ω_i is the corresponding weight (${\mathit{\omega }}_{\mathrm{1}}={\mathit{\omega }}_{\mathrm{3}}=\mathrm{0.556}$, ω₂=0.889), obtained from McClarren (2018). The new fast algorithm for approximations of H_S,SP and H_L,SP is referred to as SPRAY-GQ hereafter.

Figure 2 The distribution of occurrence frequency in percentage for GQ radius nodes: (a) the first node of latent heat flux, (b) the second node of latent heat flux, (c) the third node of latent heat flux, (d) the first node of sensible heat flux, (e) the second node of sensible heat flux, (f) and the third node of sensible heat flux. The peak frequencies are marked.

3.2 CFSv2.0-WW3 coupling system

A coupled system based on Climate Forecast System model version 2.0 (CFSv2.0) and WAVEWATCH III (WW3) is employed to evaluate and compare the effects of sea-spray-mediated heat flux parameterized by A15 and SPRAY-GQ. The CFSv2.0-WW3 has three components, the Global Forecast System (GFS; http://www.emc.ncep.noaa.gov/GFS/doc.php, last access: 18 March 2023) as the atmosphere component of CFSv2.0, the Modular Ocean Model version 4 (MOM4; Griffies et al., 2004) as the ocean component of CFSv2.0, and the WW3 (WAVEWATCH III Development Group, 2016) as the ocean surface wave component. The variables between CFSv2.0 and WW3 are interpolated and passed using the Chinese Community Coupler version 2.0 (C-Coupler2; Liu et al., 2018).

The CFSv2.0 is mainly applied for intraseasonal and seasonal prediction (e.g., Saha et al., 2014). The atmosphere component GFS uses a spectral triangular truncation of 382 waves (T382) in the horizontal, equivalent to a grid resolution of nearly 35 km, and 64 sigma–pressure hybrid layers in the vertical. The MOM4 is integrated on a nominal 0.5^∘ horizontal grid with enhanced horizontal resolution to 0.25^∘ in the tropics, and there are 40 levels in the vertical. The CFSv2.0 initial fields at 00:00 UTC of the first day for experiments were generated by the real-time operational Climate Data Assimilation System (Kalnay et al., 1996), downloaded from the CFSv2.0 official website (http://nomads.ncep.noaa.gov/pub/data/nccf/com/cfs/prod, last access: 18 March 2023) The latitude range of WW3 is 78^∘ S–78^∘ N with a spatial resolution of $\mathrm{1}/\mathrm{3}$^∘. The initial wave fields were generated from 10 d simulations starting from rest in a standalone WW3 model, forced by ERA5 10 m winds and ice concentration. The open boundary conditions of WW3 were also obtained by the global simulation of the standalone WW3 model.

In the coupling system, the WW3 obtains 10 m wind and ocean surface current from CFSv2.0 and then provides wave parameters to CFSv2.0. Several wave-mediated processes, including upper ocean mixing modified by Stokes drift-related processes, air–sea fluxes modified by surface current and Stokes drift, and momentum roughness length, are considered. Details of this system are referred to in Shi et al. (2022).

A series of numerical experiments is conducted to evaluate the effects of the two fast algorithms (A15 and SPRAY-GQ) of sea-spray-mediated heat flux on ocean, atmosphere and waves in two 56 d periods, from 3 January to 28 February 2017 and from 3 August to 28 September 2018 for boreal winter and boreal summer, respectively. For each period, two sensitivity experiments are carried out. The first is the SPRAY-A15 experiment in which A15 is used with two-way full coupling. The second is the SPRAY-GQ experiment in which SPRAY-GQ fast algorithm is used instead of A15. In addition, we also carry out another 7 d experiment using A92 (SPRAY-A92) to test the run time.

4 Results

4.1 Comparison with A92

Based on the daily global WSP10, T02, 2 m dew point temperature, surface pressure and SWH of ERA5, the daily global OISST, and the ESA monthly global salinity, H_S,SP and H_L,SP from A15, SPRAY-GQ and A92 are calculated (Fig. 3). The computational time for SPRAY-GQ is about the same as that for A15 and about 36 times less than the time for A92. Compared with A92 (the dotted black line), A15 (red) overestimates H_S,SP for low H_S,SP (< 50 W m⁻²) and underestimates H_S,SP for high H_S,SP (> 50 W m⁻²) with a root mean square error (RMSE $=\sqrt{\sum _{i=\mathrm{1}}^n(\widehat{y_i}-y_i{)}^{\mathrm{2}}/n}$, $\widehat{y_i}$ is A15 value, y_i is A92 value and n is the total number of grid points) of 3.40 W m⁻² (Fig. 3a), while A15 shows consistent overestimations with a RMSE of 2.98 W m⁻² for H_L,SP (Fig. 3b). Overall, the RMSE of A15 is about 2.69 W m⁻² for sea-spray-mediated total heat flux (${\mathrm{TH}}_{\mathrm{SP}}=H_{\mathrm{S},\mathrm{SP}}+H_{\mathrm{L},\mathrm{SP}}$; Fig. 3c). Andreas et al. (2015) derived A15 from A92 using single-radius droplets as bellwethers and wind functions, and extrapolated the wind functions at high wind speeds > 25 m s⁻¹. Since the wind speeds in the study are less than 25 m s⁻¹ (Fig. S1), the large difference between A15 and A92 is mainly due to the use of single-radius droplets. Compared with A15, SPRAY-GQ (blue) has less deviation from A92 for both H_S,SP and H_L,SP (Fig. 3a and b). The corresponding RMSEs of SPRAY-GQ for H_S,SP, H_L,SP and TH_SP are 0.83, 0.92 and 0.62 W m⁻², all significantly lower (P<0.05 in Student's t test) than those of A15.

Figure 3 Scatterplots of H_S,SP (a), H_L,SP (b) and total heat flux ${\mathrm{TH}}_{\mathrm{SP}}=H_{\mathrm{S},\mathrm{SP}}+H_{\mathrm{L},\mathrm{SP}}$ (c) estimated by fast algorithms (y axis) vs. those estimated by spectral integral in microphysical parameterization (x axis). The dotted black line is y=x. The corresponding RMSEs are marked in the upper left corner.

To test the robustness of the results, we also use WSP10, T02 and SPH of the OAFlux dataset to estimate H_S,SP and H_L,SP. As shown in Fig. 4, SPRAY-GQ has significantly (P<0.05 in Student's t test) lower deviations and RMSEs than A15, consistent with Fig. 3. Note that the values of H_S,SP and H_L,SP in Fig. 4 are larger than those in Fig. 3. This is because OAFlux only provides neutral wind speeds, calculated from wind stress and the corresponding roughness by assuming air is neutrally stratified. The neutral winds from OAFlux are larger than winds in ERA5, as indicated by previous studies (Lindemann et al., 2021; Seethala et al., 2021).

Figure 4 The same as Fig. 3, but WSP10, 2 m air temperature and 2 m specific humidity of OAFlux are used.

In addition, since it is common to derive SWH from empirical equations (e.g., Andreas et al., 2008, 2015; Andreas and Decosmo, 2002; Andreas, 1992), we also use SWH generated by empirical equations of WSP10 (Andreas, 1992) instead of ERA5 SWH to estimate H_S,SP and H_L,SP (Fig. 5). Again, the RMSEs decrease significantly (P<0.05 in Student's t test) in SPRAY-GQ compared to A15, though the RMSEs become higher for all estimates due to the enhanced biases of SWH. The difference between SPRAY-GQ and A92 is always smaller than that between A15 and A92. Next, we will evaluate and compare the two fast algorithms in an atmosphere–ocean–wave coupled system (CFSv2.0-WW3).

Figure 5 The same as Fig. 4, but SWH is derived by WSP10 instead of ERA5 SWH.

4.2 Comparison in the CFSv2.0-WW3 coupling system

To compare the computational time of different parameterizations in the large-scale modeling system, the run time of the fully coupled experiments for 7 d forecast is given in Table 1 as an example. It is shown that the run time is about the same for SPRAY-GQ and SPRAY-A15. Both experiments run about 17 times faster than SPRAY-A92.

To illustrate the numerical errors of the two fast algorithms discussed in the context of the coupled system, comparisons are made for simulated SSTs, WSP10s and SWHs against OISST and ERA5 reanalysis. The results in the first 3 d are excluded in the comparison, since the wave influences are weak at the beginning of the simulations. Overall, the WSP10s of simulations are generally in the range of 0–25 m s⁻¹ globally. At middle to high latitudes, the WSP10s generally exceed 10 m s⁻¹ (Figs. S2 and S3 of the Supplement), at which the effects of sea spray can become significant (Andreas et al., 2015, 2008).

Figure 6 The 53 d average SST (^∘C) differences between SPRAY-A15 and OISST (a; SPRAY-A15 minus OISST), the differences between SPRAY-GQ and SPRAY-A15 (b; SPRAY-GQ minus SPRAY-A15), and the time series of domain-averaged RMSE and MAE (c; 0–360^∘ E, 40–75^∘ S) in January–February 2017. The first 3 d simulation is discarded. The dotted areas are statistically significant at 95 % confidence level.

Figure 7 The 53 d average differences of total heat flux (a–c), latent heat flux (d–f), and sensible heat flux (g–i) between SPRAY-GQ and SPRAY-A15 (SPRAY-GQ minus SPRAY-A15) in January–February 2017. The direct differences indicate sea-spray-mediated heat flux differences (b, e, h), and the indirect differences indicate interfacial (bulk) heat flux differences resulting from sea spray (c, f, i). The dotted areas are statistically significant at 95 % confidence level. A positive value of flux indicates an upward direction.

4.2.1 Sea surface temperature (SST)

In the austral summer, compared with OISST, large SST biases (> 1 ^∘C or $<-\mathrm{1}$ ^∘C) of SPRAY-A15 occur in the Southern Hemisphere (SH; Fig. S4a in the Supplement), especially in the Southern Ocean. It is always a challenge to reduce the large SST biases in the Southern Ocean for climate models (e.g., Alessandro et al., 2019; Wang et al., 2014; Li et al., 2013; Bodas-Salcedo et al., 2012; Ceppi et al., 2012). In Fig. 6a, SSTs north (south) of 50^∘ S in experiment SPRAY-A15 are mainly underestimated (overestimated). The domain-averaged RMSE (0–360^∘ E, 40–75^∘ S) in experiment SPRAY-A15 increases in the first month and then levels off (solid red line in Fig. 6c), while the domain-averaged RMSE in experiment SPRAY-GQ levels off about a week earlier (solid black line in Fig. 6c). The mean RMSE in SPRAY-GQ is significantly lower than that in SPRAY-A15 (P<0.05 in Student's t test). The increased (decreased) SSTs north (south) of 50^∘ S in SPRAY-GQ compared to those in SPRAY-A15 (Fig. 6b) reduce the RMSE of SST in SPRAY-GQ. We also calculate the mean absolute error, MAE $=\sum _{i=\mathrm{1}}^n\mathrm{|}\widehat{y_i}-y_i\mathrm{|}/n$, where $\widehat{y_i}$ is simulated value, y_i is OISST data and n is the total number of grid points. The MAEs are consistent with RMSEs (dotted line in Fig. 6c). Furthermore, the mean errors, ME $=\sum _{i=\mathrm{1}}^n(\widehat{y_i}-y_i)/n$ (Fig. S5a in the Supplement), are smaller in SPRAY-GQ than SPRAY-A15.

To understand the effects of sea spray droplets on SST, we calculate the total heat flux (TH $=H_{\mathrm{S},\mathrm{T}}+H_{\mathrm{L},\mathrm{T}})$ differences between SPRAY-GQ and SPRAY-A15 (Fig. 7a). The TH differences are significantly correlated with SST differences (Fig. S4b in the Supplement), with the spatial correlation coefficient of −0.41 (P<0.05 in Student's t test). We further decompose direct and indirect effects of sea spray droplets on heat fluxes following Song et al. (2022). The direct effect (H_S,SP and H_L,SP) is induced directly by sea spray droplets, calculated from A15 (Eqs. B1–B4 of Appendix B) and SPRAY-GQ (Sect. 3.1). The indirect effect (H_S and H_L) is the heat flux variation induced by changes of atmosphere and ocean variables (including wind, pressure, humidity and temperature) caused by direct effect, estimated by subtracting H_S,SP and H_L,SP from the output heat fluxes (H_S,T and H_L,T) of experiment SPRAY-A15 and SPRAY-GQ.

Figure 8 The same as Fig. 6, but for August–September 2018 in 0–360^∘ E, 20–75^∘ N.

Figure 9 The same as Fig. 7, but for August–September 2018.

In the Southern Ocean, although direct differences of H_S,SP and H_L,SP are relatively small (< 10 W m⁻², Fig. 7b, e and h), the resulting changes of temperature and humidity lead to relatively large differences in indirect effects of H_S and H_L (Fig. 7c, f and i). Enhanced (reduced) TH_SP from ocean to atmosphere in the summer leads to increased (decreased) air–sea temperature difference and thus enhances (weakens) H_S. Meanwhile the warmer (cooler) air also causes more (less) evaporation and thus more (less) H_L. Finally, the enhanced (reduced) TH cools (warms) SST.

Figure 10 The 53 d average WSP10 (m s⁻¹) differences between SPRAY-A15 and ERA5 (a; SPRAY-A15 minus ERA5), the differences between SPRAY-GQ and SPRAY-A15 (b; SPRAY-GQ minus SPRAY-A15), and the time series of domain-averaged RMSE and MAE (c; 0–360^∘ E, 0–60^∘ N) in January–February 2017. The first 3 d simulation is discarded. The dotted areas are statistically significant at 95 % confidence level.

Figure 11 The same as Fig. 10, but for August–September 2018.

In the boreal summer, large SST biases (> 1 ^∘C or $<-\mathrm{1}$ ^∘C) of SPRAY-A15 mainly occur at middle to high latitudes of the Northern Hemisphere (NH; Fig. S6a in Supplement). Significant underestimations occur in the western and northern part of the North Pacific and at mid-latitudes of the North Atlantic, while large positive SST biases mainly occur in the eastern part of the North Pacific and at high latitudes of the North Atlantic (Fig. 8a). In experiment SPRAY-GQ, SSTs are warmer (cooler) in the previously underestimated (overestimated) regions (Fig. 8b). Therefore, the domain-averaged RMSE and MAE (0–360^∘ E, 20–75^∘ N) in SPRAY-GQ are significantly lower (P<0.01 in Student's t test) than in SPRAY-A15 after the first 3 weeks (Fig. 8c). Compared to SPRAY-A15, the overall underestimation is reduced in SPRAY-GQ (Fig. S5b). The spatial correlation coefficient between TH differences and SST differences (Figs. 9a, S6b) is −0.32 (P<0.05 in Student's t test). Consistent with the austral summer, the SST changes are related to the changes of heat flux (Fig. 9). The indirect effects of latent heat flux (Fig. 9f) play a major role in TH differences, which are modified by the direct effects (Fig. 9b, e and h). In addition, the changes of surface winds also contribute to the changes of SST. The reduced winds weaken the upper ocean mixing, the water becomes more stratified, and then the SST tends to be warmer and vice versa (Figs. S7 and S8).

4.2.2 Significant wave height (SWH) and 10 m wind speed (WSP10)

Compared with experiment SPRAY-A15, significant differences of WSP10 in SPRAY-GQ occur at middle to low latitudes of the NH (0–360^∘ E, 0–60^∘ N) in both winter and summer (Figs. S7b and S8b). As we know, satellite scatterometer and altimeter data are usually used to validate WSP10 and SWH for short-term weather forecast (e.g., Accadia et al., 2007; Djurdjevic and Rajkovic, 2008; Myslenkov et al., 2021). However, due to the spatial and temporal coverage of satellite data, we can only obtain the monthly averaged satellite data for the globe. Therefore, we compare the monthly averaged WSP10 and SWH from simulations with the corresponding satellite data (Figs. S9–S12). The comparison results (Figs. S9a and c–S12a and c) are consistent with those compared with ERA5 (Figs. S9b and d–S12b and d). From Figs. S9e–S12e, the differences of WSP10s between ERA5 and the satellite data are always less than 1 m s⁻¹, and the differences of SWHs are always less than 0.3 m. Since ERA5 provides daily data for comparison, we will use ERA5 for validation in the following.

The ME of WSP10 (SPRAY-A15 minus ERA5) is 0.28 and 0.47 m s⁻¹ in winter and summer (red in Fig. S5c and d), respectively, mainly due to the overestimations over the Pacific and the Atlantic Ocean (red in Figs. 10a and 11a), whereas in SPRAY-GQ the ME (SPRAY-GQ minus ERA5) is 0.15 and 0.33 m s⁻¹ in winter and summer, respectively (black in Fig. S5c and d). The domain-averaged RMSEs and MAEs of WSP10s increase with time in the first 2 weeks and then gradually level off (Figs. 10c and 11c). The differences of WSP10 RMSEs and MAEs between SPRAY-GQ (black) and SPRAY-A15 (red) are very small in the first 2 weeks. Afterwards, the mean values of RMSE and MAE in SPRAY-GQ are significantly lower than those in SPRAY-A15, at 95 % confidence level in both boreal winter (Fig. 10c) and boreal summer (Fig. 11c).

The simulated SWHs changes are closely related to the changes of WSP10s (Shi et al., 2022). Therefore, the differences of SWHs (Figs. 12 and 13) are consistent with those of WSP10s (Figs. 10 and 11), with overestimated (underestimated) WSP10s corresponding to overestimated (underestimated) SWHs compared with ERA5. The SWHs in SPRAY-GQ are significantly different from those in SPRAY-A15 (Figs. 12b and 13b). In winter (summer), the SWH RMSE averages for SPRAY-A15 and SPRAY-GQ are 1.31 m (0.98 m) and 1.23 m (0.87 m), and after the first 2 weeks the RMSE and MAE in SPRAY-GQ are significantly lower than those in SPRAY-A15 at 95 % confidence level in both winter (Fig. 12c) and summer (Fig. 13c).

Figure 12 The 53 d average SWH (m) differences between SPRAY-A15 and ERA5 (a; SPRAY-A15 minus ERA5), the differences between SPRAY-GQ and SPRAY-A15 (b; SPRAY-GQ minus SPRAY-A15), and the time series of domain-averaged RMSE and MAE (c; 0–360^∘ E, 0–60^∘ N) in January–February 2017. The first 3 d simulation is discarded. The dotted areas are statistically significant at 95 % confidence level.

Figure 13 The same as Fig. 12, but for August–September 2018.

The direct and indirect effects of sea spray droplets on heat fluxes can influence estimates of WSP10 and then SWH. The changes of WSP10s are related to the direct effects (H_S,SP and H_L,SP; Fig. 7b, e and h; Fig. 9b, e and h). The spatial correlation coefficients between WSP10 differences (Figs. S7b and S8b) and TH_SP differences (Figs. 7b and 9b) are 0.51 and 0.69 (P<0.01 in Student's t test) in winter and summer, respectively, because TH_SP differences can influence the sea level pressure (SLP) distribution (Figs. S15 and S16) and subsequently surface winds. For example, compared with SPRAY-A15, the decreased TH_SP of SPRAY-GQ in the Northwest Pacific in summer (Fig. 9b) leads to higher SLP and a smaller pressure gradient (Fig. S16) and thus decreased WSP10 (Fig. 11b), while the increased TH_SP in the Gulf of Alaska (Fig. 9b) leads to lower SLP and larger pressure gradient (Fig. S16) and thus enhanced WSP10 (Fig. 11b). The accelerated (decelerated) WSP10s further result in increased (decreased) interfacial heat transport (H_S, H_L) and increased (decreased) SWHs.

5 Conclusions and discussion

Based on a GQ method, we develop a new fast algorithm based on Andreas's (1989, 1990, 1992) full-size microphysical parameterization (A92) for sea-spray-mediated heat fluxes. Using global satellite measurements and reanalysis data, we found that the difference between SPRAY-GQ and A92 is significantly smaller than that between A15 and A92 (Andreas et al., 2015). To evaluate the numerical error of the SPRAY-GQ/A15 fast algorithm, we implement them in the two-way coupled CFSv2.0-WW3 system. A series of 56 d simulations from 3 January to 28 February 2017 and from 3 August to 28 September 2018 are conducted. The results are compared against satellite measurements and ERA5 reanalysis. The comparison shows that the sea-spray-mediated heat flux in SPRAY-GQ can reasonably modulate total heat flux compared with SPRAY-A15 and significantly reduce the SST biases in the Southern Ocean (middle to high latitudes of the NH) for the austral (boreal) summer and WSP10 and SWH after the first 2 weeks at middle to low latitudes of the NH for both boreal winter and summer. Overall, our fast algorithm based on GQ is applicable to sea-spray-mediated heat flux parameterization in coupled models.

To investigate the effects of spray-mediated heat flux on simulations, two 56 d experiments without sea spray effect (CTRL) in boreal winter and summer are conducted, respectively, and the differences of simulated SST, WSP10, SWH, T02 and SPH between SPRAY-GQ and CTRL are compared in Figs. S17–S21 in the Supplement. The introduction of sea spray cannot significantly reduce the global overall errors of simulations, but it leads to regional improvements (blue in Figs. S17e and f–S21e and f). For example, compared with CTRL in January–February 2017, SST MAE of SPRAY-GQ in the southeast of Australia decreases (Fig. S17e) because of warmer SST (Fig. S17c) related to reduced wind (Fig. S18c). The reduced wind here also leads to lower SWH (Fig. S19c) and thus reduced SWH overestimation (Fig. S19e). Meanwhile, SPRAY-GQ reduces MAE of T02 and SPH (Figs. S20e and S21e) by increasing temperature and moisture (Figs. S20c and S21c). The reduced errors are related to the relatively large WSP10s over the areas (Figs. S2 and S3), since the effects of sea spray become important at wind speeds larger than 10 m s⁻¹.

In addition to the variables aforementioned, the changes of simulated cloud fraction were also compared. However, the effects of sea-spray-mediated heat flux on cloud fraction are non-significant for the 2-month simulation, so the results are not shown. Besides this, the lack of other processes related to sea spray may be one of the reasons why the global overall error cannot be reduced effectively. For example, for simulated WSP10 and SWH in SPRAY-GQ, the significant overestimations in the SH still exist especially in August–September 2018 (Figs. S18 and S19 in the Supplement). As Andreas (2004) indicated, sea spray droplets also influence the surface momentum flux by injecting more momentum into the ocean from the atmosphere, which might further decrease the surface wind speed. We will consider this process in a future study.

Sea-spray-mediated heat fluxes are related to the sea spray generation function (SSGF). Based on a number of laboratory and field observations, varieties of SSGF were derived (e.g., Koga, 1981; Monahan et al., 1982; Troitskaya et al., 2018; Andreas, 1992, 1998, 2002; Fairall et al., 1994; Veron, 2015), whereas their differences can reach six orders of magnitude (Andreas, 1998). There is currently no consensus on the most suitable choice. In this study, we use SSGF of Fairall et al. (1994), recommended by Andreas (2002), to get a mean bias of 3.70 and 0.095 W m⁻² for latent and sensible heat flux, respectively (Andreas et al., 2015), consistent with recent observations of Xu et al. (2021b). However, the improved SST and other variables cannot be reliably assigned to the usage of the GQ method, due to the uncertainties of the coupled model itself and SSGF.

When wind speed is larger than 10 m s⁻¹, spray-mediated heat flux can become as important as the interfacial heat flux (Andreas and Decosmo, 1999, 2002). Particularly, even in the absence of air–sea temperature difference, the spray-mediated sensible heat flux is still present (Andreas et al., 2008). As indicated by previous studies (e.g., Garg et al., 2018; Song et al. 2022), it is necessary to superimpose the spray-mediated heat flux on the bulk formula to complete the physics of turbulent heat transfer for coupled simulation. Since the full microphysical parameterization (A92) is computationally expensive, an efficient algorithm that captures the main features of A92 can be beneficial to large-scale climate systems or operational storm models. The GQ method proposed in the study can efficiently calculate the spray-mediated heat flux and agrees better with A92 than A15. Therefore, the GQ-based spray-mediated heat flux is promising for wide applications in large-scale climate systems and operational storm models.

Appendix A: Microphysical parameterization of A92

Based on the cloud microphysical parameterization of Pruppacher and Klett (1978), Andreas (1989, 1990, 1992) proposed a parameterization of sea-spray-related heat fluxes for droplets with different radius, from formation at sea surface to equilibrium with environment; that is,

$$\begin{array}{}\text{(A1)}& {\displaystyle }{Q}_{\mathrm{S}}={\mathit{\rho }}_{\mathrm{w}}{C}_{\mathrm{ps}}\left(T_{\mathrm{w}}-T_{\mathrm{eq}}\right)\left[\mathrm{1}-\exp \left(-{\displaystyle \frac{{\mathit{\tau }}_f}{{\mathit{\tau }}_T}}\right)\right]\left({\displaystyle \frac{\mathrm{4}\mathit{\pi }{r}_{\mathrm{0}}^{\mathrm{3}}}{\mathrm{3}}}{\displaystyle \frac{\mathrm{d}F}{\mathrm{d}{r}_{\mathrm{0}}}}\right)\text{(A2)}& {\displaystyle }{Q}_{\mathrm{L}}=\left\{\begin{array}{l}{\mathit{\rho }}_{\mathrm{w}}{L}_{\mathrm{v}}\left\{\mathrm{1}-{\left[\frac{r\left({\mathit{\tau }}_{\mathrm{f}}\right)}{r_{\mathrm{0}}}\right]}^{\mathrm{3}}\right\}\left(\frac{\mathrm{4}\mathit{\pi }{r}_{\mathrm{0}}^{\mathrm{3}}}{\mathrm{3}}\frac{\mathrm{d}F}{\mathrm{d}{r}_{\mathrm{0}}}\right),{\mathit{\tau }}_{\mathrm{f}}\le {\mathit{\tau }}_{\mathrm{r}},\\ {\mathit{\rho }}_{\mathrm{w}}{L}_{\mathrm{v}}\left\{\mathrm{1}-{\left(\frac{r_{\mathrm{eq}}}{r_{\mathrm{0}}}\right)}^{\mathrm{3}}\right\}\left(\frac{\mathrm{4}\mathit{\pi }{r}_{\mathrm{0}}^{\mathrm{3}}}{\mathrm{3}}\frac{\mathrm{d}F}{\mathrm{d}{r}_{\mathrm{0}}}\right),{\mathit{\tau }}_{\mathrm{f}}>{\mathit{\tau }}_{\mathrm{r}}.\end{array}\right.\end{array}$$

Here, Q_S and Q_L are sensible heat flux and latent heat flux resulting from sea spray droplets with initial radius r₀, ρ_w is the sea water density, C_ps is the specific heat, L_v is the latent heat of vaporization of water, T_w is the water temperature, T_eq is the temperature of droplet when it reaches thermal equilibrium with ambient condition, r_eq is the radius of droplet when it reaches moisture equilibrium with ambient condition, τ_f is the residence time for droplets in the atmosphere, r(τ_f) is the corresponding radius, τ_T is the characteristic e folding time of droplet temperature, and τ_r is the characteristic e folding time of droplet radius. The detailed calculation of these microphysical quantities can be found in Andreas (1989, 1990, 1992). $\mathrm{d}F/\mathrm{d}{r}_{\mathrm{0}}$ is the sea spray generation function, which represents the produced number of droplets with initial radius r₀ (Andreas, 1992). For this term, the function of Fairall et al. (1994) was recommended by Andreas (2002). According to the review in Andreas (2002), the $\mathrm{d}F/\mathrm{d}{r}_{\mathrm{0}}$ of Fairall et al. (1994) is related on that of Andreas (1992) as

$$\begin{array}{}\text{(A3)}& {\displaystyle \frac{\mathrm{d}F}{\mathrm{d}{r}_{\mathrm{0}}}}=\mathrm{38}\times \mathrm{3.84}\times {\mathrm{10}}^{-\mathrm{6}}{U}_{\mathrm{10}}^{\mathrm{3.41}}{r}_{\mathrm{0}}^{-\mathrm{0.024}}{\left.{\displaystyle \frac{\mathrm{d}{F}_{A\mathrm{92}}}{\mathrm{d}{r}_{\mathrm{80}}}}\right|}_{U_{\mathrm{10}}=\mathrm{11}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}},\text{(A4)}& {\displaystyle }\begin{array}{rl}& {\left.{\displaystyle \frac{\mathrm{d}{F}_{A\mathrm{92}}}{\mathrm{d}{r}_{\mathrm{80}}}}\right|}_{U_{\mathrm{10}}=\mathrm{11}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}}=\\ & \left\{\begin{array}{l}{e}^{\left(\mathrm{4.405}-\mathrm{2.646}\left(\log r_{\mathrm{80}}\right)-\mathrm{3.156}{\left(\log r_{\mathrm{80}}\right)}^{\mathrm{2}}+\mathrm{8.902}{\left(\log r_{\mathrm{80}}\right)}^{\mathrm{3}}-\mathrm{4.482}{\left(\log r_{\mathrm{80}}\right)}^{\mathrm{4}}\right)},\\ r_{\mathrm{80}}\le \mathrm{15}\mathrm{\mu }\mathrm{m};\\ \mathrm{1.02}\times {\mathrm{10}}^{\mathrm{4}}{r}_{\mathrm{80}}^{-\mathrm{1}},\mathrm{15}\le r_{\mathrm{80}}\le \mathrm{37.5}\mathrm{\mu }\mathrm{m};\\ \mathrm{6.95}\times {\mathrm{10}}^{\mathrm{6}}{r}_{\mathrm{80}}^{-\mathrm{2.8}},\mathrm{37.5}\le r_{\mathrm{80}}\le \mathrm{100}\mathrm{\mu }\mathrm{m};\\ \mathrm{1.75}\times {\mathrm{10}}^{\mathrm{17}}{r}_{\mathrm{80}}^{-\mathrm{8}},r_{\mathrm{80}}\ge \mathrm{100}\mathrm{\mu }\mathrm{m}.\end{array}\right.\end{array}\end{array}$$

Here, U₁₀ is the 10 m wind and $r_{\mathrm{80}}=\mathrm{0.518}{r}_{\mathrm{0}}^{\mathrm{0.976}}$.

The total sea spray fluxes are obtained by integrating Q_S and Q_L corresponding to all r₀. Based on Andreas (1990), the lower and upper limits of r₀ are 2 and 500 µm, respectively; that is,

$$\begin{array}{}\text{(A5)}& {\displaystyle }\stackrel{\mathrm{‾}}{Q_{\mathrm{S}}}=\underset{\mathrm{2}}{\overset{\mathrm{500}}{\int }}{Q}_{\mathrm{S}}\left(r_{\mathrm{0}}\right)\mathrm{d}r,\text{(A6)}& {\displaystyle }\stackrel{\mathrm{‾}}{Q_{\mathrm{L}}}=\underset{\mathrm{2}}{\overset{\mathrm{500}}{\int }}{Q}_{\mathrm{L}}\left(r_{\mathrm{0}}\right)\mathrm{d}r.\end{array}$$

Note that $\stackrel{\mathrm{‾}}{Q_{\mathrm{S}}}$ and $\stackrel{\mathrm{‾}}{Q_{\mathrm{L}}}$ are nominal sea spray fluxes but not the actual H_S,SP and H_L,SP (Andreas and Decosmo, 1999, 2002), because there are interactions between these two terms and the microphysical functions also lead to uncertainties (Fairall et al., 1994). Therefore, $\stackrel{\mathrm{‾}}{Q_{\mathrm{S}}}$ and $\stackrel{\mathrm{‾}}{Q_{\mathrm{L}}}$ are tuned by non-negative constants α, β and γ (Andreas and Decosmo, 2002; Andreas et al., 2008, 2015; Andreas, 2003) as

$$\begin{array}{}\text{(A7)}& {\displaystyle }{H}_{\mathrm{S},\mathrm{SP}}=\mathit{\beta }\stackrel{\mathrm{‾}}{Q_{\mathrm{S}}}-(\mathit{\alpha }-\mathit{\gamma })\stackrel{\mathrm{‾}}{Q_{\mathrm{L}}},\text{(A8)}& {\displaystyle }{H}_{\mathrm{L},\mathrm{SP}}=\mathit{\alpha }\stackrel{\mathrm{‾}}{Q_{\mathrm{L}}}\end{array}$$

In Eq. (A8), the α term indicates the sea-spray-mediated latent heat flux from the top of DEL to atmosphere. Because the evaporation of droplets absorbs heat, which is provided by sea-spray-mediated sensible heat (Fairall et al., 1994), the negative α term appears in Eq. (A7). The evaporation also cools DEL and thus increases the air–sea temperature difference: therefore it contributes to a positive γ term in Eq. (A7). Different values of α, β and γ were given in Andreas and Decosmo (2002), Andreas (2003), and Andreas et al. (2008, 2015) to minimize the bias between estimations and observations of turbulent heat fluxes measured by eddy correlation. And Andreas et al. (2015) validated the most observation data, which are 4000 sets, to derive $\mathit{\alpha }=\mathrm{2.46},\mathit{\beta }=\mathrm{15.15},\mathit{\gamma }=\mathrm{1.77}$.

Appendix B: Fast algorithm of A15

Andreas (2003) and Andreas et al. (2008, 2015) developed a fast algorithm to approximate H_S,SP and H_L,SP by a characteristic radius; that is,

$$\begin{array}{}\text{(B1)}& {\displaystyle }\begin{array}{rl}{H}_{\mathrm{S},\mathrm{SP}}& =\mathit{\beta }\stackrel{\mathrm{‾}}{Q_{\mathrm{S}}}-(\mathit{\alpha }-\mathit{\gamma })\stackrel{\mathrm{‾}}{Q_{\mathrm{L}}}\\ & \approx {\mathit{\rho }}_{\mathrm{w}}{C}_{\mathrm{ps}}\left(T_{\mathrm{w}}-T_{\mathrm{eq},\mathrm{100}}\right)V_{\mathrm{s}}\left(u_{\ast }\right),\end{array}\text{(B2)}& {\displaystyle }{H}_{\mathrm{L},\mathrm{SP}}=\mathit{\alpha }\stackrel{\mathrm{‾}}{Q_{\mathrm{L}}}\approx {\mathit{\rho }}_{\mathrm{w}}{L}_{\mathrm{v}}\left\{\mathrm{1}-{\left[{\displaystyle \frac{r\left({\mathit{\tau }}_{\mathrm{f},\mathrm{50}}\right)}{\mathrm{50}\mathrm{\mu }\mathrm{m}}}\right]}^{\mathrm{3}}\right\}{V}_{\mathrm{L}}\left(u_{\ast }\right).\end{array}$$

Here, T_eq,100 is T_eq of droplets with r₀=100 µm, τ_f,50 is τ_f of droplets with r₀=50 µm, and V_s and V_L are functions of the bulk friction velocity u_∗. As indicated by Andreas et al. (2008, 2015), the characteristic radiuses of 100 µm and 50 µm for sensible and latent heat fluxes are chosen, respectively, because Q_S and Q_L show a large peak in the vicinity of these values (Fig. 1). V_s and V_L are calculated in Andreas et al. (2015) as

$$\begin{array}{}\text{(B3)}& {\displaystyle }{V}_{\mathrm{S}}=\left\{\begin{array}{l}\mathrm{3.92}\times {\mathrm{10}}^{-\mathrm{8}},\mathrm{0}\le u_{\ast }\le \mathrm{0.1480}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}\\ \mathrm{5.02}\times {\mathrm{10}}^{-\mathrm{6}}{u}_{\ast }^{\mathrm{2.54}},u_{\ast }\ge \mathrm{0.1480}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}\end{array}\right.,\text{(B4)}& {\displaystyle }{V}_{\mathrm{L}}=\left\{\begin{array}{l}\mathrm{1.76}\times {\mathrm{10}}^{-\mathrm{9}},\mathrm{0}\le u_{\ast }\le \mathrm{0.1358}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}\\ \mathrm{2.08}\times {\mathrm{10}}^{-\mathrm{7}}{u}_{\ast }^{\mathrm{2.39}},u_{\ast }\ge \mathrm{0.1358}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}\end{array}\right..\end{array}$$

Appendix C: Gaussian quadrature (GQ)

GQ is a method to approximate the definite integral of a function f(x) via the function values at a small number of specified nodes (Gauss, 1815; Jacobi, 1826). In this study we use the form of n-node Gauss–Legendre quadrature on $[-\mathrm{1},\mathrm{1}]$ as

$$\begin{array}{}\text{(C1)}& \underset{-\mathrm{1}}{\overset{\mathrm{1}}{\int }}f\left(x\right)\mathrm{d}x\approx \sum _{i=\mathrm{1}}^n{\mathit{\omega }}_{i}f\left(x_i\right).\end{array}$$

Here, x_i is the specified node, and ω_i is the corresponding weight. For n=3, $x_{\mathrm{1}}=-\mathrm{0.775}$, x₂=0, x₃=0.775, ${\mathit{\omega }}_{\mathrm{1}}={\mathit{\omega }}_{\mathrm{3}}=\mathrm{0.556}$, ω₂=0.889.

For a function g(ξ) on [a, b], Eq. (C1) can be transformed to

$$\begin{array}{}\text{(C2)}& \begin{array}{rl}\underset{a}{\overset{b}{\int }}g\left(\mathit{\xi }\right)\mathrm{d}\mathit{\xi }& =\underset{-\mathrm{1}}{\overset{\mathrm{1}}{\int }}g\left({\displaystyle \frac{b-a}{\mathrm{2}}}x+{\displaystyle \frac{a+b}{\mathrm{2}}}\right){\displaystyle \frac{\mathrm{d}\mathit{\xi }}{\mathrm{d}x}}\mathrm{d}x\\ & \approx {\displaystyle \frac{b-a}{\mathrm{2}}}\sum _{i=\mathrm{1}}^n{\mathit{\omega }}_{i}g\left({\displaystyle \frac{b-a}{\mathrm{2}}}{x}_i+{\displaystyle \frac{a+b}{\mathrm{2}}}\right).\end{array}\end{array}$$