The results from the WRF model simulation of two sequential TCs Kyarr and Maha are presented in this section. The JTWC and ERA5 estimate basic tropical cyclone data at every 6h can be observed, though in situ observations during tropical cyclone events in the Arabian Sea region are lacking. The tracks which TCs Kyarr and Maha were observed and simulated are presented in Figures 1A, B . In ERA5, WRF model simulation, the centers of Kyarr and Maha are located at the minimum point of sea level pressure. In Figures 4A, B , the tracks of Kyarr, Maha provided by JTWC, ERA5, WRF model simulation all show that TCs first moved to the northwest and then changed abruptly. In general, the TC tracks provided by ERA5 are basically consistent with those provided by JTWC. As would be observed, there is a significant difference between the WRF model simulation and ERA5, JTWC. For example, the WRF model simulation track of TC Kyarr biased to the west side of the observed track at first. Then the track was greatly biased to the east side of the observed track after it’s sudden change. Though deviating from the Maha’s observed track at the beginning, it’s track changes were basically consistent with those provided by ERA5 and JTWC afterwards. In another study, it was found a higher track deviation might be attributed to the initial positioning error (Osuri et al., 2012). In general, the WRF model simulation exhibits relatively large errors in the track of TC Kyarr. Unlike TC Maha, the WRF model simulation has a relatively high consistency with the track provided by JWTC and ERA5. From the above, in this paper, we will use the tracks provided by JTWC as the tracks of TCs Kyarr and Maha.

The tracks, intensities and translation speeds of Kyarr and Maha in October 2019 are shown in Figures 1A–G . TCs Kyarr and Maha followed a similar path across the AS from October 24 to November 06, 2019. First, Kyarr originated and reached tropical storm status on the night of 25 October 2019 near India ( Figure 1A ). Then Kyarr moved north-westwards continuously developed from 00:00 UTC to 18:00 UTC on October 26, and it quickly strengthened from Category 1 to Category 4. The Category 4 stage persisted for approximately 54 hours according to the Saffir-Simpson scale.

At noon on October 29, Kyarr made a sudden south-westerly turn and gradually weakened before dissipating into a tropical storm and finally disappearing over the western AS. Afterwards, a week after Kyarr originated, Maha originated in the sea near southern India and moved north-westwards. Later, Maha was upgraded to Category 1 on the evening of November 02 and gradually strengthened to reach its maximum intensity of Category 3 in the early morning of November 04. Approximately a day after Maha continued moving northwest as a Category 3 TC, the track of Maha changed to the east near the area where Kyarr suddenly turned. Finally, Maha gradually weakened and disappeared over the northern AS.

For the purpose of comparing the WRF model simulation results, the minimum central pressure and maximum sustained wind speed from JTWC and ERA5 were used ( Figures 1C–F ). The maximum 10 m wind speed from the WRF model simulation and ERA5 were used as the maximum sustained wind speed.

Notably, the TCs intensity indicated by ERA5 appears to be particularly weak compared to the intensity indicated by JTWC and WRF model simulation. Moreover, the WRF model simulation and ERA5 peaks (maximum sustained wind speed and minimum central pressure) were delayed relative to those of JTWC. However, according to the evolution of minimum central pressure recorded by JTWC, WRF model simulation and ERA5 for Kyarr and Maha, the minimum central pressure values of TC Kyarr were 923 hPa, 955, and 974 hPa, respectively ( Figure 1C ). The JTWC-recorded minimum central pressure for TC Maha was 959 hPa, that recorded by WRF model simulation and ERA5 were 972 and 998 hPa, respectively ( Figure 1D ). For the maximum sustained wind speed of Kyarr and Maha, JTWC maximum sustained wind speed were faster than WRF model simulation and ERA5. And the JTWC-recorded maximum sustained wind speed for TC Kyarr was 69 m/s, while that of WRF model simulation and ERA5 were 45 and 29 m/s, respectively ( Figure 1E ). Moreover, the maximum sustained wind speed of TC Maha recorded by JTWC was 54 m/s, while that recorded by WRF model simulation and ERA5 were 39 and 20 m/s, respectively ( Figure 1F ). The WRF model simulation underestimated the minimum central pressure and maximum sustained wind speed during the development stage compared to that of JTWC but was closer to ERA5.

The translation speeds of Kyarr and Maha were calculated according to the locations of the TCs centres in the time series ( Figure 1G ). Kyarr moved relatively slowly (1–4 m/s) during most stages over the AS. Maha moved relatively fast during the early stage, but after November 02, Maha moved relatively slowly and attained a translation speed (U_h) of approximately 2.5 m/s.

