The Xiadonghai embayment, the area for the present study, is near Luhuitou on the south coast of Hainan Island (Figure 1). Tides here are irregularly diurnal, with a mean tidal range of 0.85 m (State Oceanic Administration, 1999) and a weak tidal current, i.e., generally 0.1–0.2 m/s. The local wave climate is seasonal and highly influenced by typhoon activity with around 80% of waves approaching the coast from East South East (ESE) to West South West (WSW) directions. The average wave height range is 0.6-0.8 m, and the wave period range is 2.5–3.5 s (State Oceanic Administration, 1999). Storminess is a feature of the wet season from March to November. In general, ~five typhoons (at least three being very strong) hit this coast per year. A maximum summer wind velocity of 45 m/s was recorded during Typhoon Kelly (Wang et al., 1998). Typhoon Elaine on October 9, 1971, produced a maximum storm surge, reaching 3.92 m above chart datum (Wang et al., 1998). Reliable measurements and documentation for typhoon wave data are only back to the 1980s in the study area. By extrapolating the typhoon wave data for Hainan Island (22 years of observations), maximum significant typhoon wave heights are calculated as ~4.84 m every 2 years, 9.51 m every 50 years, 10.80 m every 100 years (Yao et al., 1991), and 18.7 m every 200 years (Yin, 2014).

Xiaodonghai (XDH) belongs to the Sanya national coral reef reserve. The reef flat is ~300 m wide and >1,800 m long. A back-reef moat of 1 m average depth lies between the reef crest and beach. The deepest parts of the moat reach 1–2 m but become shallow (70 cm) during low tide when the reef crest and much of the moat (trough) dry up (Zhang, 2001).

In situ investigations at the XDH study site were carried out in December 2016 during low spring tides. Most of the boulders are apparently derived from the reef front, as generally shown by their overlapping erosional marine features, barnacles on their upper surface, and other faces with a freshly cut appearance. For the boulders on the reef flat with a long axis length of >0.5 m, we measured the dimensions of long (a), intermediate (b), and short (c) axes, and the long axis directions. The positions and elevations of CCBs were measured using an RTK–GPS survey System, at a vertical and horizontal accuracy of ±1.5 and ±1 cm, respectively. In total, 1247 CBDs were measured in the field.

Small pieces were collected from 42 typical boulders for the estimation of bulk density, carried out at the laboratory using the following equation:

where ρ_b is the boulder bulk density, ρ_w is the seawater density (1.02 g/cm³), W_a is the boulder weight in the air, and W_f is the boulder submerged weight in seawater.

As CBDs are often irregular and non-rectangular in shape, their volumes (V) were estimated using the following equations to avoid over-estimation (Lau et al., 2016):

For the larger boulders (a>2 m), we also use a Riegl VZ-4000 terrestrial laser scanner to measure the boulder volume which was calculated by a direct point-to-point cloud comparison method (Xie et al., 2017).

Boulder mass was estimated as the product of volume and bulk density. The bulk density was carried out on 10 selected coral samples and calculated as 1.82 t m⁻³ based on Equation (1).

The morphology of sedimentary particles [e.g., maximum projection sphericity (MPS), oblate prolate index (OPI), disc–rod index (DRI), and roundness] is widely recognized as an important source of information on sediment provenance, transport, and dynamic environment (Table 1; Dobkins and Folk, 1970; Pettijohn et al., 1987; Illenberger, 1991; Wang et al., 2004; Chen et al., 2011). OPI is a useful index for differentiating oblate from prolate boulder shapes. MPS reflects the balance of drag and gravitational forces acting on a particle immersed in a fluid and, as such, is an important index of the behavior of particles transported by or settling in the water. DRI is effective in distinguishing rod-shaped from disc-shaped boulders.

