We use the ERA5 hourly reanalysis dataset (Hersbach et al., 2020) to investigate the submonthly air-sea interaction in the western North Pacific around the KOE region. ERA5 is the fifth generation, the latest, of atmospheric reanalysis produced by the ECMWF. It is based on the IFS Cy41r2 model with a spatial resolution of ∼0.28° or ∼31 km (TL639) and 137 vertical levels up to 1 Pa. Data assimilation is performed using the 4D-Var method on the interval of 12 h. A large number of observations from various sources including satellite radiometers and scatterometers are assimilated. The high spatial resolution and the frequent observational adjustment ensure the usefulness of ERA5 in air-sea interaction studies on the submonthly timescales. The variables involved in this work include sea surface temperature (SST), air temperature at 2 m (TA), air-sea sensible and latent heat fluxes (SHF and LHF), net surface thermal and solar radiation (STR and SSR), total precipitation (PRE), net evaporation (EVP), zonal and meridional components of 10 m wind velocity (U₁₀ and V₁₀). Surface net heat flux is calculated as the sum of SHF, LHF, STR and SSR. Here positive heat flux indicates heat release from the ocean to the atmosphere. Net freshwater flux is defined as PRE-EVP. Surface wind stress is calculated following Yelland and Taylor (1996) as τ = ρ_a(a + bU)U², where ρ_a = 1.225 kg/m³ is air density, U=U102+V102, a = 6×10⁻⁴, b = 7×10⁻⁵. The use of this formula to calculate wind stress is because of a known issue of problematic wind stress provided by the ERA5 dataset.¹

For comparison, the 6-h NCEP/NCAR Reanalysis 1 dataset (Kistler et al., 2001) is also employed. The spatial resolution is 2.5°. Evaporation is converted from latent heat flux following EVP = LHF/ρ_wλ, where ρ_w = 1024 kg/m³ is sea water density, and λ = 2.5×10⁶−2386⋅SST is the latent heat of vaporization (Brutsaert, 1982). Wind stress is directly provided by the dataset.

The General Ocean Turbulence Model (GOTM; Burchard, 2002)² version 6.0.0 is employed in this work. It is a single-column model for the most important hydrodynamic and thermodynamic processes related to vertical mixing. The model is built upon a library of traditional and state-of-the-art turbulence closure parameterizations for vertical momentum, heat, and tracer fluxes (Umlauf and Burchard, 2005), which has been adopted in a number of global or regional 3D ocean models such as ROMS, FVCOM, FESOM, GETM and NEMO. It can also handle convective adjustment due to buoyancy disturbances. It is possible to consider horizonal advection by prescribing the velocity and temperature/salinity gradient, but this is not done in this work since we aim at simulating SST variability purely driven by the local atmosphere, adjusted by vertical mixing in the water column. The model is forced by surface heat, momentum, and freshwater fluxes from ERA5. Due to the ignorance of Ekman transport and the prescribed heat flux without modification by wind-induced anomalous advection of air temperature and humidity, the effect of wind is only in its mechanical stirring of the upper ocean. However, the prescribed observational heat flux should already have the thermohaline effects of winds implicitly incorporated. In this work, the model is run on 100 evenly distributed vertical layers from the sea surface to two times the climatological winter mixed layer depth at a certain geographical location, the mixed layer depth data taken from the ECCO2 ocean reanalysis (Menemenlis et al., 2008). Initial temperature and salinity conditions are also from ECCO2. The model is integrated from October 1 to February 28 of the specific winter under study (2001/02 and 2008/09, see below), with time interval of 1 h. Only model output of December-February is used.

The KOE region has the largest submonthly (< 30 days) SST variability in the North Pacific (Figure 1A), which makes up about 20% of the total SST standard deviation (STD; Figure 1B). Correspondingly, this region also exhibits significant submonthly TA and net heat flux (NHF) variability (Figures 1C,D), which is indicative of the vigorous high-frequency air-sea interaction there. From physical perspectives, submonthly atmospheric variability can be divided into 2 bands: 2–10 and 10–30 days, associated with weather events and low-frequency oscillations (Chang and Orlanski, 1993; Chang et al., 2002). In view of this, we examine submonthly SST variability and air-sea interaction in terms of these 2 bands, respectively. Figure 2A shows the multi-year time series of STD of 2–10 and 10–30 days bandpass filtered SST variability in boreal winter (DJF) averaged over the KOE region. It is obvious that 2–10 days SST variability is very weak before 2007, but it rises abruptly in winter 2007/08. Meanwhile, 10–30 days variability dominates before 2007, and drops dramatically afterward. This is also evident from the power spectra at 145°E, 38°N in two selected representative winters, 2001/02 and 2008/09 (Figure 2B), and the spatial distribution of standard deviation on the two scales in the two selected winters (Figures 2C–F).

