Increasing trend in rapid intensification magnitude of tropical cyclones over the western North Pacific

Rapid intensification (RI) is defined to be a significant increase in tropical cyclone (TC) intensity over a short time period and is particularly difficult to simulate and predict. Consequently RI continues to pose a significant challenge for operational forecasting (Knaff et al 2018). Due to continued concerns about climate change and its impacts on TC activity, there has been an increasing focus on temporal variations related to RI activity. In particular, several studies have examined RI over the western North Pacific (WNP), where 90% of Category 4 and 5 TCs on the Saffir-Simpson Hurricane Wind Scale [e.g. 1 min maximum sustained wind (V_max) ≥ 113 kt] experience at least one RI episode, defined to be a V_max increase of ≥30 kt in 24 h, during their lifetimes. The WNP also has the highest ratio of TCs with RI to the total number of TCs of any global TC basin, on average (Wang and Zhou 2008). Super Typhoon Haiyan in 2013, which attained a peak 1 min maximum sustained wind intensity of 170 kt per Joint Typhoon Warning Center (JTWC) estimates (Chu et al 2002), is a recent example of how a TC undergoing RI (RI-TC) could become more intense in the context of global warming (Lin et al 2014).

Changes in the frequency of RI over the WNP have been discussed in several previous publications (e.g. Wang and Zhou 2008, Wang et al 2015, Wang and Liu 2016, Fudeyasu et al 2018, Kang and Elsner 2019). On interannual timescales, there is a significant relationship between RI and El Niño-Southern Oscillation (ENSO), with not only a larger average number of RI occurrences but also a greater proportion of RI-TCs to the total number of TCs in El Niño years than in La Niña years (Wang and Zhou 2008, Fudeyasu et al 2018). Wang and Zhou (2008) noted that during El Niño years, WNP TC occurrences migrated equatorward and eastward on average, where there was increased upper ocean heat content favoring RI. They also found that weakened vertical wind shear (VWS), in response to a Rossby wave caused by enhanced heating over the central equatorial Pacific, favored more RI occurrences during El Niño years. Fudeyasu et al (2018) showed that more TCs forming over the southeastern quadrant of the WNP during El Niño years increased the opportunity for these TCs to pass through the region where WNP RI typically occurred, thereby increasing the rate of RI-TCs.

There is also an apparent decadal variation in the frequency of WNP RI that was noted in Wang and Zhou (2008) and investigated in detail by Wang et al (2015). In general, more and fewer RI events occur during negative and positive phases of the Pacific decadal oscillation (PDO), respectively (Wang et al 2015). Over the main RI region, thermodynamic effects dominate decadal changes in RI occurrence over the main RI region, through increasing (decreasing) TC heat potential (TCHP) during the negative (positive) PDO phase (Wang et al 2015). During the negative PDO phase, the extremely warm water of the equatorial Pacific is advected into the main RI region by the steering of an anomalous low-level anticyclonic circulation in the subtropical gyre of the WNP.

In addition, Wang and Zhou (2008) and Kang and Elsner (2019) found no significant trend in the annual RI-TC number over the WNP during 1965–2004 and 1986–2015, respectively. However, due to the decreasing tendency in total TC occurrence, there is a significant increasing proportion of RI-TCs since the 1980s (Fudeyasu et al 2018, Kang and Elsner 2019). Kang and Elsner (2019) further suggested that global warming explained more than half of the change in the ratio of RI-TCs. Similar to what was found for RI-TCs, Wang et al (2015) and Kang and Elsner (2019) found no trend in annual RI frequency, whereas the annual ratio of RI cases to total cases exhibited a significant increasing trend (Bhatia et al 2019).

Most of the aforementioned publications have focused on the variation in the frequency of either RI events or RI-TCs over the WNP, as well as the large-scale conditions driving these changes. Fewer studies have examined the variability of WNP RI magnitude over the past several decades. Recently, several studies have noted that intensification rates of TCs have increased since the 1980s over the globe as well as in some individual TC basins (Kishtawal et al 2012, Bhatia et al 2019). Using quantile regressions for 24 h intensity changes (ΔV₂₄s) for all TC cases during 1982–2009, Bhatia et al (2019) reported that the 95th percentile exhibited the largest trends. These trends were approximately 4 kt decade⁻¹ and 3 kt decade⁻¹ over the globe and in the Atlantic, respectively. Since the 95th percentile of ΔV₂₄s happens to be the traditional threshold for identifying RI (Kaplan and DeMaria 2003, Kaplan et al 2010, Shu et al 2012, Knaff et al 2018), these positive linear trends indicate that the Atlantic and global RI magnitudes have increased during the period from 1982 to 2009. However, it is still unknown how RI magnitude has changed over the WNP. Therefore, the present study investigates the potential long-term trend in the magnitude of WNP RI and its possible relationship with changes in the large-scale environment.