Coral boulders on the Yulin Bay reef flat originated from two possible sources (Figures 2, 3, 5): (1) the active reef edge or reef slope (collectively called “reef-edge boulders”) and (2) from mid-Holocene emerged reef remnants or other emerged features (collectively called “emerged reef boulders”). When a boulder exhibited a surface appearance resembling that of the remnant reef at the locality and had at least one axis that was shorter than the elevation of the emerged reef (~1.1 m), the boulder was considered to originate from the emerged reef (i.e., “emerged reef boulders”).

To characterize the current velocity of moving boulders along the XDH site, the minimum current velocity (MCV) needed to dislodge and transport boulders were calculated using the following equations (Nandasena et al., 2011).

(sliding)

(rolling)

(lifting) For a reef-edge boulder in a joint-bounded setting, transport can only be initiated by lifting when u is the current velocity (m/s), C₁ is the lift coefficient (0.178), θ is the angle of bed slope, and μ_s is the coefficient of static friction (0.7).

The U/Th age dating of the youngest surfaces of large boulders helped to determine the deposition time. The U/Th dating analyses were conducted at a multi-collector inductively coupled plasma mass spectrometer (MC-ICPMS) follow a standard procedure performed by Edwards et al. (1987) and Cheng et al. (2013). U/Th ages were calculated using standard equations. We assume that the U/Th age date of the outermost edge of a dead coral head is the time the coral head died (Terry et al., 2016). It is assumed that growing corals were separated from the living reef surface when a boulder was transported from the submerged reef slope and deposited onto the reef flat by high-energy wave events (Yu et al., 2009; Lau et al., 2016). This assumption is supported by previous studies at various sites around the SCS (Zhao et al., 1983; Zhang, 2001; Yu et al., 2012), because “live/dead corals” have been observed presented on the reef slope and fresh coral fragments found were on the reef flat, suggesting living corals have been distributed on the reef slope. The assumption has also been widely validated by the coral dead ages in the top surface of large boulders combining with known historical extreme wave events (e.g., Yu et al., 2012; Lau et al., 2016, 2018). In the study, eight coral samples were collected based on the method of Terry et al. (2016) from the largest 12 boulders for U/Th age dating (Supplementary Figure 1). We also collected one oyster shell attached to larger boulders that were above the maximum tide line (Supplementary Figure 1) for AMS ¹⁴C dating and calibrated using IntCal13 calibration curve (Reimer et al., 2013) with are regional Delt-R of 18 ± 37 (Southon et al., 2002). Therefore, our radiocarbon ages represent minimum dates for boulder dislocation. Our age dating results are listed in Table 2.

Figure 2 shows the positions and sizes of measured CBDs at the XDH study site. The fringing reef slopes seaward at an angle of ~2° (Figure 1) and consists of dying, primarily, Porites type corals. The slightly emerged reef flat is covered by numerous dead and often inverted Porites CBDs with long axes of 0.5–4.0 m (Figures 2, 3). Boulder shapes were platy, sub-angular to rounded, or rectangular to ellipsoidal with some having sharp broken edges (Figures 2b, 3b). The average boulder volume and weight were 0.43 m³ and 0.78 t, respectively (Table 3). The largest and smallest boulder a-axis dimensions were 3.76 and 0.52 m, respectively, while 2.86 and 0.31 m for b-axes, and 2.51 and 0.14 m for c-axes. The estimated range in mass was 0.05– 28.81 t (Table 3). The volume of the largest boulder was ~16 m³ (0.7 × 3.2 × 2.9 × 2.5 m; Figure 2c) with an estimated weight of ~28.81 t. Approximately 50 and 23% of a-axis orientations of larger boulders (>5 t) were parallel (SW) and perpendicular (NE) to the reef edge, respectively. The orientations of the remaining 27% of larger boulders were SSE, including the three largest boulders observed (Figure 2).