As shown in Figure 3A , the background ocean environment on October 24 provided favourable conditions for TC generation, and the AS was covered by warm water with a high SST. During the passage of Kyarr and Maha, there was a cool trail along their tracks with a rightwards bias ( Figures 3B–H ). In particular, four distinct cold patches were detected, which are marked with black boxes and correspondingly labelled C1, C2, C3, and C4 in Figure 3 . On October 31 ( Figure 3D ), the first and second cooling patches appeared after the passage of Kyarr. Both cooling patches were located to the right of TC Kyarr’s track. In regions C1 and C2, the maximum decreases of SST were 3°C and 4.5°C, respectively ( Figure 3D ).

As seen in Figures 3E–H , the sea surface cooled again after the Maha passed through the same area 8 days later ( Figure 3H ). And two more new cooling patches emerged, regions C3 and C4, their temperature decreases of 5°C and 4.5°C, respectively. While the first cooling patch in region C1 continued to enlarge, the maximum SST drop in region C2 gradually recovered during the passage of TC Maha ( Figures 3E–H ). Thus, after the passage of Kyarr and Maha, there were three cooling regions on the sea surface. And the cooling in region C4 is significantly stronger than that in other regions.

Figure 5 shows the wind field simulated by the WRF model. The calculated wind stress curls were obtained from WRF model simulation and ERA5, respectively. A positive wind stress curl closely linked to Ekman pumping to generate upwelling and promote cold wake development (Chiang et al., 2011). The wind stress curls obtained by WRF model simulation and ERA5 calculation were all positive over the AS, but the wind stress curls according to ERA5 were relatively weak. The ERA5 spatial pattern of the wind stress curls were similar to that of WRF model simulation ( Figures 5A–C, G–I ), the magnitude was 2–6 times larger for WRF model simulation ( Figures 5D–F, I, J ). For example, at 00:00 UTC on October 27, the wind stress curls of Kyarr obtained from WRF model simulation were obviously four times stronger than that of ERA5. The maximum value of the wind stress curls from ERA5 was ~1.1×10^-5 N/m³ ( Figure 5A ) and that from WRF model simulation was~5×10^-5 N/m³ ( Figure 5D ). The maximum value of wind stress curls of Maha obtained from ERA5 was ~0.8×10^-5 N/m³ ( Figure 5I ) and that from WRF model simulation was ~4.8×10^-5 N/m³ ( Figure 5L ), at 18:00 UTC on November 04. The powerful wind of the TCs induced a series of physical processes, and the sea surface cooled, as shown in Figure 3 . As the SST decreased along the TC tracks analyzed above, it can be seen that the wind forcing according to ERA5 was too weak to produce this feature. In contrast, the WRF-obtained atmospheric forcing could produce this cold wake more realistically.

Figure 6 shows the time evolution of chl-a concentration at sea surface caused by Kyarr and Maha. The chl-a concentration was 0.5~1 mg/m³ in most areas of AS before two sequential TCs ( Figure 6A ). Unsurprisingly, the three typical regions with obvious sea surface cooling and affected by strong wind stress curls were also noted to exhibit pronounced increases in chl-a after two sequential TCs. During the passage of TC Kyarr, the chl-a concentration did not significantly increase along the TC track but showed only a small increase of 2 mg/m³ near the generation region of Kyarr ( Figure 6B ). After Kyarr and during Maha, there were three chl-a concentration regions: the region near the generation of Kyarr and regions C1 and C3. In these three regions, the chl-a concentration could reach ≥ 3.5 mg/m³ ( Figure 6C ). After the passage of Kyarr and Maha, the chl-a concentration in the region near the generation of Kyarr did not obviously increase. In region C1, the chl-a concentration gradually returned to ≤ 2.5 mg/m³. In region C3, the chl-a concentration increased to a maximum of 4.5 mg/m³, and another chl-a concentration region appeared, region C4. There was a significant increase in the chl-a concentration in this region, with a maximum value > 5 mg/m³.

In Figures 7A–F , the evolution of sea surface relative vorticity is shown. After TC Kyarr passed ( Figures 7A–C ), the three regions exhibited positive relative vorticity, regions C1, C3, and C4 reached approximately 1×10^-5 s^-1, 1.5×10^-5 s^-1, and 1×10^-5 s^-1, respectively. Interestingly, after TC Kyarr passed through, a positive relative vorticity also appeared in between of regions C3 and C4 ( Figure 7C ). However, the relative vorticity in this region was smaller than the relative vorticity in regions C3 and C4 and gradually weakened ( Figures 7D–F ). After Maha passed, the positive relative vorticity of regions C3 and C4 developed further and reached maximum values of 2×10^-5 s^-1 and 3×10^-5 s^-1, respectively. The relative vorticity of region C1 did not change significantly, but the area with positive relative vorticity in this region increased ( Figures 7D–F ). Because cyclonic eddies are closely linked to the positive relative vorticity, further research on this topic is presented in the next section.

