A long-term view of tropical cyclone risk in Australia

Natural hazard risk is assessed by leveraging, among other things, the historical record. However, if the record is short then there is the danger that risk models are not capturing the true envelope of natural variability. In the case of tropical cyclones in Australia, the most reliable observational record spans less than 50 years. Here, we use a much longer (ca. 6000-year) chronology of intense paleo-cyclones and, for the first time, blend this information with a catastrophe loss model to reassess tropical cyclone wind risk in Northeast Australia. Results suggests that the past several decades have been abnormally quiescent compared to the long-term mean (albeit with significant temporal variability). Category 5 cyclones made landfall within a section of the northeast coast of Australia almost five times more frequently, on average, over the late Holocene period than at present. If the physical environment were to revert to the long-term mean state, our modelling suggests that under the present-day exposure setting, insured losses in the area would rise by over 200%. While there remain limitations in incorporating paleoclimate data into a present-day view of risk, the value of paleoclimate data lies in contextualizing the present-day risk environment, rather than complementing it, and supporting worst-case disaster planning.


Introduction
The present-day natural catastrophe risk environment is often necessarily framed in terms of the available historical record. Risk models can be augmented using statistical and machine-learning techniques to include consideration of 'grey swan' weather events (Lin and Emanuel 2016) that may be unprecedented in terms of historical experience but foreseeable given our ability to extrapolate beyond the limits of our observations. Extrapolated or resampled extremes are, however, still reflective of the statistical properties of the observations on which they are based. Climate change presents a new challenge as it suggests that historical experience alone may no longer be a good indicator of future risk (Collins et al. 2012;Jewson et al. 2021). At a time when there is much effort focused on attributing the impacts of the recent warming trends on extremes (Stott et al 2016;Bellprat et al. 2019;Perkins-Kirkpatrick et al. 2022), it is opportune to consider the degree to which the climate system is inherently non-stationary (Obeysekera and Salas 2016;Slater et al 2021) and the risk envelope of longer-term natural climate variability may be larger than our recent experience suggests.
A powerful tool to help understand the dollar-loss impact of weather extremes is the natural catastrophe loss (CAT) model. CAT models are decision support tools widely used within the (re)insurance industry over the past several decades to help price natural hazard risk and for aggregate risk management. They are complex probabilistic models that frame risk in terms of the hazard, exposure and vulnerability. In this framework, risk is measured in terms of the fractional loss incurred on a property relative to its insured value due to damage from natural hazards. An overview of many of the approaches used in CAT modelling can be found in Grossi and Kunreuther (2005), Mitchell-Wallace et al. (2017) and Michel (2018). The hazard event set of a cyclone CAT model may comprise hundreds of thousands of synthetic tropical cyclone events with the aim of augmenting the sparse observational record using stochastic modelling techniques.
In the case of tropical cyclones in Australia, the observational record extends back to 1909, although due to varying rates and methods of data assimilation, is considered homogenous only since 1981 (BoM 2018; Courtney et al. 2021). This 40-year period samples predominantly from the positive (El Niño like) phase of the Interdecadal Pacific Oscillation (IPO). The IPO is a multi-decadal signal of climate variability in the Pacific and has been shown to exert an important influence on tropical cyclone frequency in Australia, by modulating the intensity of shorter-term regional modes of variability such as the Indian Ocean Dipole (IOD) and El Niño-Southern Oscillation (ENSO) (Power et al 1999;Grant and Walsh 2001;Magee et al. 2017). Positive (negative) IPO, El Niño (La Niña) and positive (negative) IOD are all associated with a reduction (increase) in tropical cyclone frequency along the east coast of Australia. In addition, external influences such as volcanism (Pausata and Camargo 2019;Altman et al. 2021) and solar activity (Elsner and Jagger 2008;Haig and Nott 2016) have also been suggested to drive periods of lower or higher cyclone frequency on long timescales, although their effect is less certain than shorter-term drivers (Camargo and Polvani 2019). Given the importance of low-frequency climate variability on tropical cyclone activity via the modulation of shorter-term modes, the observational record is almost certainly too short to capture the full envelope of natural variability in tropical cyclone risk in Australia.
Paleoclimate data comprise proxy observations of past extreme weather events that have left an imprint in the physical environment in some form. These imprint records can be reconstructed to infer timeseries of extreme weather events that far exceed the direct observational record in length. However, these records often only capture high magnitude 'tail' events which makes it difficult to reconcile with higher resolution observational data. Paleotempestology is a sub-field of paleoclimatology which uses geohistorical proxies to infer characteristics of past storm activity. Examples of paleotempestology proxies include overwash deposits preserved in coastal sediments, microfossils such as pollen or diatoms, or wave-and flood-generated sedimentary deposits in marine sediments, coral shingle or shore-parallel beach ridges (Bradley 2014).
The 'beach ridge method' is a paleotempestology technique used to derive multi-millennial proxy chronologies of intense paleo-cyclones and has been used along the tropical coastal areas of Australia (e.g. Nott and Hayne 2001;Nott 2003;Nott and Jagger 2013). At many coastal locations in the Australian tropics, strand plains comprising series of shoreparallel beach ridges have developed over the late Holocene period as a result of sea-level fall of approximately 1-1.5 m over this period (Lewis et al. 2012) and consequential progradation of the coastline. This means the oldest ridges are furthest landward from the present-day coastline (Fig. 1).
Sea-level rise since the end of the last glaciation (circa 20 ka BP [1000 years before present]) had overrun earlier storm deposits, meaning beach ridge paleo-data rarely extend back further than 6-7 ka BP when global sea level stabilized and locally, began to fall (Lewis et al. 2012).
Over this 6-7000-year period, it is assumed that each prominent beach ridge-over a certain elevation threshold and comprising marine-origin sediment-was produced by a paleo-cyclone. Aeolian-capped (wind-blown) ridges and ridges under 2 m AHD (Australian Height Datum) are disregarded as potentially having been formed by the prevailing trade wind-wave climate.
The elevation of each of the prominent, marine-origin beach ridges represents the minimum height at which the paleo-cyclone storm tide (i.e. astronomic tide plus meteorological storm surge component) must have reached to have produced the ridge. By dating each relevant beach ridge using Optically Stimulated Luminescence (OSL, see Bateman 2019 for a review), a ca. 6000-year chronology of paleo-cyclone storm tide heights are derived.
There has been some discussion in the literature on the validity of the beach ridge method as an indicator of paleo-cyclones, centring on the role of aeolian transport and the possibility of tsunami processes also leading to ridge development (Tamura 2012(Tamura , 2014Nott 2013). Tamura (2014) acknowledged that the ridges were likely largely deposited by wave/storm tide action. Nott (2013) further highlighted that these ridges can be initiated by Fig. 1 Beach ridge plain at Cowley Beach, south of Cairns. Darker colours represent swales between the lighter colour ridge crests; black line represents shore-parallel transect taken through ridges to derive paleocyclone chronology. Figure modified from Nott et al. (2009) aeolian activity but once they reach a height of approximately 1.5 m (AHD) aeolian processes no longer play a role.
It is also important to note that there may be multiple cyclones contained in a single ridge, and thus cyclones of lesser magnitude are not counted (Nott and Jagger 2013). We focus on the last (highest in the stratigraphy) unit to be deposited onto each ridge within a ridge plain. Hence, our records are censored, meaning we are dating only the most extreme events. Extensive observations indicate that the most extreme events are overwhelmingly depositional, and not erosional, as this is the net response of the cross-shore profile when the level of inundation is higher than the existing first ridge (Dolan and Godfrey 1973;Sallenger 2000;Nott and Hubbert 2005;Fritz et al. 2007;Nott 2011). We deliberately miss most events that occurred to develop a ridge plain. As a result, this study only focusses on the tail events. The way in which we deal with the rest of the frequency-magnitude distribution within a contemporary loss modelling framework is described in Sect. 3.1. Here, we use beach ridge chronologies compiled by Nott and Hayne (2001) and Nott (2003) at five spatially discrete coastal sites in North Queensland, Australia, to reassess tropical cyclone wind risk in this area. The physical parameters of paleo-cyclones have been inverse-modelled prior to this study from storm surge heights derived by Nott and Hayne (2001) and Nott (2003). In this previous work, hydrodynamic modelling was undertaken to determine the conditions necessary to generate a marine inundation equal to or greater than the height of the ridges. We use the model results of Nott and Hayne (2001), Nott (2003) to calibrate an Australian tropical cyclone loss model-used by the (re)insurance industry-to reproduce the derived frequency/magnitude distribution of tropical cyclones along this section of coast according to paleo-proxies and calculate property losses accordingly. We then compare these paleo-derived loss estimates to those using the contemporary observational record in the tropical cyclone model to examine present-day cyclone risk within the context of longer-term natural climate variability.

Study area and datasets
The available paleo-data cover a 400-km length of coastline in North Queensland, Australia, extending from approximately 15 km north of Port Douglas to 30 km south of Townsville (Fig. 2).
This area covers the three main population centres of Cairns (150,000), Townsville (180,000) and Mackay (130,000). The relatively high frequency of cyclone events, lowlying coastal topography and moderate exposure (with reference to the built environment) means this region makes a significant contribution to the national tropical cyclone risk. Within this area, five spatially discrete field sites were used. These were (north to south): Supporting Information further details the beach ridge record at each site and exclusions of particular ridges in the series.

Sea-level fall
Since approximately 5 ka BP, paleoclimate indicators suggest that sea level has fallen by approximately 1.5 m in Queensland relative to the present-day mean sea level (Lewis et al. 2012). At the end of the Last Glacial Maximum, approximately 20 ka BP, sea levels began to rise rapidly as temperatures warmed and ice melted-a period known as the Holocene Marine Transgression. A sea-level high stand was reached between 8 and 5 ka BP, after which hydro-isostatic adjustments (elastic response of the continental crust to varying water loads) led to a gradual sea-level fall of approximately 1.5 m to present day. A falling sea level may mean that older beach ridges were not produced by a storm surge relative to the present-day mean sea level; rather, the surge may have been relatively smaller (and the tropical cyclone relatively weaker) and the mean sea level higher. Indeed, many (but not all) of the shore-normal dune transects exhibit a gradual decrease in beach ridge height in a seaward direction ( Figure S1), suggestive of the influence of gradual sealevel fall on ridge elevations.
The beach ridge elevations at each site were adjusted to 'correct' for this effect by assuming a linear sea-level fall of 1.5 m over 5000 years. If the chronology extended beyond 5000 years, the linear trend was extended back to the oldest date in the timeseries.

Wave exposure
The results of Nott and Hayne (2001) and Nott (2003) were used in this study to build linear relationships between TC central pressure (CP) to beach ridge elevation for a range of CPs for each site. At most sites modelled in the studies referred to in Table 1, the beach ridge elevation was assumed to be equivalent to the flow depth, i.e. the astronomic tide plus storm surge. However, at exposed coastal locations, the height of storm tide inundation during a TC event is equivalent to the astronomic tide plus storm surge (inverse barometric effect and wind setup) plus wave effects.
At wave-exposed, open coast locations, wave effects can add a considerable amount on top of a storm surge. For example, Mortlock et al. (2018) found that on average waves added 16% to water levels along the Mackay coast during TC Debbie in 2017, but this varied considerably between sheltered (2%) and exposed (> 50%) locations.
At each site, the exposure to waves was assessed based on the beach ridge morphology and coastal planform orientation relative to the TC angle of approach that produced the largest surge. We make the assumption that the coastal planform orientation at the time of cyclone landfall is similar to that of today, with respect to the line drawn perpendicular from one end of the embayment to the other. This is because, despite undoubted wave climate variability during this period of time, the directional modal wave climate has a limited impact on shoreline planform at these sites because of the protection afforded by the Great Barrier Reef (GBR) matrix. Details of the wave exposure analysis are provided in Supporting Information.
At sites with little to no wave exposure, the elevation of the beach ridge was assumed to be equal to the surge flow depth (i.e. with no additional wave component on top). At sites with high wave exposure, the beach ridge height was assumed to represent the combined height of the storm surge with additional wave effects. This means that for the waveexposed sites, the surge component and TC intensity needed to be reduced.

3
The reduction factor for wave-exposed sites was estimated using numerical wave modelling undertaken by Nott et al. (2009) at Cowley Beach. Further information on the wave reduction correction is provided in Supporting Information. The tropical cyclone central pressure was then calculated based on the wave-reduced ridge heights.

Tidal uncertainty
The estimated CPs have uncertainty margins which are associated with not knowing the height of the tide during the paleo-TC event. For example, if the TC made landfall at high tide the surge component (and TC intensity) necessary to equal the beach ridge height would be smaller than if the TC had crossed the coast at low tide.
Tides follow an 18.6-year lunar tidal cycle, which can be predicted at any location from astronomic constituents. The 18.6-year tidal constituents were derived for each survey site from which the 2σ probability tidal range was defined. This tidal range represents 95% of all possible tidal states at each site. It is larger than the mean spring tidal range but does not include the highest and lowest astronomical tides (HAT and LAT, respectively). The 2σ tidal range and associated CP values for each site are given in Table 2.
The influence of tidal state on paleo-cyclone CP estimates was represented as a cosine function with an amplitude defined by the 2σ tide range and the associated CP values given in Table 2 for each site.
For each paleo-cyclone event, the cosine was randomly sampled 10,000 times at equally-spaced intervals between two consecutive high tides. The median value from these 10,000 samples was then used as our estimate of paleo-cyclone intensity. The 5th and 95th percentiles, estimated from Monte Carlo sampling, were used to represent tidal uncertainty margins around each CP estimate. This approach accounts for the random coincidence of TC landfall and tidal state, rather than assuming the tidal state is the same for all events. Other sources of uncertainty in this study are discussed in Sect. 3.3.

Average recurrence intervals
A generalized extreme value (GEV) distribution was fitted to the paleo-CP estimates at each site to estimate average recurrence intervals (ARIs) to 50,000 years (the length of the stochastic TC event set in the probabilistic loss model, CyclAUS, discussed in Sect. 2.6). The 5th and 95th percentiles obtained during the randomization of the tide represent the uncertainty of a given event. The GEV fits and ARI estimates with 95% confidence bounds for each site are given in Supporting Information. In each case, the GEV distribution was selected based on a two-sample Kolmogorov-Smirnov (K-S) test. The two-sample K-S test is a nonparametric test of the equality between an empirical and fitted cumulative distribution function (cdf) of the data. The null hypothesis is that the samples are drawn from the same distribution. The two-sample K-S test can also be modified to serve as a goodness-of-fit test between fitted distributions. Here, the maximum absolute difference between the cdfs of the empirical and fitted distributions is used as the test statistic (Massey 1951; Marsaglia et al. 2003). The fitted distribution that returns the lowest maximum absolute difference when compared with the empirical cdf, is selected as the 'best-fitting' distribution. For extrapolating extreme values, an empirical cdf cannot be used as it is finite. Extreme value distributions, by comparison, are continuous and allow for extrapolation. Weibull, Gumbel, Fréchet, Exponential and Pareto distributions were tested for (see Kotz and Nadarajah 2000, for an overview on these extreme value distributions). The ARIs at each site represent the expected value of time in years between single tropical cyclones exceeding a certain intensity at that site, or within a specific radius of the site. The radius over which these values are relevant can be estimated from the local geomorphology, based on the alongshore extension of the beach strandplain. The distance north and south of the east-facing sites-where surveying the same transect would return broadly the same ridge sequence-is indicative of the spatial extent of the homogenous paleo-TC chronology.
The site geomorphologies extend between 3 and 50 km with a combined shoreline length of approximately 100 km. The five sites cover a total shoreline length of 600 km. They are all spatially discrete, i.e. none overlap. Assuming the TC risk environment is broadly consistent in this area, we can assume that the paleo-data represent a sixth (100/600 km) of the number of events that would be expected to make landfall over the total (600 km) shoreline over the time period sampled (~ 6000 years).
This assumes that each site contains a unique set of paleo-cyclone storm tide markers in the sedimentary record, i.e. that no one single event caused a paleo-ridge at multiple adjacent sites. While there is no way of knowing this with certainty, contemporary cyclone observations suggest this would be unlikely. The distance between adjacent sites range from between 50 km (Rockingham Bay to Cowley Beach) to 170 km (Cowley Beach to Wonga Beach). Observations from TC Debbie, which crossed the North Queensland coast as a Category 3 cyclone in 2017, indicate the maximum storm surge occurred approximately 50 km south of landfall, with sites as little as 20 km adjacent to the location of surge maximum experiencing over 50% less surge (Mortlock et al. 2018). Given our neighbouring sites are at least 50 km apart, the chance of a single event leaving an extreme paleo-storm tide ridge at multiple sites is considered low.
To generate ARI estimates relevant for the whole coast between Cairns and Townsville, we pooled the five site CP estimates together (114 events), fitted a single GEV curve through all of them, then rescaled to a time interval six times longer. A new GEV curve was then fitted through the resampled data. This increased the probability of landfall occurrence from a site-by-site basis by an approximate factor of six, consistent with the total shoreline length of analysis.

Tropical cyclone loss model
CyclAUS is a stochastic tropical cyclone wind loss model used widely within the reinsurance market in Australia (e.g. Haig et al. 2014). Three modules are incorporated to simulate the TC hazard and calculate resulting losses. First, a TC generation module, which contains historical information and TC data to generate TC events that reflect the varying probable magnitude and frequency for events within the study area. Second, an exposure module which contains data on property location, age and concentration. Third, a vulnerability module, which contains functions that relate wind speed to property damage level. Vulnerability functions vary based on building design and age.
The TC generation module is based on the characteristics of approximately 300 historical TC events since the 1969/70 austral season. Monte Carlo sampling of log-normal probability distributions is used to initialize synthetic TC events at 'gates' and in 'cells' around the Australian coast, after which autoregressive modelling is used to generate a 50,000-year event set of approximately 370,000 synthetic TC tracks and parameters. The maximum 3-s gust wind speeds of each event are mapped onto a variable resolution grid with the highest resolution of 90 m reserved for the most populous areas of the coast.
In this analysis, a market portfolio representing the total insured value across the four major population centres in the study area (Mackay, Bowen, Cairns and Townsville) was used. In our use of CyclAUS here, we only modify the hazard component of the model to reflect the longer-term cyclone activity regime based on the paleo-cyclone ARI estimates. We make no modifications to the exposure or vulnerability components, to represent the present-day 'as-is' property distribution and characteristics.

Model validation
To validate the modelled representation of recent cyclone activity, the modelled frequency/ magnitude estimates from CyclAUS were compared to the recent observational record for the Cairns to Townsville area. The landfalling CP values of each TC event were extracted from the Bureau of Meteorology (BoM) historical catalogue and a GEV distribution was fitted to the tail of the observations (Fig. 3a). Fig. 3 A Comparison of landfalling TC ARI estimates based on the historical event catalogue (green circles), and extrapolation to 50,000 years using a GEV fit (black line), and ARI estimates from the CyclAUS 50,000-year stochastic event set (grey circles). B Comparison of landfalling TC ARI estimates using the CyclAUS event set (grey circles), the ~ 6000-year resampled paleo-data (red circles) and the adjusted CyclAUS event set (blue circles). Also shown are the 95% confidence intervals for the paleo-GEV curve fit at mean tide (dotted black lines) and the paleo-GEV curves for high and low tide scenarios (dotted red lines). The paleo-data begins at an ARI of 50 years. In both (A) and (B), solid black horizontal lines denote the central pressure of Category 2 and Category 5 tropical cyclones based on the Australian Cyclone Severity Scale

3
The comparison of the BoM curve to the paleo-curve was first limited to the study area (~ 600 km shoreline length, approx. Bowen to Cooktown, see Fig. 2) but it was found that a lack of historical data in this area hindered a robust extreme value analysis. The comparison was expanded to include ~ 1000 km shoreline length (approx. Mackay to Cape Melville National Park), which increased the sample size from 26 to 43 landfalling events. Figure 3a shows that there is a good comparison between the fit of the historical observations (black line) and the CyclAUS event set (grey circles) from the 5 to 50,000-year ARIs. The annual frequency of landfalling events also compares well (0.89, BoM record, and 0.88, CyclAUS). These results suggest that the model is replicating the present-day risk environment well.

The average annual loss and occurrence exceedance probabilities
To quantify the change in risk between the present-day risk environment (reflective of the past ~ 50 years of cyclone observations) and the longer-term risk environment (reflective of the past ~ 6000 years from paleo-proxies), we use the change in the average annual loss (AAL) metric, expressed in Australian Dollars (AUD). The AAL is the average annual expected insured loss due to property damage resulting from tropical cyclone winds and associated water ingress. Here, we use the aggregate AAL for the study area-that is, the AAL expected across all of the modelled addresses. The AAL within the model is calculated as the sum of all cyclone events that produce a loss, divided by the number of synthetic years in the simulation (in this case, 50,000). By adjusting the frequency of cyclone events within the model based on the paleoclimate analysis, the number of loss-producing events changes between the present-day/longer-term model runs, resulting in a differential in the AAL. This differential is the primary metric we use for expressing the change in risk on the present-day exposure environment between our short and long-term view of the tropical cyclone hazard.
In addition to the AAL, we also compare losses at various occurrence exceedance probabilities (OEPs) out to an ARI of 1000 years. An OEP is the probability that the largest loss each year exceeds a certain amount of loss and the relationship to ARI is given by OEP = 1/ ARI. For example, an OEP of $5 million at an ARI of 100 years means that there is a 1% chance each year of a single cyclone event causing a property damage loss greater than $5 million.
The AAL is a useful metric for assessing the technical insurance premium, in other words, the cost (and affordability) of insurance to the insured. OEPs can be used to understand disaster management requirements or critical infrastructure planning for government, or capital retention and risk transfer requirements for (re)insurance. Simply put, change in ("delta") AAL is a metric for mean shift changes, whereas delta OEP at high ARIs is a metric for changes in tail risk.

Comparison of present-day and long-term cyclone probabilities
Recurrence estimates from the present-day CyclAUS event set were compared to those derived from the paleo-data (Fig. 3b). After sensitivity testing, an adjusted event set was generated that better matched the paleo-data for this area (blue circles, Fig. 3b).

3
The event set was adjusted by varying intensity modelling within CyclAUS to produce a curve with frequency/intensity characteristics resembling the paleo-record. Intensity modelling refers to the rate at which TCs intensify in the model domain (i.e. the drop in central pressure per time step) after they are initiated.
It was found that increasing the rate of TC intensification per timestep by up to 100% (0% increase at ≤ − 20 hPa changes up to a 100% increase at 0 hPa change) produced the best fit to the paleo-data. At very long ARIs (1000-50,000 years), the modified CyclAUS event set falls below the paleo-curve but is still within the bounds of tidal uncertainty (Fig. 3b). At these very low probability ranges, differences in the tail of the curve are unlikely to significantly impact AALs because there are very few of them in the event set.
In Fig. 3b, the paleo-data begin at an ARI of 50 years. This is because the ARI surveyed at each study site was around 300 years (Table 1), which equates to an ARI of around 50 years for the whole study area due to the multiple sites. Recurrence estimates < 50 years, therefore, are not based on paleo-evidence. The impact of the paleo-adjusted event set on expected property damage losses is discussed in Sect. 3.2.
We also compared the TC landfall probabilities using the observational (~ 50-year) record, to those using the paleo-data (~ 6000-year) (Fig. 4). The lack of data in the observational record (26 events in total, 8 which follow an extreme value distribution) means there is very low confidence in fitting and extrapolating to high ARIs. In addition, extrapolation of extreme events is generally considered reliable for periods up to three times the data period (Harper 1996) which restricts the use of the observational record to the 150-year ARI level. By comparison, the larger amount of available paleo-data provides a better confidence in the extrapolation, although the uncertainty associated with not knowing the state of tide at the time of the paleo-cyclone landfall remains a significant limitation.
In Fig. 4, the data have been extrapolated forward (observational) and back (paleo) to the TC Cat. 5 threshold to compare estimates at this level. The observational record suggests that Cat. 5 TCs can be expected to make landfall somewhere within the study area, on average, once every ~ 90 years (based on the past ~ 50 years of data). The paleo-record indicates that, over the late Holocene period (past ~ 6000 years), Cat. 5 TCs have made landfall in this area with much greater frequency-on average, once every ~ 20 years. It remains unknown whether this is due to an increase in TC intensity or TC frequency. For Fig. 4 Comparison of landfalling TC occurrence estimates in North Queensland using the observational ~ 50-year record and the ~ 6000-year resampled paleo-record. The GEV fit and 95% confidence intervals (CI) are shown for both datasets. The Cat. 5 TC threshold on the Australian Cyclone Severity Scale (929 hPa) is shown to compare ARIs at this level 1 3 instance, a greater number of intense tropical cyclones making landfall may be due to an increase in the intensity of all tropical cyclones, an increase in the total number of tropical cyclones forming, or both.
In summary, this work supports the assertion of Nott (2003) and Nott et al. (2007) that the short (~ 50-year) observational record under-samples the long-term landfall probability of high intensity TCs in North Queensland. This is consistent with other studies that suggest that the past 50 years have been unusually quiescent for TC activity in the context of the past few centuries to millennium (Callaghan and Power 2011;Haig et al. 2014).

Changes in expected property damage losses
The paleo-adjusted event set was run through CyclAUS to compare estimates of insured property damage losses within the study area with the present-day estimate. The model was run on a market portfolio, which is representative of the total present-day insured exposure in the area and does not represent the results of any single insurer. Since we are not concerned with absolute values here, and to further preserve anonymity, we only report the estimated relative change in AAL and OEP, when using the paleo-adjusted view of risk compared to the present-day view of risk (Table 3).
Results show that there is a significant increase in the level of modelled insured losses for the study area when incorporating the ~ 6000-year paleo-cyclone data, compared to the present-day estimate based on the recent observational record. Table 3 shows there is a ~ 250% increase in the AAL estimate for this area. The difference in losses is most pronounced at shorter ARIs, from around + 280% at the 10-year ARI to around + 200% at the 1000-year ARI.
These large differences reflect the substantial downward shift in the TC hazard recurrence curve when fitted to the paleo-data (Fig. 3b), difference between blue and grey circles). For frequently occurring natural hazards, any change to the AAL is usually driven in the most part by changes at lower ARIs (i.e. less intense, but more frequent events). This is because the AAL is the integral of all losses in the event set, and losses typically accumulate faster at lower ARIs because there are more of them in the event set, than at the tail, where a few extreme events do not add a significant amount to the AAL. Since these smaller-scale events are not preserved in the paleo-record, the curve was extrapolated back from evidence of the most intense paleo-cyclones that left a lasting marker on the coastal landscape.
This approach assumes that changes in the frequency of high magnitude landfalling TCs are synonymous with changes in lower magnitude events. However, both the TC frequency and magnitude distributions can change simultaneously with climate forcing. For example, projections of near-future changes in TC activity are for an increase in intensity but decrease in frequency of cyclogenesis in the Southwest Pacific region (e.g. Knutson et al. 2020).
While there is a significant level of uncertainty associated with this comparison, it is consistent with the existing body of literature (e.g. Callaghan and Power 2011;Haig et al. 2014) that suggests that we are currently experiencing an abnormally quiescent period of TC activity along the Australian east coast compared to the longer-term average. This suggests that present-day loss estimates may be towards the lower end of what could be reasonably expected over the longer term within the bounds of natural climate variability experienced over the past 6000 years.