Boulders were rarely seen at distances 0–30 m from the reef edge. Approximately 90% of boulders were distributed within the zone of 50–120 m, with a mode at about 82 m from the reef edge (Figures 2, 4). A wide range in boulder weight was also observed from <1 t to >28 t. However, the majority of weights were <2 t, with an estimated mean weight of 0.78 t (Table 3). Along the shore, the majority of boulders were distributed on the reef platform, but few large boulders (>1 m) were found in the back-reef moat. Larger boulders were usually found shorter distances to the reef edge (Figure 4). The boulder size distribution showed a Gaussian curve (Figure 3) as observed in previous studies (Deguara and Gauci, 2017). This indicates that boulders were most likely deposited by one type of wave (Figure 3; Barbano et al., 2010), because typhoon wave heights tend to exponentially decrease landward after breaking (Goto et al., 2009, 2010). The negative correlation between boulder size (or weight) and transport distance suggests that the morphology of the coast has great effects on patterns of boulder distribution. The roughness of the uneven reef platform caused by surface pits and detachment scarps inhibits landward transport (Naylor et al., 2016). However, a seaward fining trend (Figure 4) was also found at distances from 0–65 m, implying that backwash flow may play a key role in transporting boulders in this zone. Few large boulders at distances >160 m from the reef edge, indicating past high-energy waves that struck Hainan Island, lacked sufficient force to displace boulders (with long axes > 1 m) farther inland (Figures 2–4). In addition, the carbonate boulders show uneven distribution along with the reef platform and are mostly concentrated in the middle area. This should be mainly related to hydrodynamic field characteristics (e.g., storm energy) and morphology (e.g., reef margin slope) (Wang et al., 2021). Initial terrain modulates the characteristics of the hydrodynamic field and leads to the strengthening of hydrodynamic force in certain areas (Wang et al., 2021).

Boulder shape parameters varied significantly at different distances from the reef edge (Table 4), falling within various categories (Figure 5). As shown in Table 4, MPS decreases from 1.91 to 0.91 as the distance from the reef edge increases from 50–80 to 110–150 m. This indicates that the number of spherical boulders gradually diminishes landwards. Clear decreasing trends were similarly found in OPI (0.26 to −0.64) and DRI (0.87 to 0.64) at distances increasing from 50–80 to 110–150 m, respectively. This means that the numbers of prolate and disc-shaped boulders increase landwards. This is probably because disc-shaped boulders settle more slowly than other boulder shapes (Flemming, 1964), and are therefore more likely to be carried farther from the reef edge (Wang et al., 2004). Boulder shapes distributed at our XDH study site, therefore, showed a reasonably flat trend.

Figure 6 shows the lower limit of wave-induced current velocities during a typhoon event that can initiate large boulder transport for each transport mode. At the XDH study site, sliding or rolling is considered the most probable transport mode, since 97% of observed boulders were carried <0.5 m above sea level (Table 3). This is supported by the numerical results (Figure 6), indicating that minimum current velocities (MCVs) of 2.41–5.71, 2.32–5.50, and 4.86–14.74 m/s are sufficient to initiate boulder movement by rolling, sliding, and lifting, respectively. Although some large boulders may be transported by lifting, rolling is considered the most likely mode of transport for the majority of boulders because the elevation of CBDs is almost near the sea level and the inverted orientation of numerous boulders, which satisfies the condition proposed by Terry et al. (2016). Based on our assumption, we determined 2.41–5.71 m/s was the range required to transport these boulders (Figure 6; Table 5). In particular, boulders may be moved step by step in a series of many waves during typhoons, with estimated velocities representing the minimum current velocities experienced.