Figure 8 depicts the temporal evolution of the SSHA before, during and after the passage of Kyarr and Maha. At first, no SSHA existed before the two sequential TCs ( Figure 8A ). Then a SSHA began to appear during the passage of TC Kyarr( Figure 8B ). Meanwhile, a cold eddy was first identified by eddy detection algorithms in the area during the passage of TC Kyarr ( Figure 8B ). Then, as Kyarr continued to advance, a second SSHA appeared on the sea surface. And a second cold eddy was also identified ( Figure 8C ). After the passage of Kyarr, the sea surface height in C1, C3 and C4 regions decreased by -10 cm, -15 cm and -15 cm, respectively. And a third cold eddy was identified at location of the turn in Kyarr’s track ( Figure 8D ). In the three cold eddy areas, the sea surface geostrophic currents gradually strengthened and showed the characteristics of cyclone ( Figures 8B–D ). Moreover, during the passage of Maha, negative SSHAs were further reinforced. In region C1, the sea surface height decreased by approximately -15 cm. In regions C3 and C4, the sea surface height extensively decreased, and both reached -25 cm. The sea surface geostrophic currents formed more intense cyclonic features in these regions, and the highest velocities were found around the edges of the three cold eddies ( Figures 8E, F ). Finally, the shape of the three cold eddies stabilized on November 07 ( Figure 7F ).

After the two TCs passed through on November 07, three stable cold eddies were identified by eddy detection algorithms. The centre of cyclonic eddy e1 in region C1 was located at 68.125°E, 16.875°N and had a radius of up to 69 km and an area of up to ~15000 km². The centre of cyclonic eddy e2 was located at 65.375°E, 18.375°N and had a radius of up to 50 km and an area of up to ~7850 km². The centre of cyclonic eddy e3 was located at 68.125°E, 16.875°N and had a radius of up to 41 km and an area of up to ~5280 km². Each eddy radius is defined as the average distance between the eddy boundary and the centre of the eddy identified by the eddy detection algorithms. The long forcing time of strong wind during the passage of Kyarr and Maha are shown in Table 2 . The maximum forcing time of each eddy was all much longer than the geostrophic adjustment time during the passage of Kyarr. However, during the passage of Maha, the maximum forcing time was all slightly below the maximum adjustment time in these three cold eddy regions.

Figures 2A–E show the variations in temperature profiles from five Argo floats. Figure 2F shows the temporal and spatial variations in the five selected Argo floats. According to Figure 2F , Argo floats were located near the tracks of Kyarr and Maha. And Argo floats 2902202, 2902205, and 2902272 were concentrated near e1. Argo float 2902175 was somewhat far from e1, approximately 100 km from the edge of e1, and Argo float 2902210 was located approximately 50 km from the edge of e2. The five Argo floats recorded changes in the ocean at different times during the passage of Kyarr and Maha. For example, Argo floats 29002210 and 2902272 recorded changes in the interior of the ocean caused by TC Kyarr, and Argo floats 2902205, 2902202, and 2902175 recorded changes in the interior of the ocean caused by TC Maha.

Furthermore, similar to the Argo float observations, the spatial series of the HYCOM model output results were added to the temperature profile to demonstrate that HYCOM model can reproduce ocean conditions. However, the upper ocean temperature changes captured by the HYCOM model output after the passage of TC Kyarr were more affected than the temperature changes recorded by Argo floats, and the temperature decreased even more dramatically ( Figures 2A–E ). In general, the HYCOM model output results were generally consistent with the Argo float measurements.

In this study, the upper ocean is defined as depths of 0–200 m. To investigate the thermal response of the upper ocean, the heat content change (ΔH ) was calculated from a pair of temperature profiles from each Argo float (Zedler et al., 2002).

Here, ρ0 is the density of seawater, and Cpw is the specific heat of seawater, so ρ0Cpw = 4.1 MJ°C^-1 m^-3. Integration is performed over depths of 0 to 200 m. ΔT is the change in sea water temperature before and after the TCs, as shown in Figure 2 for the Argo float temperature profile of the time series. A negative value indicates heat loss and a positive value indicates heat gain within a specified depth range. As shown in Figure 9 , the ΔH calculated from the Argo profiles demonstrated a significant difference.

Argo floats 2902175 and 2902202 showed that the upper 200 m of the ocean lost heat along the temperature profiles, while Argo floats 2902210, 2902205, and 2902272 showed that heat was gained in the subsurface ocean (50-100 m). And Argo float 2902210 showed that a subsurface ocean heat increase could reach 13 MJ°C^-1 m^-3, Argo float 2902205 showed a subsurface ocean heat increase of approximately 6 MJ°C^-1 m^-3. Argo float 2902272 showed a subsurface ocean heat increase of 2 MJ°C^-1 m^-3. Moreover, Argo float 2902210 showed an ocean heat gain of approximately 180 m. In previous sections, it was found that the passage of Kyarr led to the generation of three cold eddies. The interior of the cold eddy is an unstable thermal structure. Therefore, the heat flux exchanges inside cold eddies need to be further studied.