The reason for the abrupt changes in SST variability is because ERA5 undergoes an update of SST boundary forcing in September 2007 from the 5-daily HadISST dataset to the daily OSTIA dataset (Hersbach et al., 2020). The temporal resolution before the update is too crude to resolve 2–10 days variability, but it becomes resolvable after the update, hence the standard deviation rises. The drop of 10–30 days standard deviation indicates that the variability on this scale before the update is too strong, probably subject to the aliasing artifact of the unresolvable higher frequencies. In the following we mainly focus on winter 2008/09 when investigating the characteristics of air-sea interaction on the two scales, and use winter 2001/02 as an opportunity to assess the potential role of different levels of SST variability on air-sea interaction. It is confirmed that selecting other representative years before and after the SST resolution shift does not change the main conclusions qualitatively. Moreover, the spatial resolution is changed from 0.25 to 0.05° along with the temporal resolution shift, but the 0.05° resolution oversamples the model grid resolution of ∼0.28°. The change of spatial resolution is, therefore, minimal.

Local, near-surface air-sea heat exchange can be described using the following theoretical model (Bishop et al., 2017):

where T_a and T_o are the air temperature and the SST, respectively. ε_a and ε_o are atmospheric and oceanic stochastic forcing terms conforming to the normal distribution. α, γ_a and ω_a are characteristic frequencies for thermal exchange, damping, and stochastic forcing of the atmosphere. Likewise, β, γ_o and ω_o are for the ocean. These 6 frequency parameters are all positive by definition, and their reciprocals are termed the corresponding timescales. Air-sea heat flux is represented by β(T_o−T_a). Incoming shortwave radiation and net thermal radiation are affected by cloud cover, not directly related to air-sea temperature difference, and are therefore not counted in the heat flux term. Latent heat flux changes SST but not surface air temperature and should therefore be included in (2) but not in (1), whereas sensible heat flux appears in both equations. This theoretical model is an extension to the classic stochastic climate model established by Hasselmann (1976) and Frankignoul and Hasselmann (1977) in which the SST variability is completely driven by atmospheric forcing (ω_o = 0).

Bishop et al. (2017) solved the above model (1 and 2) numerically and argued that in case of pure atmospheric forcing, simultaneous correlation between SST and heat flux is close to zero, but simultaneous correlation between SST tendency (dT₀/dt) and heat flux is negative and large. On the other hand, in case of pure oceanic forcing (ω_a = 0), SST-heat flux correlation is a high positive value, but SST tendency-heat flux correlation is very small. See their Figures 1A,B for details. However, these arguments were obtained based on monthly-mean time series and a set of predefined values for the parameters α, β, γ_a, γ_o, ω_a and ω_o. For submonthly scales, it is unknown how the model would behave. We therefore solve the model following the method of Bishop et al. (2017), but with hourly interval and parameters more suitable for the submonthly scales. The resultant lead-lag correlation between SST/SST tendency and heat flux in pure oceanic and atmospheric forcing cases are shown in Figures 3A,B. The values of the parameters can be found in the figure caption. It is inferable that the ocean-forcing case is still characterized with strong positive SST-heat flux correlation at zero lag, reflecting the small thermal inertia of the atmosphere. In the atmosphere-forcing case, however, there exists a peak of medium-strength negative correlation between SST and heat flux when the latter leads by about a day. The SST tendency-heat flux correlation is strong and negative when heat flux is in lead of an hour which is hardly visible from the figure. This indicates that the heat obtained from the atmosphere quickly forces the SST tendency, but it takes the SST about a day to respond to this forcing. Simultaneous SST-heat flux correlation is weaker than, but of the same sign as, the peak correlation in this case.

