The best track data of the typhoon, including its center location, pressure at sea level, and maximum sustained wind speed at 6-h intervals, were provided by the Joint Typhoon Warning Center (JTWC). SST data were obtained from the Tropical Rain Measuring Mission’s Microwave Imager (TMI), which has a daily temporal resolution and a spatial resolution of 0.25°×0.25°. The microwave measurements applied by TMI SST data are not contaminated by clouds; thus, this product has been widely used in previous studies for assessing the SST response to TCs (e.g., Vincent et al., 2012a; Mei and Pasquero, 2013). Argo profile data, including measurements of temperature and salinity, were obtained from the Coriolis Global Data Acquisition Center of France (Argo, 2000). Sea level anomaly (SLA) data were obtained from the Archiving, Validation, and Interpretation of Satellite Oceanographic (AVISO) and had the same temporal and spatial resolutions as the SST data. The SLA data are based on radar altimeters, which are also not limited to cloud-free regions. The SLA data were used to detect eddies following the methods of Zhang et al. (2018), who modified the algorithm from Chaigneau et al. (2009) and Chelton et al. (2011).

The coupled ocean-atmosphere-wave-sediment transport (COAWST) modeling system (Warner et al., 2008; Warner et al., 2010) was applied with oceanic and atmospheric components of the Regional Ocean Modeling System (ROMS, version 3.7) (Shchepetkin and McWilliams, 2005) and the Weather Research and Forecasting (WRF) model (version 4.0.3) (Skamarock et al., 2008), respectively. The Model Coupling Toolkit was used to transform the variables between the different components in the COAWST modeling system (Larson et al., 2005). In particular, the WRF model transferred sea surface wind stress, latent heat, sensible heat, shortwave radiation and longwave radiation data to the ROMS model. Meanwhile, the WRF model received the SST, which was continuously updated by the ROMS model. The time interval of exchange between the two models was 120 seconds. The COAWST modeling system has been widely used in many studies to investigate air-sea interaction processes during TCs (Zhao and Chan, 2016; Prakash and Pant, 2017; Wu et al., 2018).