The nine ages of the 12 largest CBDs (Tables 2, 3) range from 1721 BC to AD 1756 (except modern age samples), with accumulated age uncertainties typically varying between 0.6 and 15% of the absolute age. Five samples were dated within the last 2,000 years, with four young boulders dated to the last 1,000 years. The earliest documented large wave events on Hainan Island can be traced back to AD 982. This allows us to compare historical and geological records over the past millennium. The age of XDH-1001 (AD 1756 ± 14) may correspond to large typhoon events in AD 1751/1753 (Table 2). The devastating typhoon of 1751 is recognized as causing one of the worst floods on the east coast of Hainan Island. Descriptions record that “Strong wind, powerful waves overflow onto the shore and Wanning city was flooded, numerous people died” (Chen, 1995). Similarly, the severe typhoon of 1753 was also extremely violent, producing strong waves that swept across Hainan Island, especially in the south (Chen, 1995). Descriptions record the following conditions: “Extraordinarily strong winds and rainstorm, typhoon flood onto the village, hundreds of great trees and crops blow down in Sanya city; more that 50% houses destroyed, 27 people were drowned by typhoon flood in Wanning city; buildings destroyed, trees and crops all broken, in Haikou city.” The age of XDH-933 (AD 1436 ± 90) may correspond to the catastrophic typhoon of AD 1524 (Table 2). The typhoon ravaged the southeast coast of Hainan Island and is considered the greatest calamity during the Ming dynasty (Chen, 1995). The description records that “In August, extraordinarily strong wind and rainstorm, typhoon floods over the land for a distance of more than 10 km, several ships washed up on land more than 1 km, more than 80% houses destroyed, countless livestock and people were drowned, never seen or heard of before by old people in Wanning city” (Chen, 1995). Given the errors on the U/Th dating results, the large typhoon in AD 1524 is consistent with our measured boulder ages. Moreover, some large boulders have been moved by modern typhoons according to the survey of local residents, such as the No. 648 boulder was deposited during landfalling typhoon Pamela over the Hainan Island in 1972 CE. This result was also consistent with the No. 648 boulder U/Th age.

Table 3 shows that the largest boulder detached from the reef edge had an estimated mass of 28.81 t. As specific signatures have not yet been developed to differentiate categorically between typhoon and tsunami boulders, it is important to adopt multiple factors (e.g., coastal geomorphology, boulders spatial distribution, and history of typhoon activities) to determine the reliability of typhoon reason (Goto et al., 2010; Cox et al., 2012; Nakamura et al., 2014; Lau et al., 2016). Previous investigations showed that maximum boulder weights displaced from the reef edge by typhoon flood inundation are usually ≤ 100 t (Goto et al., 2009), although more recent studies have identified larger transported clasts up to ~200 t (May et al., 2015). These numerous data indicate that typhoon flood inundation plays a significant role in producing boulder deposits in the region.

Morphometric investigations show that carbonate boulders with different sizes have been transported at varying distances from the reef edge (Figures 2a, 4). An inverse pattern between boulder size and transport distance has been identified at our XDH study site (Figure 4), a type of distribution often observed elsewhere (Goto et al., 2009, 2010; Deguara and Gauci, 2017) and associated with typhoon wave deposits (Goto et al., 2010). Goto et al. (2010) and Barbano et al. (2010) suggested that boulder distributions following multiple Gaussian curves, rather than falling into a single group, would be associated with different types of waves. Thus, the single distribution at our XDH site is evidence for one wave type. The flat trend of the boulder distribution (Table 4) might also be produced by a single wave type. Furthermore, the predominant SSW to WSW alignment of long axis orientations for larger boulders (>5 t) (Figure 2) is consistent with the strong impact of typhoon waves reaching the shoreline (State Oceanic Administration, 1999). In addition, the observed landward-fining trend of boulders sizes (Figures 3, 4) requires the action of multiple waves (Goto et al., 2009; Deguara and Gauci, 2017), indicating that coastlines of the study region have been exposed to frequent high-energy events. Taking all these morphometric characteristics into consideration, we reasonably think that typhoon events are responsible for the spatial distribution of boulders. This is also supported by the result of analyzing the theoretical distribution pattern of storm boulders on the XDH site (Wang et al., 2021). In addition, there are few historical records of tsunamis on the Hainan Island coastline compared with the northern South China Sea (Liu et al., 2007; Mak and Chan, 2007). In addition, the measured dates of large CBDs during the past 1,000 years are all consistent with historical catastrophic typhoons (Table 2; Chen, 1995; above section). Furthermore, the only observed tsunami event occurred on January 4, 1992, generated by an Mw3.7 earthquake offshore from southwestern Hainan. The maximum tsunami wave recorded in Yulin Bay was only 0.78 m high (Ye et al., 1994).