It is further confirmed that the result of Figure 3A is robust regardless of the values of the parameters, but Figure 3B, the pure atmospheric-forcing case, depends on γ_a and γ_o. The sensitivities of peak and simultaneous SST-heat flux correlations in the atmospheric-forcing case are then tested against wide ranges of γ_a and γ_o, and are shown in Figures 3C,D (note that in the figures 1/γ_a and 1/γ_o are shown). As 1/γ_o increases, peak SST-heat flux correlation becomes weaker, but remains statistically significant, and levels out for 1/γ_o ≥ 100 d. The SST response time, i.e., the time lag of the peak, increases slightly, but is mostly less than 2 days. The simultaneous correlation is always smaller than the peak value but of the same sign, weakening with increasing 1/γ_o and losing statistical significance for 1/γ_o ≥ 20 d. With increasing 1/γ_a, peak and simultaneous correlations both strengthen and remain significant, and the SST response time is almost always smaller than 2 days. Because the simultaneous correlation of SST-heat flux is always of the same sign as the peak value, and because its sign reverses in the pure atmospheric-forcing and pure oceanic-forcing cases, theoretically it can be used as a criterion for the direction of air-sea thermal forcing. For instance, Gulev et al. (2013) analyzed the simultaneous correlation between monthly SST and turbulent heat flux in the midlatitude North Atlantic ocean interior and found that the correlation is positive and strong on interdecadal timescales but negative and of medium strength on shorter timescales. Thereby they concluded that it is the ocean driving the atmosphere on interdecadal scales and the other way around on shorter scales, supporting the conjecture of Bjerknes (1964). Our solution using hourly interval confirms that this criterion can be extended to such small scales as submonthly, but only if the model itself is accurate. SST tendency-heat flux correlation is also sometimes used, for instance in Wu et al. (2015) who studied 10–20 and 30–60 days air-sea interaction in tropical and subtropical western North Pacific.

The pure atmospheric and oceanic forcing cases analyzed so far are extreme scenarios. In reality ω_a and ω_o should both be non-zero, corresponding to a mixture of forced and intrinsic SST variability. The 6 parameters must also exhibit spatial, temporal, and scale dependence. We next examine the characteristics of air-sea interaction on submonthly timescales in the framework of the theoretical model, but with objectively obtained parameters based on observations. To do so, using ERA5 SST and 2 m air temperature filtered for the 2–10 and 10–30 days scales, we estimate the parameters β and γ_o for each grid point in the research domain by regressing local values of dT_o/dt against T_a−T_o and T_o. The residual has normal distribution by definition, and its standard deviation is assigned to ω_o. Similarly, α, γ_a and ω_a are determined for each grid point. The parameters are then substituted to the model which is solved numerically for T_a and T_o. Results show that the spatial pattern of ω_o on the two scales (not shown) well resemble the standard deviation of SST in Figures 2E,F, indicating potential control of intrinsic oceanic perturbation on SST variability. The other parameters exhibit large ranges including negative values, resulting in unphysical temperatures over most of the domain (not shown). It is thus inferable that the theoretical model is problematic in the real ocean on submonthly timescales, most likely due to the ignorance of potentially important processes on these scales, e.g., wind mixing, convective mixing etc.

Next, we calculate SST-heat flux correlation directly from observations in ERA5. Here the heat flux is the net heat flux including all four components, i.e., sensible, latent, thermal radiation, and shortwave radiation. The simultaneous SST-heat flux correlation on 2–10 days timescale is shown in Figure 4A, with positive correlation found over much of the domain, except for the western subtropics and in the Japan Sea. In the KOE region the correlation is weak and does not exhibit obvious time-lag dependency (Figure 4B). This could probably mean that atmospheric and oceanic forcings are on the same order. On 10–30 days timescale, correlation in the KOE region becomes stronger and negative (Figure 4C), but the maximum negative correlation is found when heat flux leads by 2 days (Figure 4D). This behavior resembles the theoretic SST-heat flux correlation in the atmospheric-forcing case (Figure 3B). However, the SST tendency-heat flux correlation in Figure 4D has a negative peak when heat flux delays, in contrast to the hour-scale lead time by the heat flux as shown in Figure 3B. This indicates that the observed SST-heat flux interaction is indeed more complex than predicted by the theoretic model. A more sophisticated tool is thus required to explicitly assess the forced and intrinsic SST variability.

The single-column ocean model GOTM is employed to evaluate the local SST variability forced by atmospheric perturbations. The model is forced by air-sea heat fluxes (sensible, latent, thermal radiation, and shortwave radiation), freshwater fluxes (precipitation-evaporation), and momentum fluxes (wind stress). Air-sea coupling is apparently not possible in this model, but this is deliberately chosen to obtain purely forced SST variability rather than mixed forcing and response that are hard to isolate as in the theoretical model and the observations. The model is run separately for each grid point 1.5° apart in the domain. Unlike the theoretical model which only considers the ocean surface, oceanic dynamic and thermodynamic processes in the vertical dimension, including turbulent mixing and convection, are fully captured in the model. The modeled SST variability, therefore, is forced purely, directly, and locally by the overlying atmosphere, adjusted by mechanical and buoyant mixing in the water column. On the other hand, the part of SST variability missing in the model is attributed to intrinsic local and advective oceanic processes. A set of experiments are carried out using the GOTM model, each with different combination of forcing terms. The FULL experiment, e.g., is driven by all of the 6 forcing terms, whereas the BUOY experiment uses time-varying buoyancy fluxes and constant wind stress. Experiments other than the FULL experiment are called the sensitivity experiments. See Table 1 for a summary of all the experiments.