The spatial domain for both the WRF and ROMS models was 116°E to 150°E and 9°N to 33°N ( Figure 1A ), with a horizontal resolution of 0.05×° 05°. Thirty-six and 40 vertical sigma levels were set in the WRF and ROMS models, respectively. The topography in the ROMS model was based on ETOPO1 (https://www.ngdc.noaa.gov/mgg/global/global.html). The WRF Single-Moment 6-class scheme was chosen as the microphysics scheme. Longwave and shortwave radiation was computed via the Rapid Radiative Transfer Model and Dudhia scheme. The planetary boundary layer scheme and cumulus parameterization used in the WRF model were the Yonsei University scheme and Kain-Fritsch scheme, respectively. A generic length-scale mixing closure was applied in the ROMS vertical mixing parameterization scheme. Third-order upwind horizontal advection and fourth-order centered vertical advection were adopted for the temperature and salinity simulations.

The initial and boundary data of the WRF were from the National Centre for Environmental Prediction Final Analysis. To identify preexisting mesoscale cyclonic eddies, the initial conditions, e.g., temperature, salinity, current and water level, for the ROMS were obtained from the Hybrid Coordinate Ocean Model (HYCOM) and Navy Coupled Ocean Data Assimilation (NCODA) system (Cummings, 2005). HYCOM is a primitive equation ocean general circulation model with a resolution of 0.08° in the zonal direction and a nonuniform resolution in the meridional direction, i.e., 0.08° from 40°S to 40°N and 0.04° in the remaining regions, which can resolve mesoscale eddies. Through the NCODA system, the assimilation of available satellite altimeter observations, SST, and in situ vertical temperature and salinity profiles simulates the preexisting mesoscale eddies consistently with the observations. The lateral boundary information in the ROMS was also obtained from the HYCOM.

ROMS was initialized to start on July 5, 2013, with a simple setup of a 3-day model time (July 5-7) before the formation of Typhoon Soulik as a tropical storm on July 8. WRF was initialized to start on July 8 and was instantly coupled with the ROMS model. The control run simulation was conducted by incorporating near real-time eddies from the HYCOM, which was referred to as EXP_HYCOM.

To examine the effect of preexisting CEs on the upper ocean response to Typhoon Soulik and the underlying mechanism, we conducted an experimental simulation where the preexisting CEs were not included. The oceanic initial fields on July 5, 2013, were replaced by seven-year (2010-2016) averaged climatology fields in July from HYCOM, which was referred to as EXP_CLIMAT.

To diagnose the oceanic response, a temperature budget analysis was performed based on the following equation (Wei et al., 2014):

where the left-hand term is the temperature change rate and the right-hand terms are advection (1), horizontal diffusion (2), and vertical diffusion (3). In the ROMS model, the vertical boundary condition at the sea surface (z = ζ(x, y, t)) for temperature is K h ∂ T ∂ z =   Q ρ 0 C p , where the surface heat flux is contained in the vertical diffusion (vertical mixing) term. Each term on the right-hand side of equation (1) can be integrated over time to obtain its accumulated contribution to the change in temperature since 0000 UTC on July 8.

The remote sensing of SST during Typhoon Soulik ( Figure 3 ) revealed that the SST near its track was mostly higher than 29°C on July 6 and 7, 2013. The SST slightly dropped (by less than 1°C) to the north of the track on July 8, and significant cooling occurred later on July 10. Although sea surface cooling took place along the entire track, a detailed examination found that the cooling was not homogeneously distributed but presented a multicore structure, identified as several spatial peaks of sea surface cooling, with a minimum SST of less than 24°C. During the period from July 10 to July 12, the maximum wind speed of Typhoon Soulik gradually decreased in association with a stable translation speed, i.e., approximately 6 m/s. Wei et al. (2017) determined that large-scale sea surface cooling occurred during the passage of Typhoon Soulik with a prognostic algorithm where only wind was incorporated; however, multicore sea surface cooling was captured by a machine learning algorithm where ocean parameters, e.g., SST, sea surface height, and temperature at 100 m depth, were incorporated. Considering the inconsistent patterns of typhoon intensity and translation speed in relation to the spatial features of sea surface cooling, the multicore sea surface cooling was attributed to other factors, e.g., oceanic conditions.

The SLA before, during and after Typhoon Soulik ( Figure 4 ) captured several eddies on July 6, 2013, before the formation of the typhoon. Three CEs, located between 125°E and 139°E, stayed at roughly the same positions before, during and after the passage of the typhoon. The maximum sea surface cooling induced by Typhoon Soulik ( Figure 5 ) was the minimum value in the SST anomaly, defined as the difference with that at each location during its passage on July 8. The cold core features became more prominent in the anomalous field, and their locations were consistent with those of the three cold eddies (CE1, CE2, and CE3) ( Figure 5 ). Thus, the cold cores induced by Typhoon Soulik were associated with the preexisting CEs, and the spatial peaks with the maximum sea surface cooling within the three cold eddies, referred to as the cooling centers, were marked as C1, C2 and C3 ( Figure 5 ). The SST anomalies at the C1, C2 and C3 cooling centers were -8.10°C, -7.95°C and -6.60°C, respectively.

The upper ocean response to Typhoon Soulik was simulated by the coupled atmosphere and ocean model. The sea level on July 8 in the EXP_HYCOM simulation ( Figure 6A ) included eddies, and the spatial pattern of the CEs was highly consistent with the AVISO data ( Figure 4 ). Similarly, the maximum sea surface cooling for each grid in EXP_HYCOM ( Figures 6C ) showed a multicore structure similar to that of the remote sensing of SST ( Figure 5 ), indicating that the cooling associated with the CEs was successfully reproduced by the model. Furthermore, the cooling centers in the model (M1, M2, M3) that occurred within the CEs were highly consistent with those in the altimetry data (C1, C2, C3) ( Figure 5 ). The maximum sea surface cooling at the M1, M2 and M3 cooling centers was -6.73°C, -7.35°C and -8.31°C, respectively.

The sea level on July 8 in the EXP_CLIMAT simulation did not capture the occurrence of CEs ( Figure 6B ), and there was no multicore sea surface cooling. The maximum sea surface cooling showed a continuous band pattern ( Figure 6D ) with a weaker value (-4.85°C) compared with that in EXP_HYCOM (-8.31°C) and in the observation (-8.1°C). In EXP_CLIMAT without the CEs, the maximum SST cooling at the same positions as M3, M2 and M1 were only -2.11°C, -4.68°C and -4.48°C, respectively. As the CEs were not present in the initial field, the continuous band of the SST cooling response in EXP_CLIMAT indicated that the three cold cores induced by Typhoon Soulik in EXP_HYCOM were related to the preexisting CEs. Considering that the sea surface cooling at M2 in EXP_HYCOM was closer to the corresponding observations than those at M1 and M3, the following discussion focuses on the dynamic mechanisms of CE2.

The mixed layer depth and the strength of stratification beneath it were found to be the two major oceanic parameters influencing the wind-induced cooling amplitude (Vincent et al., 2012b). The vertical temperature profiles at M2 showed similar mixed layer depths before the passage of Typhoon Soulik in EXP_HYCOM and EXP_CLIMAT ( Figure 7 ), but the temperature gradient below the mixed layer was larger in EXP_HYCOM than in EXP_CLIMAT. A large temperature gradient in a CE makes vertical mixing and upwelling more efficient, facilitating sea surface cooling (Wu et al., 2007; Zheng et al., 2010). Consistently, the colder water entrained from the upper thermocline was demonstrated to enhance sea surface cooling when there was a large temperature gradient at the base of the mixed layer (Wu et al., 2007).

A region around M2, which was 200 km along the track and within 200 km from the TC track on both the right and left sides (see the black box in Figure 6 ), was chosen for further dynamic analysis. Typhoon Soulik passed the region at approximately 0100 UTC on July 11, 2013. The maximum sea surface cooling occurred 18 and 20 hours after the typhoon passage in EXP_HYCOM and EXP_CLIMAT, respectively. The oscillation period of the temperature anomaly was close to the inertial period at M2, i.e., 30 hours, indicating the existence of near-inertial oscillations (Price, 1983; Gill, 1984; Shay and Elsberry, 1987). To eliminate the influence of the near-inertial oscillations on the contribution of mixing and advection to SST anomalies, the sea surface cooling was obtained by averaging the SST anomalies during one inertial period from 3 to 33 hours after the typhoon passage ( Figures 8A, B ). The SST anomaly without CE2 in EXP_CLIMAT was relatively uniform along the track ( Figure 8B ), which was in agreement with previous research where the composite method was applied along the track (e.g., Mei and Pasquero, 2013). In contrast, several semicircular SST anomaly contours (black contours) between 100 km and 200 km to the right of the track in EXP_HYCOM illustrated a prominent along-track variation in the SST anomaly, probably leading to a cooling core. The SST anomaly contours to the right of the track had the same shapes as the outermost zeta contour (red dashed curve) in EXP_HYCOM ( Figure 8A ), demonstrating that the along-track variation was associated with the preexisting CE2.

Both experiments showed rightward-shifted patterns of sea surface cooling ( Figures 8A, B ), which is a well-known feature of the SST response to typhoons in the Northern Hemisphere. The clockwise rotation of the wind stress tends to resonate with the wind-induced inertial currents, enhancing vertical mixing and resulting in a rightward strengthening of surface cooling (Price, 1981). This pattern was prominent in EXP_CLIMAT ( Figure 8D ), where the maximum sea surface cooling induced by mixing (-4.4°C) appeared 50 km to the right of the track. However, there was no obvious signal in EXP_HYCOM ( Figure 8C ); instead, the maximum sea surface cooling induced by mixing (-6.0°C) appeared 40 km to the left of the track. The mean sea surface cooling within 200 km that was induced by mixing in EXP_HYCOM was 1.5°C higher than that in EXP_CLIMAT. Thus, the rightward-shifted cold wake in EXP_HYCOM was attributed to horizontal advection rather than mixing.

The SST anomaly and corresponding contribution of a long-track advection in EXP_HYCOM showed a similar cambered shape on the right side of the TC track ( Figures 8A, E ), indicating that the cooling pattern in EXP_HYCOM was mainly caused by a long-track advection. However, the contribution of a long-track advection was relatively uniform along the track in EXP_CLIMAT and was negative (positive) on the right (left) side ( Figure 8F ). The advection-induced cambered SST anomaly in EXP_HYCOM was consistent with the outermost zeta contour ( Figure 8E ), indicating the impact of the geostrophic eddy current. The contribution of the cross-track advection in both experiments showed the same pattern: positive in the central area and negative in the surroundings. In contrast with the result in EXP_CLIMAT, EXP_HYCOM also presented a more important contribution of cross-track advection due to the influence of CE2, although along-track advection was more prominent. The contribution of vertical advection was negligible in both experiments for the change in SST (not shown), which was consistent with the results of Lin et al. (2017).

Overall, the contributions of mixing, along-track advection and cross-track advection to the SST anomaly were all enhanced by CE2 ( Figure 8 ), resulting in an intensified rightward-shifted cold wake during Typhoon Soulik. The enhanced SST anomaly that was induced by mixing mainly occurred in the central area, i.e., within 50 km of the track of the typhoon, without the rightward-shifted signal. Thus, eddy-related advection played a dominant role in the rightward-shifted cold wake in EXP_HYCOM.

A temperature budget analysis at M2 was conducted to further quantify the contributions of mixing and advection to the rightward-shifted cold wake ( Figures 9 , 10 ). In both EXP_HYCOM and EXP_CLIMAT, strong sea surface cooling, e.g., a temperature anomaly of more than -1°C, appeared immediately before the passage of the typhoon and rapidly deepened to approximately 50 m within several hours ( Figures 9A, B ). The sea surface cooling was -6.94°C in EXP_HYCOM ( Figure 9A ), while the contributions of vertical mixing and advection were -4.24°C and -2.71°C, respectively ( Figures 9D, G ). Thus, vertical mixing dominated the SST response, with a contribution of 61% to sea surface cooling, which is in agreement with the results of previous studies (e.g., Price, 1981; Lin et al., 2017). The accumulated contribution of total advection to the surface temperature anomaly gradually increased as Typhoon Soulik approached, reaching its maximum on July 11, right after the typhoon passed M2. Advection contributed 39% to sea surface cooling; in contrast, horizontal diffusion had a negligible contribution (not shown).

The sea surface cooling in EXP_CLIMAT was -4.26°C ( Figure 9B ), and the contributions of vertical mixing and advection were -4.30°C and 0.06°C, respectively ( Figures 9E, H ). Compared to EXP_HYCOM, vertical mixing played a more important role in EXP_CLIMAT by contributing all of the sea surface cooling. Without preexisting CEs, advection was driven by typhoon-generated currents, which suppressed cooling by introducing a small positive temperature anomaly.

Under the influence of the typhoon, the sea surface cooling in EXP_HYCOM was 2.68°C greater than that in EXP_CLIMAT ( Figure 9C ). The vertical mixing in EXP_HYCOM actually caused a smaller temperature anomaly in sea surface cooling compared to that in EXP_CLIMAT ( Figure 9F ). Instead, the total advection in EXP_HYCOM induced a much larger decrease in SST (-2.71°C) than that (+0.06°C) in EXP_CLIMAT ( Figure 9I ). The advection in EXP_HYCOM was induced by a combined impact of typhoon-generated currents and eddy-related currents, and the latter processes were not incorporated in EXP_CLIMAT. Thus, advection due to eddy-related currents played a dominant role in enhancing sea surface cooling by introducing an additional negative temperature anomaly at M2.

The representativeness of the cooling pattern at M2 was further assessed over an area of 1°, i.e., 111 km in the along-track direction and 111 km in the cross-track direction, centering around M2. During the passage of the typhoon, the area-averaged sea surface temperature anomaly in EXP_HYCOM (-6.27°C) was -2.20°C lower than that in EXP_CLIMAT (-4.07°C). Correspondingly, the enhancement of sea surface cooling due to vertical mixing (-4.09°C and -3.31°C) and advection (-2.18°C and -0.76°C) contributed 35% and 65% in EXP_HYCOM and EXP_CLIMAT, respectively. The advection still played a dominant role in the eddy-induced enhancement of the area-averaged sea surface cooling at M2, indicating that M2 can represent the cooling core area related to the eddy-induced dynamics.

The advection effect − ( u ∂ T ∂ x + v ∂ T ∂ x + w ∂ T ∂ z ) consists of three terms: − u ∂ T ∂ x (along-track advection), − v ∂ T ∂ x (cross-track advection), and − w ∂ T ∂ z (vertical advection), and their contributions at M2 are shown in Figure 10 . In EXP_HYCOM, along-track advection (-2.72°C) contributed a majority of the temperature anomaly (-2.71°C) induced by advection ( Figure 10A ). The cross-track advection contributed a negligible role of 0.01°C ( Figure 10D ). In EXP_CLIMAT, cross-track advection caused a positive temperature anomaly (0.89°C) that was mostly offset by along-track advection (-0.84°C) ( Figures 10B, E ). In both EXP_HYCOM and EXP_CLIMAT, the contribution of vertical advection to sea surface cooling was negligible due to the small vertical temperature gradient and vertical current in the mixed layer ( Figures 10G, H ). In comparison, upwelling caused an obvious negative temperature anomaly in the subsurface at an approximately 100 m depth due to the large vertical temperature gradient in the thermocline ( Figures 10G, H ). However, this negative temperature anomaly was mostly counteracted by the positive temperature anomaly that was induced by vertical mixing ( Figures 9D, E ).

A comparison of the three advection terms between EXP_HYCOM and EXP_CLIMAT is shown in Figures 10C, F, I . At the forcing stage (-0.5 to 0.5 days, Sanford et al., 2011), along-track advection was almost the sole advection process that directly induced a negative temperature anomaly ( Figure 10C ). The cross-track advection enhanced sea surface cooling after the forcing stage ( Figure 10F ). The along-track and cross-track advection terms contributed -1.88°C (68%) and -0.89°C (32%), respectively, to the enhancement of the average sea surface cooling. The difference in vertical advection was insignificant since the process offered a limited contribution to sea surface cooling in both EXP_HYCOM and EXP_CLIMAT.

The averaged along-track and cross-track components of the current at M2 from July 8 to July 10, 2013, were both negative (in the moving direction of the typhoon or to the left of the typhoon track) with values of -0.24 m/s and -0.05 m/s, respectively ( Figures 11A, B ). For a CE, negative along-track and cross-track components of the current indicate that it is located to the northwest of the eddy (Yang et al., 2013; Zhang et al., 2018). The corresponding along-track temperature gradient ∂ T ∂ x is negative, and the cross-track temperature gradient ∂ T ∂ y is positive. Specifically, the average ∂ T ∂ x at M2 was -5.7×10^-6°C/m, and the average ∂ T ∂ y was 5.5×10^-6°C/m in EXP_HYCOM during this period. Therefore, the along-track advection − u ∂ T ∂ x was negative, and the cross-track advection − v ∂ T ∂ y was positive, which is in agreement with the model results at M2 before the approach of the typhoon ( Figure 10A ). Since the magnitudes of the along-track current and temperature gradient were higher than their corresponding cross-track values, along-track advection was dominant before July 10, resulting in a negative temperature anomaly due to eddy advection. The advection-induced temperature anomaly averaged from 0000 UTC on July 8 to 0000 UTC on July 10 at M2 in EXP_HYCOM was -0.082°C.

On July 10, a cold wake appeared to the east of M2 ( Figure 3 ) that increased the magnitude of the negative along-track temperature gradient ∂ T ∂ x at M2. Together with a strengthened along-track negative velocity ( Figure 11 ), it induced a negative along-track advection effect − u ∂ T ∂ x , resulting in enhanced sea surface cooling at M2.

In this study, the along-track current located to the northwest of CE2 and to the right side of Typhoon Soulik also produced forward advection of the cold wake. Thus, the along-track advection enhanced the typhoon-induced sea surface cooling at M2 by strengthening the westward current and increasing the negative along-track temperature gradient. After July 10, although both the along-track velocity and temperature gradient fluctuated with time, the accumulated contribution of a long-track advection to the surface temperature anomaly was always negative ( Figure 10 ), persistently enhancing sea surface cooling.

Typhoon Soulik drove a negative ∂ T ∂ y at M2, associated with a positive cross-track component of the current immediately after its passage. Thus, the cross-track advection effect, − v ∂ T ∂ y , was positive ( Figure 10E ). In contrast, the preexisting CE2 drove positive ∂ T ∂ y and negative cross-track current components at M2 that weakened the negative temperature gradient and the positive current, respectively. As a result, the cross-track advection effect − v ∂ T ∂ y decreased (enhanced cooling) due to the effect of CE2 in EXP_HYCOM ( Figure 10D ), which was not reflected in EXP_CLIMAT ( Figure 10E ).

Although a larger temperature gradient was identified below the mixed layer at M2 within the preexisting CE in EXP_HYCOM ( Figure 7 ), vertical mixing did not drive a larger temperature anomaly for sea surface cooling compared with that in EXP_CLIMAT. In contrast, advection played a more important role in the enhancement of sea surface cooling. Considering the stable location of the eddy during the passage of the typhoon, the advection was mainly attributed to the eddy-related anticlockwise current instead of the eddy movement, as demonstrated in Ning et al. (2019).

To determine why the contribution of mixing to the surface temperature anomaly at M2 was smaller in EXP_HYCOM than in EXP_CLIMAT, the contribution of advection to the subsurface temperature anomaly and the vertical temperature gradient due to the preexisting CE was examined to understand the subsequent contribution of mixing to the surface temperature anomaly. When the advection of the cold wake associated with the preexisting CE enhanced the sea surface cooling at M2, it also affected the subsurface temperature anomaly ( Figures 9G, I ). In the subsurface layer (70-130 m), the accumulated temperature anomaly induced by advection was initially positive and rapidly became negative in both EXP_HYCOM and EXP_CLIMAT due to the passage of Typhoon Soulik. Vertical advection played a dominant role in the temperature anomaly in this layer ( Figures 10G, H ), while the contributions of a long-track and cross-track advections were very small ( Figures 10A, B, D, E ).

The vertical advection tended to strengthen the vertical temperature gradient by inducing subsurface cooling, which enhanced the mixing-induced sea surface cooling ( Figures 10G, H ). This was in agreement with the results of previous studies that more than half of the storm-core sea surface cooling was neglected if the effect of upwelling was not considered for slow-translating TCs (Yablonsky and Ginis, 2009; Chiang et al., 2011). However, there was little difference in the value of subsurface cooling induced by upwelling in EXP_HYCOM and EXP_CLIMAT ( Figures 10G–I ), indicating that CEs made little contribution to the change in the vertical temperature gradient induced by vertical advection.

The along-track advection reduced the vertical temperature gradient, as it induced obviously greater sea surface cooling compared with that in the subsurface ( Figures 10A, B ). Moreover, as the sea surface cooling induced by a long-track advection in EXP_HYCOM was more significant than that in EXP_CLIMAT, the reduction in the vertical temperature gradient induced by a long-track advection was larger in EXP_HYCOM than that in EXP_CLIMAT, indicating the modification effect of the preexisting CE on the vertical temperature gradient by a long-track advection.

The sea surface cooling induced by cross-track advection was larger than that at the subsurface in EXP_HCYOM, indicating a decrease in the vertical temperature gradient induced by cross-track advection ( Figure 10D ). However, there was only a small difference in the temperature anomaly between the surface and subsurface in EXP_CLIMAT ( Figure 10E ), implying a negligible contribution to the vertical temperature gradient. Therefore, compared to the result without CEs, preexisting CEs in EXP_HCYOM could decrease the vertical temperature gradient by cross-track advection.

Overall, preexisting CEs tended to reduce the vertical temperature gradient by horizontal advection, which may offset the larger temperature gradient at M2 within the preexisting CEs, resulting in less sea surface cooling induced by mixing with a larger temperature gradient.

Considering that the effect of eddy-related advection on the SST response was related to the relative locations of the mesoscale ocean eddies to the typhoon track, the eddy-related current could lead to forward advection of the typhoon-induced cold wake when the along-track component of the eddy-related current was consistent with the moving direction of a typhoon. The scenario was also valid for the effect of anticyclonic oceanic eddies; for example, when a preexisting oceanic anticyclonic eddy was located to the right of a typhoon track, the along-track current in the fourth quadrant of the eddy led to forward advection of the cold wake, enhancing sea surface cooling (Yablonsky and Ginis, 2013). On the other hand, Ning et al. (2019) found that CE-induced horizontal advection weakened sea surface cooling and increased subsurface warming. Due to the complex oceanic dynamics during the passage of typhoons and the coeffect of mesoscale eddies, further study is required for more cases with three-dimensional observations.