Numerical models of current velocity indicate that MCVs of 2.41–5.71 m/s are required to transport the observed boulders (Figure 6). Similarly, on Lanyu Island off southern Taiwan, Nakamura et al. (2014) calculated MCVs of 1.4–16.9 m/s for CBDs found at distances up to 120 m from the shore. In comparison, investigations of transported carbonate boulders on Ko Larn Island in southern Thailand by Terry et al. (2015) showed that boulder transport occurred at MCVs of 2.7–7.1 m/s, whereas CBDs carried up to 4.7 m above mean sea level on Ko Samui Island indicate past MCVs of 2.3–8.6 m/s (Terry et al., 2016). In the Philippines, Kennedy et al. (2017) revealed how the passage of Super Typhoon Haiyan in 2013 produced extreme high-energy waves that were able to displace a giant boulder of volume 83 m³ and mass 208 t, the transport of which required MCVs up to 8.9–9.6 m/s (May et al., 2015). Our MCV calculations from carbonate boulder deposits at XDH are therefore consistent with data from other typhoon-impacted coastlines across tropical Asia–Pacific region. Combining such evidence, it appears that the accumulation of carbonate boulders at our XDH site is the result of high-energy large typhoons.

Coral boulders provide key insight into the relationship between basin-wide large typhoon activity and paleoclimate change (Figure 7). The relationship between boulder ages and paleoclimate variability may be revealed by examining the records from different sites. In combination with the previous findings (Zhou et al., 2019b), we propose that there have been seven major periods of large typhoon periods (T1–T7) on the coast of Hainan Island (Figure 7), i.e., 1800–1500 BCE (T7), 1200–900 BCE (T6), 50–120 CE (T5), 550-800 CE (T4), 900-1,000 CE (T3), 1350–1900 CE (T2), and 1910–2000 CE (T1).

In the northwestern Pacific and South China Sea regions, modern observations and paleotempestological studies suggest that ENSO plays a significant role in typhoon longitudinal variations, with more movement routes directed toward Korea and Japan during El Niño times, but toward the northern South China Sea during La Niña times (Chan, 1985, 2005; Woodruff et al., 2009; Chen et al., 2012; Zhou et al., 2019a,b). However, little is known about the relationship between the major typhoon activities and climate dynamics at centennial-millennial scales. Comparisons between the patterns of the two routes suggest a general inverse correlation (Figures 7A,C,D). For example, the seven periods with intense, frequent typhoons in Hainan Island (T1–T7) match the periods of more La Niña-like times at centennial-millennial scales. In contrast, more active typhoons in Japan and Korea (except T3 and T7) occur mainly during El Niño times (Woodruff et al., 2009; Yang et al., 2017; Figures 7A,H). This confirms that the longitudinal variations of the intense typhoons are mainly regulated by ENSO patterns over the past 4,000 years (Zhou et al., 2019b).

However, two intense typhoon periods (T4 and T6) do not follow the pattern of ENSO, indicating that other climate drivers may be involved. Previous studies and modern observations both suggest that TCs variations in the WNP and SCS also greatly regulated by the thermal state of the western Pacific warm pool (WPWP) (Yue et al., 2019; Zhou et al., 2019a,b). Figures 7C,H show that almost all typhoon periods in the SCS match the warmer WPWP. This correlation between intense typhoon periods and warm sea surface temperature (SST) of WPWP phases has been observed by modern observations (Chen and Huang, 2008). During the warmer years of the WPWP, more TCs would be formed in the western part of WNP, facilitating TCs landfall in the SCS (Chen and Huang, 2008).