The value of paleoclimate data for tropical cyclone risk estimation
Paleo-records can be useful for providing a better estimate of the 'worst-case scenario' than observational data because very long records are more likely to sample very rare, catastrophic events with long recurrence intervals. Paleo-records can also be useful for providing analogues of past 'hyperactive' and 'inactive' periods of TC activity, as this study has shown.
Here, the beach sedimentary record was not of a sufficient temporal resolution to subsample periods of TC hyper-and inactivity. Instead, we compared TC recurrence estimates using the full ~ 6000-year paleo-record to the past ~ 50 years of observations. While this highlighted, importantly, that the short historical record under-samples the landfall probability of high magnitude TCs compared to the long-term 'average', a direct comparison such as this assumes stationarity in the climate system and an invariant probability distribution over this ~ 6000-year period. However, because TCs are a manifestation of constantly shifting climatic conditions, the spectrum of recurrence interval risk is also nonstationary (Frappier et al. 2007). This suggests that if the same analysis was performed for the same area over a different time period, results would be different.
For paleo-data to be useful for contemporary risk modelling, they should capture information on both the time-varying intensity and frequency of tropical cyclones and ideally, be of sufficient temporal resolution to resolve the full magnitude spectrum of events. In practice, it is rare to obtain a record that is both of high temporal resolution and that provides magnitude/frequency information. While the beach ridge method provides magnitude/frequency information, it does so at a low temporal resolution (~ 300-year event intervals). This means it only provides direct information on the tail of the distribution (the high magnitude, infrequent events). Figure 4 shows, for example, that the paleo-data only record the long-term occurrence of Cat 5 events landfalling in this area, and no information on Cat 4 s and below.
However, if we assume that the same physical relationships observed today are relevant for the late Holocene period in general, we can use the Cat 5 paleo-cyclone data to infer the rest of the distribution. As shown in Fig. 3b, we can fit the remainder of the distribution by generating a long-term stochastic TC event set using a tropical cyclone loss model that best matches the paleo-data available. This approach potentially negates the need for explicit paleo-proxy information on low-magnitude, more frequent TC events and makes coarse resolution paleodata more readily applicable for modelling the entire risk spectrum.
Lastly, a limitation of using the beach ridge method for risk modelling is the multiple assumptions-and high uncertainty-associated with deriving CP estimates from the ridge elevations. Similar sources and magnitudes of uncertainty also exist in other paleoclimate analyses. We can summarize the sources of uncertainty in this study as follows: 1. The possibility that paleo-cyclone deposits have not been fully isolated from aeolian or paleo-tsunami processes in the beach ridge timeseries (see Tamura 2012Tamura , 2014Nott 2013). 2. The margin of error of the Optically Stimulated Luminescence (OSL) beach ridge dating method (which ranges from ± 10 to ± 400 years). 3. The homogeneity of sea-level fall over the past ~ 6000 years in the North Queensland area, and the contribution of sea-level variations on beach ridge heights. 4. The contribution of wave effects to total water levels at each site and during each paleo-TC. 5. The tidal state at the time of paleo-TC landfall. 6. The paleo-TC parameters (radius of maximum winds, forward speed, track angle relative to coast and distance of landfall from site).
As can be seen in Fig. 3b, tidal uncertainty has a large effect on the central pressure estimation (± 23 hPa, see Table 2). The underlying assumption of this method-that the beach ridge height is equivalent to the paleo-cyclone water level-can be regarded as a conservative estimate of the storm intensity because observations (Sallenger et al. 2006;Nott et al. 2009;Mortlock et al. 2018) have shown that TC water levels often exceed the coastal ridge that is produced by deposition of sediments after inundation. The selection of TC parameters in the hydrodynamic modelling by Nott et al. (various, see Table 1) also provides the most conservative estimate of the storm intensity required to generate the indicative level of inundation.
Even after all the limitations of paleoclimate proxies are considered, the question still remains as to whether we should be accounting for the longer-term variability if it does not reflect the present-day and near-future risk environment. For example, if periods of enhanced TC activity during the course of the late Holocene were driven by external forcing such as low solar activity (Elsner and Jagger 2008;Haig and Nott 2016) and enhanced volcanism (Pausata and Camargo 2019;Altman et al. 2021), this may not be a relevant analogue for the present. It is important to acknowledge these and other factors given the magnitude of the changes in loss and the implications of this for the insurance industry and beyond.
Instead, the real value of paleoclimate information for natural hazard risk assessment may lie in 'worst-case' scenario planning, and-perhaps more academically-contextualizing the current risk environment.