The remainder of this paper is organized as follows. Section 2 describes the data used and the methodology employed to identify different categories of RI events. Section 3 discusses the spatiotemporal variations in different RI categories as well as the environmental factors driving these variations. The paper concludes with a summary in section 4.

The original 6-hourly WNP TC best track data from the JTWC during 1979–2018 used in this study are provided by the International Best Track Archive for Climate Stewardship (IBTrACS) (v04r00; Knapp et al 2010), including TC central position and V_max. To reduce the uncertainty in detecting weak TCs such as tropical depressions (Klotzbach and Landsea 2015), we only consider TCs with lifetime peak intensities of at least 34 kt. Similar to previous studies (Kaplan and DeMaria 2003, Kaplan et al 2010, Shu et al 2012, Knaff et al 2018, Shimada et al 2020), an RI event is defined as ΔV₂₄ of at least 30 kt for over-water TCs that were tropical in nature, which is approximately the 95th percentile of all ΔV₂₄s for all TC cases. Note that RI events defined here are counted in 6 h intervals when a 24 h intensification rate of at least 30 kt occurs. Consequently, a TC may undergo more than one RI event during its lifespan. Figure 1 displays the annual numbers for all TCs, RI-TCs and RI events during 1979–2018. On average, around 27 TCs form over the WNP each year, with ∼12 TCs experiencing at least one RI event during their lifetimes. Considering that, on average, 54 RI events occur per year, there are ∼5 RI events for the average RI-TC. Additionally, there is a significant increasing trend in the number of RI events per RI-TC during 1979–2018, with a rate of 0.02 events RI-TC⁻¹ yr⁻¹ (p = 0.03).

Zoom In Zoom Out Reset image size

Shimada et al (2020) reported that there was an increase in the number of WNP RI-TCs from 1987 to 2018, which was likely caused by the temporal inhomogeneity of the intensity-estimating technology applied in the Japan Meteorological Agency (JMA) best track data. In order to gauge the impact of this inhomogeneity on temporal variations related to RI, we compare the annual number of RI-TCs and RI events in the JTWC data with 3-hourly Advanced Dvorak Technique-Hurricane Satellite-B1 data (ADT-HURSAT; Kossin et al 2013) during 1982–2009. The ADT-HURSAT data are more homogeneous and have been widely applied as a reference for trends related to TC intensity (Kossin et al 2013, 2014, Bhatia et al 2019). Figure 1 displays the annual numbers of all TCs, RI-TCs and RI events for JTWC. All three of these indices as calculated from JTWC data significantly correlate with the same timeseries from ADT-HURSAT, illustrating the reliability of RI events defined through JTWC data. Additionally, consistent with the trends calculated using the Dvorak-based RI events in Shimada et al (2020), the long-term tendencies are 0.03 and 0.41 yr⁻¹ for the annual numbers of RI-TCs and RI events, respectively, which are not significant at the 0.05 level based on a Student's t-test.

In addition to the definition of RI events, a threshold of 50 kt is applied to categorize RI events into different groups, which is the 99th percentile of all ΔV₂₄s for all WNP TC records. Individual RI events are then classified as weak-moderate RI (30 kt ≤ ΔV₂₄ < 50 kt) and strong RI (ΔV₂₄ ≥ 50 kt). The magnitude of RI directly refers to ΔV₂₄ for individual RI events (Balaguru et al 2018).