In Figure 5A we show the ratio between the STDs of 2–10 days SST simulated by the FULL experiment and observed in ERA5. High ratios indicate higher fraction of forced SST variability than intrinsic oceanic variability, and vice versa. Since this is a forced model, it is expected that the simulated SST variability would be stronger than in a coupled model, because the SST cannot feedback to the fluxes and damp them. The STD ratio, therefore, should be interpreted as an upper bound of atmospheric-forced fraction. Statistically, the STD ratio equals the square root of the reciprocal of the F-statistics which is defined as the variance ratio between two samples. Lower (higher) STD ratio than a critical value, therefore, indicates that the modeled and observed SST has significantly different (same) variance. It is found that on the submonthly timescale, forced SST variability dominates in the western subtropics and some spots to the east and north of the KOE, in accordance with recent findings of Yu et al. (2020). In the KOE, forced variability explains a very small fraction of the total variability, suggesting the dominant role of intrinsic SST variability there. On the 10–30 days scale, as shown in Figure 5B, the percentage of forced variability in the KOE region remains very low (<10%), but outside the KOE it is even higher (> 90%). In Figures 5C,D correlations between observed and modeled SST are examined. High and significant correlations are found in the subtropical band (25–30°N) on both timescales, roughly corresponding to the regions with high standard deviation ratio, indicating that the modeled SST there is close to observation in terms of both magnitude and phase. This is the relatively quiescent region between the strong jets of the Kuroshio and its extension and the weaker Subtropical Countercurrent, both abundant in mesoscale and submesoscale processes (Yang et al., 2013; Qiu et al., 2014, 2020; Zhou et al., 2021). GOTM’s ability in simulating the forced SST variability is thus justified in regions where oceanic advection and high-frequency perturbations and are weak. We are hence able to conclude in a quantitative manner that in the strong current regions of the KOE, the dominance of oceanic internal forcing is valid on timescales as short as 2–10 days, thus extending the results of Smirnov et al. (2014) and Bishop et al. (2017). However, it should be noted that the finding based on Figure 5 contradicts Figure 4, which showed a hint of either mixed forcing or pure atmospheric forcing in the KOE according to the theoretical model. This further disqualifies the theoretical model in tackling air-sea interaction on such small timescales.

To have a better idea on which processes are the most important, the results of the sensitivity experiments are useful. First, the SHF experiment driven solely by sensible heat flux variations exhibits different patterns of forced SST variability than in FULL (Figures 6A,B), but the percentage is still low in the KOE, except for in a small area very close to the Japanese coast. The inclusion of the latent heat flux in the THF experiment results in a STD ratio pattern resembling the FULL experiment (Figures 6C,D), indicating the importance of latent heat flux. Adding thermal radiation to THF, the NSF experiment keeps the same pattern of explained standard deviation but with reduced power (Figures 6E,F), presumably due to the damping effect of thermal radiation on SST. The NSF experiment includes the same air-sea fluxes as the theoretical model, only with the addition of vertical processes in the mixed layer. The difference between the two, therefore, reveals that only considering the sea surface is not sufficient in air-sea heat exchange on submonthly scales, and that subsurface mechanical and buoyant mixing plays a crucial role, especially in areas with deep mixed layer. Note that here the mechanical mixing is independent of wind stress. The NHF (Figures 6G,H) and BUOY (Figures 6I,J) experiments look very similar to NSF, indicating the negligible contribution from solar radiation (cloud cover) and net precipitation. Comparing the BUOY and FULL experiments, it can be found that the presence of wind stress forcing enhances the forced SST variability slightly in the KOE region and significantly in regions outside the KOE. Based on these results, three important processes have thus been identified, i.e., latent heat flux, subsurface heat exchange, and wind stress forcing (only for outside KOE). The latter two are not included in the theoretical model.