Conclusions
Results here and those of others (Callaghan and Power 2011;Haig et al. 2014;Goodwin et al. 2015) indicate that the observational period is unusually quiescent for tropical cyclones in Australia in the context of the past several millennia. As a result, the recurrence probabilities of cyclones differ considerably. For example, the observational record suggests that a Cat 5 event has a ~ 90-year recurrence interval on the North Queensland coast, whereas the paleo-record indicates that it has a ~ 20-year average period of recurrence.
The lowest recorded central pressure in Australia is 900 hPa, attained by TC Gwenda in 1999 and TC Inigo in 2003 (BoM 2021). However, recent work suggests that TC Mahinawhich made landfall on the Cape York Peninsula in Far North Queensland in 1899 (about 500 km north of our study sites)-may have reached a central pressure as low as 880 hPa . During this event, a 13 m storm tide was observed, and 300 lives were lost. This makes it the largest recorded storm surge in the world, and potentially the most intense tropical cyclone on record in the Southern Hemisphere.
It is no coincidence that only a very small portion of disaster spending (3% in Australia, 4% in the US) goes towards mitigation (Coppel and Chester 2014;Cigler 2017). The rest (i.e. > 95%) is spent on disaster recovery. The cycle of disaster, reaction and complacency has been well studied by economists and psychologists alike (Meyer and Kunreuther 2017). Protective actions, whether by individuals or governments, are usually designed to be adequate to the worst disasters experienced because, simply put, images of an even worse disaster do not come easily to mind (Kahneman 2011).
Paleoclimate data can help risk managers re-imagine the worst case and mitigate with this upper limit in mind. Our study suggests that a plausible worst-case scenario for North Queensland would be a landfalling TC in the Cairns area with a coastal-crossing central pressure of 870 hPa.
While the paleo-cyclone record may be useful in painting grew swans a lighter shade of grey, there is still some way to go before they can be incorporated into contemporary risk modelling. The first limitation is the cascading uncertainties in the process of inverse modelling cyclone parameters from paleo-proxy observations. The second limitation is the inability of paleo-cyclone data to capture a chronology of both small, frequent events and the large, infrequent events. For example, the beach ridge data used here had an average spacing of 300 years at each site. To address this issue, we adjusted the intensity modelling in the cyclone loss model such that the probability distribution matched the paleo-data at the tail, therefore obtaining inferred information on the 'small, frequent' portion of the curve. However, this is less than ideal because we assume that the shift in the probability distribution is due to a change in intensity, and not frequency. For instance, a greater number of intense tropical cyclones making landfall may be due to an increase in the intensity of all tropical cyclones, an increase in the total number of tropical cyclones forming, or both.
In summary, this study represents the first published account of blending paleoclimate data with contemporary tropical cyclone risk modelling. We find that while there are some critical limitations of incorporating paleo-data into a present-day view of risk, the value of paleoclimate data lies in contextualizing the present-day risk environment, rather than complementing it, and supporting worst-case disaster planning.