It should be noted that the increased typhoon activity occurred both in the northern SCS and eastern WNP (Korea, Japan, East China Sea) during the period of T3 and T7, suggesting that intense typhoon activity in the SCS and WNP should involve other climatic factors (Figures 7A,C). Observation records indicate that the active phase (northward) of Intertropical Convergence Zone (ITCZ) in the WNP could generate more TCs than that of inactive phases (southward) in the WNP (Ma and Chen, 2009; Cao et al., 2012). Therefore, it is hypothesized that more or less TCs and more northerly or southerly TCs genesis location likely leads to more typhoons move to both the northern and southern Northwestern Pacific during the strength ITCZ (weak ITCZ), respectively, as indicated by Chen et al. (2019). The multi-sites in the SCS and WNP over the past 1,500 years provide us the opportunity to test the role of ITCZ in governing latitudinal large typhoon activity in the centennial to millennial time-scale (Figure 7). Comparisons of the large typhoon proxy records generally support the hypothesis model between Northern WNP (the regions of Korea, Japan, and the East China Sea) and Southern WNP (Hainan Island, Taiwan Island, southern SCS, and Thailand) in typhoon activity over the past 1,500 years. Increased typhoon intensity observed in coastal of northern WNP (Korea; East China Sea; Figures 7A–F) and southern WNP during the T3 and T7 periods generally match the strength (northward) ITCZ. In contrast, the relative quiescence conditions in Northern WNP (Figures 7A,B; Korea; Japan; East China Sea) during the T2, T4, and T6 periods match the period of increased large typhoons at the southern WNP (Hainan Island; Taiwan Island; SCS and Thailand Gulf; Figures 7C–F), with the weak (southward) ITCZ phases. Historical records from the eastern coast of Asia also support an ITCZ-driven inverse typhoon activity pattern in WNP (Chen et al., 2019). Thus, similar to observations and historical studies (Cao et al., 2012; Chen et al., 2019), the CBD records from the XDH site support the pattern of more large typhoons move to Hainan Island during weak ITCZ phases. However, frequent large typhoons north move to the East China Sea during the weak ITCZ phase (T1), which should be associated with the thermal state of Kuroshio current. Because of the warmer SST of Kuroshio current (Cai et al., 2012) post 1900, it may have strengthened the intensity of TCs when TCs in the East China Sea pass over the Kuroshio (Yu et al., 2019).

Our findings from carbonate boulder deposits at Xiaodonghai provide evidence of the potential hazard exposure to extreme typhoon events in this area and surrounding low-lying coasts. Frequent incidences of people and infrastructure being inundated by large typhoon events are known on Hainan Island over the past 1,000 years (Chen, 1995). In southern Hainan, the rapid development of coastal tourism since the late 20th century has greatly promoted the aggregation of a large number of facilities in low-lying coastal areas and the extension of small fishing villages. Now the resident population (e.g., in the Sanya city) is about 15 times larger than that of 1750 and dramatically increased since the 1970s (Hainan Procincial Bureau of Statistics, 1992). With the development of Hainan Island, new infrastructure (e.g., harbors) and the number of residential and commercial buildings facing the threat of inundation is growing considerably. Instrumental records show that the ITCZ position has a gradually equatorward movement since 1949 CE (Ma and Chen, 2009), which means that more TCs would be a movement to the southern WNP, especially the northern SCS (Hainan Island). Moreover, the warm state of the WPWP will strengthen the SST of the WNP and SCS with future global warming, and so the north expanded warmer SST of WPWP is to guide more intense typhoons toward Hainan Island (Figures 7C,G; Yue et al., 2019; Zhou et al., 2019b). Furthermore, the risk of catastrophic typhoons will increase in the future according to the sea-level rise predictions of global climate models (IPCC, 2012). Our results indicate that the storm boulder records can be used to obtain information on typhoon intensity and frequency and to examine the projections of future large typhoons under changing climate.