Favorable large-scale environmental factors are critical for the occurrence of RI, including warm sea surface temperatures (SSTs), high ocean heat content (OHC), a moist low-middle troposphere and low vertical wind shear (VWS) (Shu et al 2012, Fudeyasu et al 2018). Two atmospheric environmental factors, 850–200 hPa VWS and 600 hPa relative humidity (RH), are calculated using monthly mean data from the European Centre for Medium-Range Weather Forecasts (ECMWF) Interim Re-Analysis (ERA-Interim) with a horizontal resolution of 1.5° latitude × 1.5° longitude (Dee et al 2011). Monthly mean SST data over a 1° × 1° grid are obtained from the Hadley Centre Sea Ice and Sea Surface Temperature dataset (HadISST; Rayner et al 2003). TCHP is a measure of OHC that is warmer than 26 °C (DeMaria et al 2005) and is calculated using monthly subsurface temperature profiles from the control member in the ECMWF Ocean Reanalysis System 5 (ORAS5; Zuo et al 2019) with a resolution of 1° latitude × 1° longitude. HadISST and ERA-Interim provide the primary forcing fields for ORAS5 since 1979 (Zuo et al 2019). In addition, the monthly Pacific Decadal Oscillation (PDO) index, represented by the leading PC of monthly SST anomalies from HadISST over the North Pacific Ocean, is obtained from the NOAA Earth System Research Laboratory's Physical Sciences Division (PSD) from 1948 to 2018. The ECMWF Ocean Reanalysis System 4 (ORAS4; Balmaseda et al 2013) data over a 1° × 1° grid are applied to compute TCHP over a longer time span (1958–2017).

The significance levels (p) of correlation coefficients (r), linear trends and the differences in means between two samples are all estimated using a two-tailed Student's t test.

3.1. Changes in RI magnitude and RI numbers for different categories

Figure 2(a) displays linear trends in individual quantiles of ΔV₂₄ over the WNP from 1979 to 2018, which are derived by the least squares for individual quantiles of ΔV₂₄ as a function of year. TC best track data during this period, which are primarily estimated from satellite observations, are more reliable for a climatological analysis than data in earlier years (Chu et al 2002). Higher and lower quantiles of ΔV₂₄ generally yield positive and negative trends, respectively, with the slope typically becoming greater with increasing quantile. Unlike the Atlantic and global ΔV₂₄ cases (Bhatia et al 2019), the downward trends in the lower quantiles are not all significant over the WNP, suggesting that the weakening rate of WNP TCs has generally remained unchanged. There are significant upward trends in the 75%, 90% and 95% quantiles, however, highlighting the increasing frequency of TC cases with larger ΔV₂₄. Similar to what was found in Bhatia et al (2019), the 95% quantile shows the largest trend at 0.23 kt yr⁻¹, although the slope is relatively small compared with what Bhatia et al (2019) found for the globe and for the Atlantic.

Zoom In Zoom Out Reset image size

Because the 95% quantile of ΔV₂₄ approximately corresponds to the threshold that has historically been used for defining RI (Kaplan and DeMaria 2003, Kaplan et al 2010, Shu et al 2012, Knaff et al 2018), a positive slope in this quantity indicates an increasing trend in the intensification magnitude of TCs undergoing RI in the WNP. During 1979–2018, this rate has increased by 0.10 kt yr⁻¹ (p< 0.01) as displayed in figure 2(b). Note that the mean RI magnitude (defined as the 95% quantile of ΔV₂₄) is 39.0 kt over the WNP, which is greater than that over the Atlantic (29.8 kt reported in Balaguru et al 2018, where RI referred to ΔV₂₄ with at least 25 kt). The annual RI frequency does not significantly relate to RI magnitude over the period from 1979 to 2018, with a correlation between these two parameters of 0.19 (p= 0.25). The proportion of RI cases to the total number of TCs also does not significantly correlate (r= 0.30, p= 0.06). This means that the variation of RI magnitude is not directly related on annual timescales to RI frequency over the WNP.

We next categorize all RI events into two groups and find that the annual frequency of strong RI events from 1979 to 2018 exhibits somewhat different trends from that of weak-moderate RI (figure 2(c)). There is not a significant relationship between the annual numbers of strong RI events and weak-moderate RI events (r= 0.23, p= 0.15). Weak-moderate RI frequency slightly increases at a rate of 0.16 events yr⁻¹ (p= 0.39), which leads to an insignificant trend in the total RI number as also reported in Wang et al (2015) and Kang and Elsner (2019), due to the large proportion of weak-moderate RI cases to the total number of RI cases (∼82%). However, there is a significant increasing tendency in the frequency of strong RI, with a slope of 0.25 events yr⁻¹ (p < 0.01). These distinct trends in annual numbers of different RI categories induce a significant increasing ratio of strong RI cases to the total number of RI cases (0.32% yr⁻¹, p = 0.01), further leading to the RI distribution shifting to larger ΔV₂₄.