Further, the usability of correlation as an indicator for air-sea interaction is verified with the cross-correlation between the simulated SST and the net heat flux. Because the modeled SST variability is purely forced by the atmosphere, we search for the first peak (positive or negative) of the cross-correlation with increasing lead time by the heat flux. It is found that for the vast majority of the grid points in the domain, peak correlation is found at lead time of 1 day for the 2–10 days timescale, and 3 days for the 10–30 days timescale in every experiment (not shown). The spatial pattern of the peak correlation in the NSF experiment is shown in Figures 7A,B. On the 2–10 days timescale, negative correlation is found everywhere in the research domain, but on the 10–30 days scale, positive correlations occur. The possibility of positive delayed correlation even in such pure atmospheric forcing conditions is in sharp contrast to the theoretical model which requires negative delayed correlation (Figure 3B). This again suggests the deficiency of the theoretical model on these scales due to the ignorance of ocean vertical heat exchange. THF and NHF experiment results resemble the NSF (not shown), but the BUOY experiment exhibits different patterns (Figures 7C,D). The FULL experiment is similar with BUOY but with stronger positive correlations (Figures 7E,F), highlighting the role of wind stress mechanical forcing. It is therefore conclusive that the existence of vertical heat exchange and wind stress forcing in the real ocean can alter the sign of SST anomalies driven by atmospheric perturbations, thus making the SST-heat flux correlation not an appropriate tool for diagnosing the direction of air-sea interaction.

The previous generation of ECMWF reanalysis dataset, the ERA-Interim (Dee et al., 2011), went through enhancements of SST spatial resolution (but not temporal resolution), which has been shown to exert significant impacts on the overlying atmosphere (Zhang et al., 2020). However, potential effects of the refinement of SST temporal resolution in ERA5 (Figure 2A) are still unknown. Here, we take the view that the changed temporal resolution provides the possibility to examine the role of different levels of SST temporal variability on air-sea interaction. The above theoretical, observational, and modeling analyses are thus repeated on winter 2001/02 to compare with winter 2008/09. Figures 8A,C show the simultaneous correlation between SST and heat flux, and Figures 8B,D show the cross-correlation at 145°E, 38°N. The correlation on the 2–10 days scale is still positive but very weak, indicating mixed atmospheric and oceanic forcing. On the 10–30 days scale, unlike the medium strength negative correlation in winter 2008/09 (Figures 4C,D), the correlation is now positive and strong in the KOE. The cross-correlations of 10–30 days SST-heat flux and SST tendency-heat flux are in qualitative accordance with the oceanic-forcing case in the theoretical model (Figure 3A). The change from the atmospheric-forcing regime when 10–30 days SST variability is weak to the oceanic-forcing regime now indicates the sensitivity of the correlation tool to strength of SST disturbance. GOTM model results of the FULL experiment shown in Figures 8E,F suggest that atmospheric-driven SST variability accounts for a larger (smaller) fraction of the total variability on 2–10 days (10–30 days) timescales over the entire domain, most probably due to the weaker (stronger) total variability compared to winter 2008/09, since the forced variability is slightly weaker in winter 2001/02 on both scales (not shown). In the KOE region, however, the percentage explained by forced variability is still very small, favoring dominance of intrinsic SST variability on the two scales. Delayed correlations between SST and heat flux still exhibit positive values in the purely forced model simulations (Figures 8G,H), again disqualifying correlation as an appropriate indicator for air-sea interaction on these scales.

The same analysis is performed again using the 2.5°-resolution NCEP/NCAR Reanalysis 1 dataset for winter 2008/09. SST in the NCEP/NCAR data is from the Reynolds Analysis (Smith et al., 2008) with time interval of 1 week. The 2–10 days SST variability, therefore, is poorly resolved in this dataset (Figure 9A), much like the situation of ERA5 before 2007. However, although stronger than the 2–10 days band, the 10–30 days variability (Figure 9B) is also underrepresented since a large fraction of 10–30 days SST variability is associated with mesoscale perturbations which require a much finer spatial resolution (Zhou et al., 2021). The strength of purely forced SST variability simulated by the GOTM-FULL experiment (Figures 9C,D) resembles the ERA5 results, again suggesting the importance of mixed layer processes in suppressing the direct surface response to atmospheric forcing on the submonthly scales, regardless of the spatial scale. Delayed correlations between NHF and SST (Figures 9E,F) are now negative everywhere, conforming to the pure-atmospheric driven scenario. Given the importance of wind stirring revealed above in section “General Ocean Turbulence Model Results,” this is indicative that subgrid-scale (i.e., mesoscale) anomalies of wind stress might be responsible for the reversal of NHF-SST correlation.