Some of the nondimensional numbers defined in previous studies help to reveal the characteristics of TCs and their interactions with the upper ocean. The external parameters that characterize the ocean environments and the two sequential TCs are exhibited and listed in Table 3 . The calculation results for each TC are presented in Table 4 . For each TC, the nondimensional storm speed S is O(1), where S is the ratio of the local inertial period to the TC residence time (Price et al., 1994). The strong wind stresses of TCs caused changes in the interior of the ocean that were basically consistent with the local inertial period. The Burger number B reveals the degree of pressure coupling between the mixed-layer currents and the thermocline currents (Price et al., 1994). For each TC case, a relatively large Burger number can be seen, and it is expected that the pressure coupling is very pronounced during the TC passage. The translation speeds of TCs Kyarr and Maha were only slightly greater than c (the gravest mode internal wave). Each TC had an O(1) Mach number. Previous studies have shown that under these conditions, significant upwelling occurs directly beneath TCs (Price et al., 1994). And it was caused by Ekman pumping, which is closely linked to positive wind stress curls ( Figure 5 ). Therefore, the calculations of these nondimensional numbers for each TC case show that two slow and strong TCs have considerable impacts on the upper ocean interior during their passage (Price et al., 1994).

Figures 10A–C show the changes in vertical temperature profiles and the depth of mixed-layers before, during and after Kyarr and Maha in the upper ocean at the three cyclonic eddy centres. At each eddy centre, the initial MLD was approximately 40 m, and the MLT was approximately 29°C. Then, the MLT and MLD changed sharply due to the influence of the passage of TCs. The cooling of the MLT was accompanied by the deepening of the MLD; the MLT decreased by at least 3°C, and the MLD deepened to at least 55 m at each eddy centre after Kyarr passed. For example, after the passage of Kyarr, the MLT in the central area of e2 decreased the most sharply, by approximately 5°C ( Figure 10B ), the MLD deepened by approximately 75 m. After the passage of Kyarr in e3, the MLD deepened the most, up to 80 m, and the MLT dropped to approximately 26°C ( Figure 10C ). After Maha passed through the three cold eddies, it also caused the MLD in each eddy centre to deepen and the MLT to decrease, but none of the changes were as pronounced as those after Kyarr passed through. The differences in the variations in MLD and MLT at the e1, e2, and e3 centres were mainly due to the different characteristics of each TC when they were passing these cold eddy fields. Based on the nondimensional characteristics of each TC calculated in section 4.8, two slow-moving TCs can cause strong vertical mixing and upwelling, causing the cooling of the MLT and deepening the MLD. What’s more, the sea water temperature at 150-200 m remained virtually unchanged during the whole period.

The passage of TCs greatly influenced the depth integrated heat (DIH) content of the upper ocean. Figures 10A–C show the temperature profiles of each eddy centre. The maximum temperature in each eddy centre response to each TC occurred approximately 48 h after the TCs passed. Therefore, T0(z) was the 24-h average during the temperature response to TCs at each eddy centre before the TCs (on October 24). And T1(z) was the 48-h average during the temperature response to TCs at each eddy centre after the TCs. The DIH anomaly could be obtained from the temperature profiles by integrating the temperature anomaly (T1(z)–T0(z)) over an appropriate depth range and multiplying it by the density and specific heat of sea water. Thus, following Zedler (Zedler, 2002), the computation of the DIH anomaly is as follows:

Integration was performed over depths of z1 to z2 (z is positive downwards). To examine the separate contributions between different layers of the upper ocean, three depth ranges were chosen for integration: a shallow depth range as the MLD, a middle layer that was the MLD bottom to the 150 m depth, and a bottom layer that was the 150 m to 200 m depth. The shallow range effectively quantified the MLD cooling occurring after the passage of TCs and was therefore denoted DIH_ML; the middle layer was denoted DIH_MD, and the bottom layer was denoted DIH₂₀₀.

The depth ranges and the corresponding DIH values of each eddy centre for the TCs are summarized in Table 5 . The two TCs were associated with heat loss from the three different layers, but the largest loss occurred in the middle layer after Maha, where DIH_ML = -1220.0 MJ m^-2 in the e1 eddy centre, and another large middle layer heat loss in the Maha case similar to that observed for the e1 eddy centre occurred when DIH_ML = -1192.3 MJ m^-2 over the depth range [53 m, 150 m]. These different DIH values may be partly due to the difference in vertical mixing and upwelling intensity in the upper ocean caused by TCs Kyarr and maha as they passed through the central region of each cold eddy.