To highlight the impact of the change in the strong RI number, the mean RI magnitude (${\overline {\Delta V} _{{\text{RI}}}}$) in each year is first decomposed as:

$$\begin{equation}{\overline {\Delta V} _{{\text{RI}}}} = \frac{{{n_{{\text{strong}}}} \times {{\overline {\Delta V} }_{{\text{strong}}}} + {n_{{\text{weak}} - {\text{moderate}}}} \times {{\overline {\Delta V} }_{{\text{weak}} - {\text{moderate}}}}}}{{{n_{{\text{strong}}}} + {n_{{\text{weak}} - {\text{moderate}}}}}}\end{equation} \tag{ 1 }$$

where ${\overline {\Delta V} _{{\text{strong}}}}$ (${\overline {\Delta V} _{{\text{weak}} - {\text{moderate}}}}$) denotes the average magnitude of strong (weak-moderate) RI events, while ${n_{{\text{strong}}}}$ (${n_{{\text{weak}} - {\text{moderate}}}}$) refers to the number of strong (weak-moderate) RI events. We then assess the individual importance of these four variables in determining the trend of ${\overline {\Delta V} _{{\text{RI}}}}$. As in Camargo et al (2007), we recalculate ${\overline {\Delta V} _{{\text{RI}}}}$ using the long-term means for three of the variables and annually-varying values for the fourth variable. This procedure is repeated for each of the other three variables. Figure 3 shows that the change in the strong RI number is the primary contributor to the increase in RI magnitude during 1979–2018, causing a significant increasing trend of 0.08 kt yr⁻¹ (p < 0.01) when the three other variables remain unchanged. In contrast, there are no significant trends in RI magnitude if the influence of the strong RI number increase is removed. Along with the strong correlation between RI magnitude and strong RI number (r = 0.83, p < 0.01), we conclude that the increasing number of strong RI events has primarily driven the increasing magnitude of total RI over the WNP.

Zoom In Zoom Out Reset image size

Figure 4 displays the spatial distributions of mean RI occurrences during 1979–2018 in different categories, as well as their differences between 1979–1998 and 1999–2018. Similar patterns in occurrence numbers are shown for the total number of RI events, as well as when these events are sub-divided into weak-moderate and strong RI events. The most frequent occurrences for all RI categories is east of the Philippines (figures 4(a), (c) and (e)), which is consistent with the main RI region defined in previous studies (Wang and Zhou 2008, Shu et al 2012, Wang et al 2015). However, there are substantial differences in the occurrence changes in different RI categories. For weak-moderate RI events, the occurrence differences between 1979–1998 and 1999–2018 display a northwest-southeast dipole pattern, with the maximum increase and decrease occurring around 20°N, 135°E and 15°N, 145°E, respectively (figure 4(d)). A similar feature is also seen in the occurrence change of the total number of RI events (figure 4(b)), since weak-moderate RI events comprise the large majority of the total number of RI events (∼82%). Note that the two sub-periods considered here correspond to different PDO phases (Wang et al 2015). On average, weak-moderate and total RI occurrences have migrated poleward and westward since 1979, associated with the PDO phase shift from positive (∼1979–1998) to negative (∼1999–2018). This feature is also reported by Wang et al (2015), who found that the western part of the WNP generally had more favorable large-scale oceanic and atmospheric conditions in negative PDO phases than in positive PDO phases.

Zoom In Zoom Out Reset image size

Substantial increases in strong RI occurrences are found over most of the WNP (figure 4(f)). The maximum in this increase is located around 15°N, 135°E, which corresponds to the region with the highest occurrence of strong RI events in the long-term climatology (figure 4(e)). Strong RI occurrences have decreased slightly east of 155°E. Consequently, the occurrence difference pattern of strong RI events is different from that of weak-moderate and total RI events, with strong RI events having a west-east dipole structure. This difference pattern suggests that there are potential differences in the mechanisms of large-scale environmental conditions influencing the changes in strong and weak-moderate RI occurrences over the WNP.

3.2. Environmental conditions affecting strong RI occurrences

In order to identify possible environmental conditions responsible for the variations in strong RI occurrences, we investigate changes in four oceanic and atmospheric variables: SST, TCHP, 600 hPa RH and 850–200 hPa VWS. These factors are similar to those examined in Wang et al (2015). Instead of examining the whole WNP, we investigate two sub-regions that have exhibited different changes in strong RI occurrences: Region A (10°∼20°N, 120°∼150°E) where there has been a substantial increase in strong RI and Region B (12.5°∼22.5°N, 155°∼170°E) where there has been a slight decrease in strong RI. There has been significantly greater SSTs over almost the entire WNP in 1999–2018 than in 1979–1998 (figure 5(a)), which has been primarily linked to anthropogenic global warming since 1979 (Chan and Wu 2015). As would be expected from the PDO phase change from positive to negative, the strongest SST increases are located over the subtropical WNP. Moreover, there are higher (lower) TCHPs west (east) of around 170°E over the tropical WNP in 1999–2018 than in 1979–1998 (figure 5(b)). The maximum increases in TCHP are centered over Region A, with a significant increase of ∼35 kJ cm⁻² from 1979–1998 to 1999–2018, on average. There are also significant increases in TCHP over Region B, although these increases are of a smaller magnitude (∼21 kJ cm⁻²). Changes in oceanic conditions (SST and TCHP) favor an increase in strong RI occurrences over both Regions A and B.

Zoom In Zoom Out Reset image size

Differences in 600 hPa RH between 1979–1998 and 1999–2018 are displayed in figure 5(c). The maximum increases and decreases have occurred over the equatorial WNP and the East Asian continent, respectively. The differences in 600 hPa RH averaged over Regions A and B are increases of only 1.8% and 1.7%, respectively, which are not statistically significant. This implies that strong RI occurrences are not sensitive to small changes in mid-level moisture.

Unlike other large-scale variables, 850–200 hPa VWS differences exhibit opposite signs over Regions A and B (figure 5(d)). Over Region A, the change in VWS during 1979–2018 is not significant, with a small reduction of −0.1 m s⁻¹ from 1979–1998 to 1999–2018. This indicates that variations in atmospheric conditions have had little impact on the change in strong RI occurrences over Region A, while the increasing strong RI frequency has primarily been a result of the warming ocean. In contrast, there has been a significant increase of 1.4 m s⁻¹ in VWS from 1979–1998 to 1999–2018 over Region B, which likely suppresses the occurrence of strong RI events in this region. The slight decrease in strong RI occurrences over Region B are likely linked to the RI-suppressing influence of increased VWS offsetting the RI-favoring influence of the warming ocean.

These same environmental conditions induce distinct patterns in the occurrence differences in strong and weak-moderate RI events. The more weak-moderate RI events north of 20°N correspond to significantly greater SST and larger 600 hPa RH (figures 5(a) and (c)), indicating that warmer SSTs and a moist mid-level atmosphere both favor the development of weak-moderate RI. In the same region, although TCHP and VWS exhibit a slight increase and decrease, respectively, their changes are not statistically significant (figures 5(b) and (d)). Large-scale environmental factors modulating weak-moderate RI occurrences are different from those influencing strong RI occurrences. Note that the region with significantly reduced weak-moderate RI events south of 20°N happens to be the area with increased strong RI occurrences. This increase is likely caused by the increased RI magnitude that turns some weak-moderate RI events into strong RI events.

We have shown that increases in TCHP over Region A are the predominant cause of the increase in strong RI occurrences over the WNP. On decadal timescales, there is a significant positive correlation (r= 0.94, p= 0.01) between 11 year running averages of the annual strong RI number and TCHP during 1979–2018 (figure 6). There is a significant increasing trend in TCHP in Region A from 1979 to 2018, with a rate of 0.61 kJ cm⁻² yr⁻¹ (p< 0.01), likely fueling the increasing frequency of strong RI events. Since the PDO phase shifts from positive to negative during the period investigated here, there is a significant downward tendency in the PDO index (−0.03 standard deviations yr⁻¹, p= 0.04) during 1979–2018. Further analysis shows that the long-term trends in the annual strong RI number and TCHP primarily result from the significant difference in the mean states between 1979–1998 and 1999–2018. The mean strong RI number and TCHP in 1979–1998 (positive PDO phase) are 7 and 70.7 kJ cm⁻², respectively, which are significantly smaller than the 12 and 86.0 kJ cm⁻² values observed in 1999–2018 (negative PDO phase). There are no significant trends in the number of strong RI events during 1979–1998 (0.30 yr⁻¹, p= 0.17) and during 1999–2018 (0.27 yr⁻¹, p= 0.20). The trends in TCHP are also not significant during 1979–1998 (0.30 kJ cm⁻² yr⁻¹, p= 0.51) and during 1999–2018 (0.03 kJ cm⁻² yr⁻¹, p = 0.95). In addition, decadal variations in strong RI occurrences and in TCHP are inversely correlated with the PDO index between 1979 and 2018.

Zoom In Zoom Out Reset image size

Given the short record examined here, it is difficult to determine whether the increasing trend in strong RI occurrences is induced by the transition between PDO phases, anthropogenic global warming or a combination of both factors. One method for overcoming this difficulty is to analyze a longer period including several PDO phases. However, TC intensity estimates are less reliable before the 1970s. These intensity estimates were often overestimates in the WNP prior to the use of meteorological satellites (Chu et al 2002). We examine decadal variations in TCHP (derived from the ORAS4 data) over Region A compared with the PDO since the late 1950s (figure 6). There is a persistent increasing trend in TCHP over Region A regardless of PDO phase changes, possibly linked to anthropogenic global warming. Consequently, we hypothesize that the number of strong RI events may have also increased over a longer time period.

In this study, the long-term trend (e.g. from 1979 to 2018) in RI magnitude of TCs over the WNP and the environmental contributors to the observed RI changes are investigated. A TC may undergo more than one RI event during its lifespan, and there are ∼5 RI events for the average RI-TC. There is a significant increasing trend of 0.10 kt yr⁻¹ (p< 0.01) in RI magnitude of WNP TCs, which is also demonstrated in significant positive tendencies in higher quantiles of ΔV₂₄. These variations in RI magnitude also evidence themselves in changes in the number of strong RI episodes which have increased by 0.25 events yr⁻¹ (p< 0.01) from 1979 to 2018. In contrast, changes in weak-moderate RI events show a slight increasing (but insignificant) trend and are also insignificantly correlated with variations in RI magnitude. Because of the significant increasing trend in strong RI numbers and the insignificant trend in weak-moderate RI numbers, there is an increasing proportion of strong RI events to total RI events, which induces the increasing RI magnitude of WNP TCs.

When the period is divided into two halves (1979–1998 and 1999–2018), there are distinctly different spatial patterns in the occurrence differences of weak-moderate and strong RI events. For weak-moderate RI events, there is a northwest-southeast dipole pattern in the occurrence difference, which means that the average position for RI events migrated poleward and westward from 1979–1998 to 1999–2018. As reported in Wang et al (2015), this shift is likely linked to the transition of the PDO phase from positive to negative. In contrast, the difference in strong RI occurrences exhibits a west-east dipole structure, with substantial increases and slight decreases in 1999–2018 compared with 1979–1998, occurring west and east of around 155°E, respectively. After analyzing differences in large-scale environmental conditions, we find that the significantly increasing strong RI occurrences over Region A (10°∼20°N, 120°∼150°E) are primarily induced by the warming ocean (increasing SST and TCHP), whereas changes in atmospheric variables (RH and VWS) play only a minor role. The slight decrease in strong RI events over Region B (12.5°∼22.5°N, 155°∼170°E) is likely linked to the RI-unfavorable influence of stronger 850–200 hPa VWS partially offsetting the RI-favorable influence of the increasing SST and TCHP.

Due to the relatively short period with reliable TC intensity estimates [e.g. the satellite era (since 1979)], it is difficult to determine whether the increasing RI magnitude of WNP TCs is caused by the transition between PDO phases, anthropogenic global warming, or a combination of these factors. By comparing variations in TCHP and the PDO since the late 1950s, TCHP has been shown to persistently increase regardless of PDO phase. Given the dominant role played by the warming ocean in the increasing number of strong RI events, we hypothesize that the increasing trend in strong RI frequency and RI magnitude are likely primarily fueled by global warming.

We would like to express our sincere thanks to three anonymous reviewers for their helpful comments on an earlier manuscript. This work was jointly funded by the National Natural Science Foundation of China (61827901) and the National Key Research and Development Program of China (2018YFC1507305). Klotzbach would like to acknowledge financial support from the G. Unger Vetlesen Foundation.

All data used in this study are freely available online. The data that support the findings of this study are available from the corresponding author upon reasonable request.