Guillaume Dodet
email: guillaume.dodet@univ-brest.fr

Bruno Castelle

Gerd Masselink

Tim Scott

Mark Davidson

France Floc'h

Derek Jackson

Serge Suanez

Beach recovery from extreme storm activity during the 2013/14 winter along the Atlantic coast of Europe

The storm sequence of the 2013/14 winter left many beaches along the Atlantic coast of Europe in their most eroded state for decades. Understanding how beaches recover from such extreme events is essential for coastal managers, especially in light of potential regional increases in storminess due to climate change. Here we analyze a unique dataset of decadal beach morphological changes along the west coast of Europe to investigate the post-2013/14-winter recovery. We show that the recovery signatureis site-specific and multi-annual, with one studied beach fully recovered after two years, and the others only partially recovered after four years. During the recovery phase, winter waves primarily control the timescales of beach recovery, as energetic winter conditions stall the recovery process while moderate winter conditions accelerate it. This inter-annual variability is well correlated with climate indices. On exposed beaches, an equilibrium model showed significant skill in reproducing the post-storm recovery and thus can be used to investigate the recovery process in more details.

Introduction

Sand and gravel beaches may undergo dramatic erosion and recession during sequences of extreme storm wave events [START_REF] Ferreira | The role of storm groups in the erosion of sandy coasts[END_REF], leaving them in a state of morphological dis-equilibrium. A phase of 'recovery' towards pre-storm sediment volume is then a natural morphodynamic response to this depleted state [START_REF] Brenner | Characterizing storm response and recovery using the beach change envelope: Fire Island, New York[END_REF]. Because the rates of recovery depend on the magnitude of the storm-induced changes, the subsequent hydrodynamic conditions, the sediment availability, and the geological setting, predicting the time until full recovery is achieved (if ever) is challenging. The current predictions of climate change indicate an acceleration of sea level rise [START_REF] Cazenave | Sea level rise and its coastal impacts[END_REF]) and a poleward shift of midlatitudes storms [START_REF] Tamarin | The poleward shift of storm tracks under global warming: A Lagrangian perspective[END_REF], which will likely increase extreme water levels [START_REF] Vousdoukas | Global probabilistic projections of extreme sea levels show intensification of coastal flood hazard[END_REF] and winter wave intensity in several regions of the world, particularly in the southern hemisphere [START_REF] Semedo | Projection of Global Wave Climate Change toward the End of the Twenty-First Century[END_REF]. Hence, addressing the timescales of beach recovery to extreme storm winters, such as the 2013/14 winter, can provide a measure of coastal resilience in a changing climate.

Beach recovery from severe storms has been shown to spread over years to decades [START_REF] Morton | Stages and Durations of Post-Storm Beach Recovery, Southeastern Texas Coast[END_REF][START_REF] Houser | Sensitivity of post-hurricane beach and dune recovery to event frequency[END_REF]Castelle et al., 2017a). Since beach morphodynamics are often characterized by a significant seasonal signal [START_REF] Aubrey | Seasonal patterns of onshore/offshore sediment movement[END_REF][START_REF] Masselink | Seasonal changes in beach morphology along the sheltered coastline of Perth, Western Australia[END_REF]; [START_REF] Davidson | A behavioral template beach profile model for predicting seasonal to interannual shoreline evolution[END_REF], the long-term recovery signature is often hard to detect within the shorter-term fluctuations [START_REF] Thom | Behaviour of beach profiles during accretion and erosion dominated periods[END_REF][START_REF] Stephan | Long, Mid and Short-Term Evolution of Coastal Gravel Spits of Brittany, France[END_REF]. Therefore, high-frequency monitoring of beachmorphology over long time periods is crucial to understand better storm recovery [START_REF] Turner | A multi-decade dataset of monthly beach profile surveys and inshore wave forcing at Narrabeen, Australia[END_REF]. Unfortunately, such monitoring programmes are scarce, and the few available data sets have been used mostly to characterize extreme storm responses [START_REF] Scott | The extreme 2013/2014 winter storms: Beach recovery along the southwest coast of England[END_REF][START_REF] Barnard | Extreme oceanographic forcing and coastal response due to the 2015â€"2016 El Niño[END_REF], investigate the parameters controlling beach morphological changes [START_REF] Yates | Equilibrium shoreline response of a high wave energy beach[END_REF] and develop semi-empirical equilibrium models able to reproduce these morphological changes [START_REF] Davidson | A behavioral template beach profile model for predicting seasonal to interannual shoreline evolution[END_REF][START_REF] Yates | Equilibrium shoreline response: Observations and modeling[END_REF][START_REF] Splinter | A generalized equilibrium model for predicting daily to interannual shoreline response[END_REF]. However, ongoing field monitoring programmes in France and UK have recently shed some lights on the key mechanisms involved during post-storm recovery [START_REF] Scott | The extreme 2013/2014 winter storms: Beach recovery along the southwest coast of England[END_REF]Castelle et al., 2017a;[START_REF] Burvingt | Climate forcing of regionally-coherent extreme storm impact and recovery on embayed beaches[END_REF]. [START_REF] Scott | The extreme 2013/2014 winter storms: Beach recovery along the southwest coast of England[END_REF], investigated the morphological changes at three contrasting sites in SW England during the two years that followed the extreme 2013/14 winter. They found that the recovery mechanisms and timescales were highly dependent on the site characteristics, and that high-energy wave events were essential for the recovery of sediments. [START_REF] Burvingt | Climate forcing of regionally-coherent extreme storm impact and recovery on embayed beaches[END_REF], found that for a number of very similar beaches in SW England, recovery from the 2013/14 storm was regionally-coherent, multi-annual (>3 years), and mainly controlled by winter-wave conditions. Castelle et al. (2017a) investigated how the beach-dune system of an exposed site in SW France recovered from winter 2013/14 and found that only after 1.5 year the beach-dune system almost fully recovered to its pre-winter volume. These site-specific recovery rates highlight the need to conduct studies at broader scales, including different beaches, in order toinvestigate the key parameters that control the recovery timescales.

During the 2013/14 winter, a highly unusual sequence of extratropical storms crossed the North-East Atlantic region. This winter was the most energetic winter along the Atlantic coast of Europe since at least 1948 (Masselink et al., 2016a), and most of the west European coastline was severely impacted [START_REF] Castelle | Impact of the winter 2013-2014 series of severe Western Europe storms on a double-barred sandy coast: Beach and dune erosion and megacusp embayments[END_REF][START_REF] Blaise | Review of winter storms 2013-2014 on shoreline retreat dynamic on Brittany coast[END_REF][START_REF] Masselink | The extreme 2013/14 winter storms: hydrodynamic forcing and coastal response along the southwest coast of England[END_REF][START_REF] Autret | A comprehensive hydro-geomorphic study of cliff-top storm deposits on Banneg Island during winter 2013â€[END_REF]. Although winter waves are known to be well correlated with the North Atlantic Oscillation (NAO) index at high latitudes [START_REF] Bacon | A connection between mean wave height and atmospheric pressure gradient in the North Atlantic[END_REF][START_REF] Dodet | Wave climate variability in the North-East Atlantic Ocean over the last six decades[END_REF][START_REF] Bromirski | Wave power variability and trends across the North Atlantic influenced by decadal climate patterns[END_REF], this exceptional winter was not associated with a particularly high NAO. [START_REF] Castelle | A new climate index controlling winter wave activity along the Atlantic coast of Europe: the West Europe Pressure Anomaly[END_REF] computed a new climate index based on the sea level pressure gradient between Ireland and the Canary Islands: the West Europe Pressure Anomaly (WEPA). They showed that the 2013/14 winter was associated with the highest WEPA over 1948-2016, which reflects an intensified and southward shift of the sea level pressure difference between the Icelandic low and the Azores high, driving severe storms that funnel high-energy waves toward western Europe southward of 52 N. This high WEPA was linked to a recent increase in inter-annual variability, which was also observed in the winter-wave variability along western Europe [START_REF] Castelle | Increased Winter-Mean Wave Height, Variability, and Periodicity in the Northeast Atlantic Over 1949â€[END_REF].

In this paper, we investigate the post-2013/14 winter recovery of five beaches along the west coast of Europe; these are the same beaches for which the 2013/14 storm response was reported in Masselink et al. (2016a). Our objectives are threefold: 1) to obtain insight into the time scale of recovery for this extreme event for the different locations; 2) to explain the difference in observed recovery time scales by identifying the key factors involved; and 3) to determine extent to which extreme erosion and recovery processes can be modeled using present equilibrium models.

Methods

Wave Modeling

Two wave model hindcasts were used in this study. First, a large-scale and low-resolution model was used to characterize the wave climate in the North-East Atlantic, and more particularly the N-S differences in the wave forcing along the west coast of Europe. For this purpose, the spectral wave model WAVEWATCH III V4.18 (WW3, Tolman, 2014) was implemented on a 0.5 resolution grid covering the North Atlantic Ocean and forced with the 6-hourly wind fields of the NCEP/NCAR Global Reanalysis Project [START_REF] Kalnay | The NCEP/NCAR 40-Year Reanalysis Project[END_REF] from January 1948 to December 2017. Time series of significant wave heights ( ) were extracted at three deep-water locations (shown in Figure 1): north west of Ireland (10.0 W ; 56.0 N), in the Bay of Biscay (7.0 W ; 47.0 N), and west of Portugal (11.0 W ; 40.0 N). Details of the model setup and validation of the simulations with wave buoy observations can be found in itetmasselink-extreme-2016-1. Second, a smaller scale, high-resolution model was used to simulate the wave conditions close to the breaking point at each study site. Indeed, the offshore wave conditions simulated with the 0.5 model were not necessarily representative of nearshore wave conditions at some of the sheltered study locations. For this purpose we used a WW3 hindcast (1992-2017) implemented on an unstructured grid with a resolution increasing from 10 km offshore to 200 m in the coastal region extending from north of Spain to south of Ireland (Boudière et al., 2013). This model has been extensively validated with directional buoy and satellite altimeter and showed excellent skill, with correlation coefficients of morethan 0.94 and root-mean square errors less than 0.2 m for the whole set of validation points (Boudière et al., 2013; [START_REF] Ardhuin | Numerical Wave Modeling in Conditions with Strong Currents: Dissipation, Refraction, and Relative Wind[END_REF]. Model outputs were extracted for each study site at a distance less than 6 km from the coast, in water depths of 20-35 m. Seasonal means were computed for the winter (DJFM) and for the spring-summer-autumn (AMJJASON).

Study Sites and Beach Volumes

Five beaches along the Atlantic coast of Europe were surveyed on a monthly basis for more than 10 years. This data set represents one of the most complete series of beach profiles along western Europe. The location of the study beaches are shown in Figure 1, and the morphological characteristics of the study sites are given in Table 1. Additional information on the survey methods can be found in Masselink et al. (2016a). Since Slapton Sands displayed a strong alongshore variability in beach profile evolution, two representative beach profiles were analyzed separately, corresponding to the middle (SP10) and northern end (SP18) of the beach.

The extension of this data set to November 2017 was used to investigate the morphological recovery of the beaches four years after the exceptionnal 2013/14 winter. For this purpose, the beach volume above mean sea level (V) was computed for each site, with no upper limit set except at Perranporth where data was not collected for elevations higher than 3 m above MSL. Beach volume V, which therefore includes the dune system at Vougot, Porsmilin and Truc Vert, was assumed to provide an accurate and integrated measure of the beach system change (see left-hand panels of Fig. 4 in Masselink et al. (2016a)). Note that the upper beach of Porsmilin is topped by an artificial embankment that culminates at 6.2 m above MSL. Then, thebeach volume changes (

) were divided into four components: 1) beach volume change caused by the long-term trend computed over the period prior to the 2013/14 winter; 2) seasonal signal, computed from the detrended signal as the average annual difference between the maximum and minimum beach volume ( , where i is a yearly increment and N is the number of years in the time series); 3) 2013/14 winter response, computed as the difference in beach volume prior to and after the 2013/14 winter; and 4) post-2013/14 winter recovery, computed as the difference in beach volumes between April 2014 and November 2017. Note that the long-term trend and the seasonal contribution were only computed over the time period prior to the 2013/14 winter to ensure these signals were not affected by the 2013/14 winter storm response. Rates of beach volume changes (dV/dt) were computed for the winter and spring-summer-autumn from the observations closest in time to December 1 and April 1. When no observations were available within two weeks before or after these dates, the corresponding dV/dt was not computed. In the remaining of the paper, the percentage of recovery is computed as the beach volume change associated with the post-storm recovery relative to the beach volume change associated with the 2013/14 winter response, as defined above (components 3 and 4).

Beach Equilibrium Modeling

To assess whether an equilibrium-based model can be used to forecast beach recovery to an extreme storm event, the ShoreFor model [START_REF] Davidson | A simple equilibrium model for predicting shoreline change[END_REF] was applied. This semi-empirical model predicts beachface erosion when the wave conditions are more energetic than the equilibrium conditions (computed as a weighted average of past wave conditions) and vice-versa, and the magnitude of change is proportional to the incident wave power and degree of disequilibrium. The model has two free parameters that require calibration: a disequilibrium term and a linear trend term. The linear trend-term crudely accounts for all processes other than wave-driven cross-shore transport, including longshore sediment transport processes. The reader is referred to [START_REF] Davidson | A simple equilibrium model for predicting shoreline change[END_REF] for a full description of the model. For all sites, the model was calibrated wth the period of observations prior to the 2013/14 winter, and validated during the remaining period that includes the 2013/14 winter storm response and the subsequent 4-year recovery period. The model skill was assessed with the correlation coefficient (R) between observed (x) and simulated ( ) beach volumes, the root-mean-square error (RMSE), and the root-mean-square error normalized by the observed variance prior to the 2013/14 winter (NRMSE). Because records with a significant linear trend, possibly induced by longshore transport processes or other net source/sinks of sediments, sometimes showhigh model skill solely attributable to the linear trend component in the model, the model skill was also assessed using the Brier Skill Score (BSS), which allows comparison of the model residuals with a suitable baseline ( ), taken here as the linear trend component of the model. The BSS is computed as follows:

(1)

A positive BSS indicates an improvement relative to the baseline, and values greater than 0.0, 0.3, 0.6, 0.8 are typically described as 'poor', 'fair', 'good', 'excellent', respectively [START_REF] Van Rijn | The predictability of cross-shore bed evolution of sandy beaches at the time scale of storms and seasons using process-based Profile models[END_REF][START_REF] Sutherland | Evaluating the performance of morphological models[END_REF]. Note that this modelling approach does not resolve long-shore transport processes and is thus expected to show poor skills when applied to environments dominated by long-shore transport.

Results

Modeled Wave Conditions

The wave conditions simulated with the regional model over the period 2002-2017 for the north-west of Ireland, the Bay of Biscay, and west of Portugal are shown in Figure 1. A clear seasonal signal characterizes the three time series, with winter-mean much larger than spring-summer-autumn-mean (56 % greater on average, and up to 120 % greater locally). Moreover, the winter-mean values display strong inter-annual variability ( = 0.12 on average, where is the standard deviation, and is the long-term mean ), whereas the spring-summer-autumn-mean values display much lower inter-annual variability ( = 0.06). The consequence of these fluctuations is that, contrary to spring-summer-autumn means, the winter-mean may differ significantly from one year to another. For instance, the largest winter-mean in the Bay of Biscay and west of Portugal occurred during the 2013/14 winter, and they were approximately 35 % greater than the long-term mean winter . During the following winter, wave conditions were moderate in the Bay of Biscay and west of Portugal, but obtained their maximum north of Ireland. These trendswere inverted during the 2015/16 winter as the winter-mean was very large in the Bay of Biscay and west of Portugal, but moderate north of Ireland. The most recent 2016/17 winter was moderate in all three regions. This inter-annual variability of winter-mean was shown to be significantly correlated with the WEPA index southward of 52 N (Castelle et

Beach Recovery from the 2013/14 winter

Figure 2 shows the complete time series of beach volume changes for the six beach profiles (left-hand column), and the relative contributions of the long-term trend, seasonal signal, 2013/14 winter response and post-2013/14 winter recovery (right-hand column). Contrasting behaviors are observed. First, the most exposed sites, Perranporth and Truc Vert, suffered unprecedented erosion during the 2013/14 winter. Yet, after two years, Truc Vert had fully recovered, while Perranporth had only recovered 70 % after four years. The major difference in these recovery rates occurred during the year 2015. From early February to mid-December 2015, the beach volumes at Truc Vert increased steadily and the beach recovered more than 80 % of the sediments lost during the 2013/14 winter within a span of 10 months (see Castelle et al., 2017a, for details). At Perranporth, the beach was in a recovery phase for a shorter period of time -from late-March to November 2015 -regaining only 40 % of the sediments lost during the 2013/14 winter. This contrasting response can be directly related to the difference in wave conditions in January, March, November and December 2015 that were particularly stormy at Perranporth ( was 60 % higher than the annual mean at Perranporth and only 30 % higher at TrucVert). Porsmilin was also in its most eroded state after the 2013/14 winter, but after two years the beach had recovered by almost 80 %. This fast recovery was fostered by the relatively calm wave conditions during the 2014/15 winter that did not cause much erosion at this sheltered site. The beach volumes at Vougot are dominated by a decreasing long-term trend.

Although the coastal dune retreated by more than 5 mduring the 2013/14 winter, the sediment remained in the intertidal zone and the beach volume actually increased slightly. After four years, the dune had prograded back by approximately 3 m. At Slapton Sands, the central (SP10) and east (textitSP18) profiles showed opposite behaviors as a result of beach rotation processes. An additional factor that could explain the difference in recovery rates is the difference in tidal range. Large tidal range cause shorter time period during which sediment transport can contribute to accretion within the upper intertidal profile and subsequently longer morphological response times. However, no clear conclusion on this process was drawn from our data set, since both slow (Perranporth) and fast (Porsmilin) responses were observed on these macrotidal beaches.

To investigate the relationship between beach dynamics and incident wave conditions, the rates of beach volume changes (dV/dt) during the winter season and during the spring-summer-autumn season are compared to the respective seasonal wave energy anomalies, i.e., the deviation of the season-mean wave energy from the long-term (1992-2017) annual mean wave energy (Figure 3). Overall, dV/dt displays much greater variability during the winter season than during the rest of the year. At Perranporth, Pormsilin and Truc Vert, the winter-mean variability of dV/dt is clearly controlled by the wave conditions (0.58 < < 0.65). The near-zero intercept of the linear trends indicates that the beach profile is close to equilibrium when the winter-mean E is close to the long-term yearly mean . Although winter wave conditions are associated mostly with erosive conditions, low winter-mean Ecan cause beach accretion. For instance, during the 2009/10 winter, the wave conditions were particularly calm north of 50 N, due to a very low NAO and a modest WEPA, and the sand volume at Perranporth increased by 26 m /m. For the spring-summer-autumn season, correlations between dV/dt and wave energy anomalies are much lower and mostly insignificant at the 95 % level. One reason is that dV/dt cannot progressively increase when Etends towards zero; very low energy waves contribute less to onshore sediment transport than low to moderate energy waves, hence limiting recovery [START_REF] Hoefel | Wave-Induced Sediment Transport and Sandbar Migration[END_REF]Fernández-Mora et al., 2015). At Slapton Sands profiles SP10 and SP18, the winter-mean dV/dt is also strongly controlled by the wave conditions; however, the correlations have opposite signs as a result of beach rotation. [START_REF] Wiggins | The role of multi-decadal climate variability in controlling coastal dynamics: re-interpretation of the 'lost village of Hallsands[END_REF] showed very high correlations between beach volume changesand the directional wave power at Slapton Sands, and the insignificant correlations for the spring-summer-autumn season are probably because the beach changes were mostly controlled by the wave direction and not by the wave energy. At Vougot, thereis no correlation between dV/dt and the wave conditions. Indeed, the behavior of the beach-dune system is severely impacted by the presence of a jetty at the north-eastern end of the beach. Since its construction in 1974 the beachhas continually lost sediment, independent of the wave conditions [START_REF] Suanez | Monitoring of recent morphological changesof the dune of Vougot Beach (Brittany, France) using Differential GPS[END_REF].

Finally, the beach volume changes were compared with the results of the beach equilibrium model ShoreFor [START_REF] Davidson | A simple equilibrium model for predicting shoreline change[END_REF] to assess the amount of variance attributable to cross-shore sediment transport and to antecedent wave conditions. Theanalysis described in the previous paragraph treats each year independently, while ShoreFor accounts for antecedent wave conditions. It is therefore expected to explain more of the variability in the beach volume at the cross-shore dominated beaches through the disequilibrium term than a simple model based on a linear correlation between dV/dt and the mean wave height. Figure 4 shows the observed versus simulated beach volume changes using the ShoreFor model, as well as the error metrics textitRMSE, NRMSE, R, and BSS. Inspection of this figure reveals that Perranporth and Truc Vert have low NRMSE (<5 %) and 'excellent' BSS, indicating that ShoreFor is able to reproduce fairly well the storm response and subsequent recovery, and this variance is mostly induced by cross-shore processes. With a and a 'fair' BSS, ShoreFor results are moderate, and beach volume changes at Porsmilin can also be considered as dominated by cross-shore processes. Conversely, the negative BSS scores at Slapton SP18 and Vougot indicate that the model performs worse than predictions based on the long-term trend only. At Slapton SP10, both R and the BSS are relatively high; however, the very large NRMSE (270.4 %) reveals that some significant processes are ignored by the model. Hence, different modelling approaches must be applied at Vougot and Slapton Sands to simulate volume changes, including extreme storm response and post-storm recovery.

Discussion and Conclusions

The analysis of decadal morphological changes of beaches along the Atlantic coast of Europe exposed to a pronounced seasonal wave climate revealed that the dynamics of beaches are controlled by processes operating over a variety of time scales. In decreasing order these time scales are: long-term trends (decade), post-storm recovery (years), seasonal changes (months), and storm response (days). Total beach dynamics represent the sum of these components and for different beaches the relative contribution of each of these components varies significantly, making beach volume predictions challenging and site-specific. Moreover, beach recovery is conventionally thought to be a process that occurs during the calm summer months. However, although beachesdo recover during the spring-summer-autumn period at modest and relatively steady rates (not much inter-annual variability), it is the energetic winter conditions that primarily control the time it takes for beaches to recover from extreme erosion. Highly energetic winters stall or even reverse the recovery process, whereas calm winters continue the recovery process. Therefore, climate indices such as NAO and WEPA, which are known to explain a significant part of the inter-annual variability of winter wave conditions in the North-East Atlantic [START_REF] Dodet | Wave climate variability in the North-East Atlantic Ocean over the last six decades[END_REF][START_REF] Castelle | A new climate index controlling winter wave activity along the Atlantic coast of Europe: the West Europe Pressure Anomaly[END_REF], are well correlated with the recovery process. For instance, the most exposed sites Perranporth and Truc Vert required calm winter conditions to recover from the 2013/14 winter erosion, which correspond to negative values of WEPA. This was the case for the 2014/15 and 2016/17 winters (Figure 1), during which these beaches showed relatively small rates of volume changes (Figure 3). The recovery of these beaches could have been accelerated if the 2015/16 winter, which was characterized by a high WEPA value, had not caused severe erosion and slowed down the recovery process (Figure 3). At Slapton Sands, easterlies have been shown to foster beach recovery following storm erosion by (southwesterly) Atlantic storms, and these are promoted in this region by negative NAO values [START_REF] Wiggins | The role of multi-decadal climate variability in controlling coastal dynamics: re-interpretation of the 'lost village of Hallsands[END_REF]. The systematic positive NAO winters that followed the 2013/14 winter, and the prevailing southwesterly wave conditions, limited beach recovery at this site.

Predicting long-term beach morphological change is of great importance to coastal managers. While process-based morphodynamic modeling systems are valuable tools to simulate the morphological impact of single storm events (e.g. [START_REF] Mccall | Two-dimensional time dependent hurricane overwash and erosion modeling at Santa Rosa Island[END_REF][START_REF] Almeida | Storm overwash of a gravel barrier: Field measurements and XBeach-G modelling[END_REF], their computational cost prevents their application to multi-annual or even inter-annual morphological changes. In contrast, beach equilibrium models are computationaly cheap and can be applied for investigating long-term morphological changes (e.g. [START_REF] Yates | Equilibrium shoreline response of a high wave energy beach[END_REF][START_REF] Splinter | A generalized equilibrium model for predicting daily to interannual shoreline response[END_REF]. For cross-shore transport dominated sites, the ShoreFor model calibrated with topographic data prior to the winter 2013/14 and forced with nearshore wave conditions simulated with a high-resolution model showed significant skills in reproducing the strong erosion caused by the extreme 2013/14 winter and the recovery phase that followed at Truc Vert and Perranporth. Not surprisingly, ShoreFor shows poor skill at sites where longshore processes and resulting beach rotation signal dominate shoreline variability, i.e. at Vougot and Slapton. At Porsmilin, due to the low elevation of the artificial embankments at the top of the upper beach, a significant fraction of the sediment lost during the winter 2013/14 was deposited further inland during washover events. We believe this may explain why the model failed in reproducing accurately the volume changes during and after the winter. Semi-empirical models combining the equilibrium-based behaviour owing to variability in incident wave energy with longshore processes are scarce and still under development [START_REF] Vitousek | A model integrating longshore and cross-shore processes for predicting long-term shoreline response to climate change[END_REF][START_REF] Robinet | On a reduced-complexity shoreline model combining cross-shore and alongshore processes[END_REF]. Although out of scope of this study, the further development of these models will extend the domain of applicability of shoreline change models, making it possible to address coastal vulnerability and resilience in the context of climate change. [START_REF] Mentaschi | Global changes of extreme coastal wave energy fluxes triggered by intensified teleconnection patterns[END_REF] analyzed projection of extreme wave energy fluxes under a high emission scenario (Representative Concentration Pathways 8.5) and showed a significant increase in the 100-year return level of wave energy fluxes for the southern hemisphere and for some regions of the northern hemisphere. It is very likely that such changes in the wave climate will significantly impact beach morphodynamics at both event scales (storm response) and long-term scales (post-storm recovery), which will require accurate predictions forimplementing coastal adaptation strategies.

providing the wind field reanalysis and the wave hindcast, respectively. The topographic and wave data that were used for the analysis are provided as supplementary material of this paper.

  al., 2017b) and with the NAO index further north[START_REF] Bacon | A connection between mean wave height and atmospheric pressure gradient in the North Atlantic[END_REF][START_REF] Dodet | Wave climate variability in the North-East Atlantic Ocean over the last six decades[END_REF][START_REF] Bromirski | Wave power variability and trends across the North Atlantic influenced by decadal climate patterns[END_REF]. This dependence on NAO and WEPA indices is confirmed by our results, with the highest (respectively lowest) NAO during the 2014/15 (respectively 2009/10) winter correlating with the maximum (respectively minimum) north of Ireland for this winter, and the two highest WEPA during the 2013/14 and 2015/16 winters correlating with the maximum in the Bay of Biscay and west of Portugal for these winters. Correlation coefficients between the winter-mean and the NAO and WEPA indicesare shown on Figure1.

Figure 1 .

 1 Figure 1. (left) Location map of the Atlantic coast of Europe showing the offshore bathymetry (greyscale), virtual wave buoys (pink diamonds), beach study sites Perranporth (PP), Slapton Sands (SP), Vougot (VG), Porsmilin (PM), Truc Vert (TV) (red circles), and weather stations used to compute the NAO (green circles) and WEPA indices (yellow squares). The white contour line represents the 1000 m isobath. (right) Time series of NAO and WEPA indices (top panel), and time series of raw (grey line), 3-month filtered (black line), winter-mean (blue diamond) and spring-summer-autumn mean (red triangles) significant wave height at the virtual buoys 1, 2 and 3 (bottom 3 panels). The dashed rectangle indicates the 2013/14 winter and the 4-year recovery period that followed. Squared correlation coefficients ( ) between winter-mean and NAO and WEPA indices are provided for each virtual buoy.

Figure 2 .

 2 Figure 2. (left) Time series of beach volume at the five study sites (with two profiles shown for Slapton Sands), with the beach volume set to zero on December 1 2013. The dashed blue line represents the long-term trend over the period prior to the 2013/14winter, the red line represents the 2013/14 winter response, and the green line represents the recovery period. For Vougot and Truc Vert the evolution of the location of the dune foot (grey line) is also shown. (right) Absolute values of the volume change associated with the long-term trend (blue), seasonal variability (white), 2013/14 winter response (red), and recovery period (green).

Figure 3 .

 3 Figure 3. Beach volume changes during winter (circles) and summer-spring-autumn (squares) versus the wave energy anomaly (computed as the deviation of the season-mean wave energy from the long-term (1992â€"2017)annual mean wave energy), with colors indicating years. The squared correlation coefficients between beach volume changes and wave energy anomaly are given for winter (black) and spring-summer-autumn (grey). Linear regressions for winter (dashed light grey) and summer-spring-autumn (dasheddark grey) are plotted when the correlation is statistically significant at the 95 % level (in that case is written in bold).

Figure 4 .

 4 Figure 4. Comparison between ShoreFor model results and observations. Statistical errors are given for the validation period (post-2013/14 winter), and include the root-mean-square error (RMSE), the root-mean-square error normalized by the observed variance (NRMSE), the correlation coefficient (R) and the Brier Skill Score (BSS), with the long-term trend (dash black line) used as the baseline. Note that the calibration period at Perranporth (PP) does not include the four initial measurements, because their low sampling frequency induce error in the model calibration.

Table 1 .

 1 Summary of beach site characteristics. tan is the intertidal slope and MSR stands for mean spring tide range

	Name	Region		Exposure Hinterland	(mm) tan MSR (m)
	Perranporth Cornwall,		W Exposed	Dunes	0.35 0.015 4.5
		UK
	Slapton Devon,		SE	Lagoon	2-8	0.1 4.3
		UK	Semi-sheltered
	Vougout Brittany,		NNW	Dunes	0.2-0.3 0.03 8.5
		France		Semi-exposed
	Porsmilin Brittany,		S	Seawall,	0.32 0.05 5.7
		France		Semi-shelterd	Marsh
	Truc Vert Aquitaine,	W Exposed	Dunes	0.4	0.025 3.9
		France

Acknowledgments

This work was supported financially by the French Naval Oceanographic and Hydrographic Department (Shom) through the research program PROTEVS (12CR6), and by the "Laboratoire d'Excellence" LabexMER (ANR-10-LABX-19-01) program. BC was funded by the SONO (ANR-17-CE01-0014) and CHIPO (ANR-14-ASTR-0004) projects of the Agence Nationale de la Recherche (ANR) and by the AST "Evenements extremes" of the Observatoire Aquitain des Sciences de l'Univers (OASU). GM and TS were funded by the NERC BLUE-coast project (NE/N015525/1). The French sites Vougot, Porsmilin and Truc Vert are monitoring sites of the Service National d'Observation (SNO) DYNALIT labelled by CNRS-INSU and part of the French Reseach Infrastructure ILICO. The authorswould like to thank the very many colleagues, postdocs, field technicians and PhD students involved in the beach profile data collection. We acknowledge NCEP/NCAR and Ifremer institutes for


Marvin K Nakayama

Bruno Tuffin

Sufficient Conditions for a Central Limit Theorem to Assess the Error of Randomized Quasi-Monte Carlo Methods

  L'archive ouverte pluridisciplinaire

INTRODUCTION

Analyzing a stochastic model frequently involves computing the mean performance µ. Often, µ can be expressed as an integral of a function h over an s-dimensional unit hypercube [0, 1] s for some fixed s ≥ 1. Such integrals for s > 1 are typically analytically intractable, leading to the use of numerical methods, including simulation. As these techniques incur error, we should give a measure of the error.

Monte Carlo (MC) estimates µ via random sampling [START_REF] Asmussen | Stochastic Simulation: Algorithms and Analysis[END_REF]. Repeatedly feeding independent and identically distributed (i.i.d.) uniformly distributed random vectors on [0, 1] s into integrand h produces i.i.d. outputs, which are averaged to yield the MC estimator. The method affords simple error estimation through a confidence interval (CI). Based on a central limit theorem (CLT), a CI uses the sample variance to provide a computable (probabilistic) measure of the MC error. But as the sample size n (i.e., number of evaluations of h) grows, the CI and the MC estimator's root-mean-square error (RMSE) shrink at a slow rate n -1/2 ; adding another digit of precision requires a 100-fold increase in n.

To obtain a more efficient estimator, quasi-Monte Carlo (QMC) replaces the i.i.d. uniforms driving the MC method with n deterministic points from a low-discrepancy sequence (e.g., a lattice or digital net), designed to more evenly fill [0, 1] s than a typical random sample; see [START_REF] Niederreiter | Random Number Generation and Quasi-Monte Carlo Methods[END_REF] and [START_REF] Lemieux | Monte Carlo and Quasi-Monte Carlo Sampling[END_REF]. When the integrand h has bounded Hardy-Krause variation, the Koksma-Hlawka inequality (e.g., Section 2.2 of [START_REF] Niederreiter | Random Number Generation and Quasi-Monte Carlo Methods[END_REF] shows that the QMC error decreases as O(n -1 (ln n) s ) as n → ∞, better than the rate at which MC's RMSE shrinks. While theoretically useful, the Koksma-Hlawka inequality has limited practical value as its bound is not easily computed and is often quite loose.

Randomized QMC (RQMC) suggests a way to obtain a computable error bound: randomize the QMC points r ≥ 2 i.i.d. times and build a CI from the sample variance of the resulting r i.i.d. estimators; e.g., see [START_REF] Tuffin | Randomization of Quasi-Monte Carlo Methods for Error Estimation: Survey and Normal Approximation[END_REF], Section 6.2 of [START_REF] Lemieux | Monte Carlo and Quasi-Monte Carlo Sampling[END_REF][START_REF] L'ecuyer | Randomized Quasi-Monte Carlo: An Introduction for Practitioners[END_REF]. For a given (large) computation budget of n integrand evaluations, we specify the number m of points used from each randomized sequence so that mr ≈ n. To choose such an allocation (m, r), a common rule of thumb recommends taking r small (e.g., 10 ≤ r ≤ 30) so that m is correspondingly large to benefit from QMC's superior convergence rate.

The RQMC CI's validity implicitly assumes that the RQMC estimator obeys a Gaussian CLT. When m → ∞ and r is fixed, [START_REF] Loh | On the Asymptotic Distribution of Scrambled Net Quadrature[END_REF] establishes a CLT that covers only a computationally prohibitive form of RQMC, limiting its practical use. But more generally, a Gaussian limit is not guaranteed; e.g., randomly shifting a lattice leads to non-normal limits (as m → ∞ for fixed r) [START_REF] L'ecuyer | On the Distribution of Integration Error by Randomly-Shifted Lattice Rules[END_REF]. Thus, while intuitively appealing, the CI lacks rigorous theoretical justification for most RQMC methods.

Our paper addresses these shortcomings. We provide sufficient conditions on both h and (m, r) that ensure the RQMC estimator obeys a CLT, as well as an asymptotically valid CI (AVCI). We focus on the setting where both m, r → ∞ since a Gaussian limit may not hold for fixed r. We will show tradeoffs in our restrictions on h and (m, r): more stringent limitations on h lead to looser constraints on (m, r). But in all cases, the RQMC RMSE shrinks faster than for the corresponding MC estimator with sample size n = mr.

The rest of the paper unfolds as follows. Section 2 builds our study's basic framework. We present general conditions that yield a CLT and AVCI in Sections 3 and 4, respectively. Section 5 provides simpler sufficient conditions for a CLT or AVCI, and gives graphical comparisons of the alternative restrictions. Concluding remarks are in Section 6. All formal proofs appear in [START_REF] Nakayama | Sufficient Conditions for Central Limit Theorems and Confidence Intervals for Randomized Quasi-Monte Carlo Methods[END_REF].

NOTATION AND FRAMEWORK

For an integrand h : [0, 1] s → ℜ on the unit hypercube of fixed dimension s ≥ 1, the goal is to compute

µ = [0,1] s h(u) du = E[h(U)],
where random vector U ∼ U [0, 1] s with U [0, 1] s denoting a uniform distribution on [0, 1] s , ∼ means "is distributed as", and E represents the expectation operator. We can think of h as a (complicated) simulation program that transforms s i.i.d. 1-dimensional uniform random numbers into observations from specified input distributions, which are then used to produce an output of the random performance of a stochastic system, so µ is its mean. We next explain how to apply MC, QMC, and RQMC to estimate µ.

Monte Carlo

With MC, we generate n i.i.d. copies U 1 ,U 2 , . . . ,U n of U ∼ U [0, 1] s , and compute µ MC n = ∑ n i=1 h(U i )/n as the MC estimator of µ. Let ψ 2 ≡ Var[h(U)], with Var[•] the variance operator, and assume that 0 < ψ 2 < ∞. The MC estimator is unbiased (i.e., E[ µ MC n ] = µ), as are all the estimators of µ that we consider, so

RMSE µ MC n = ψ √ n . (1)
The MC estimator obeys a Gaussian CLT Billingsley 1995, Theorem 27.1), where ⇒ denotes convergence in distribution, and N (a, b 2 ) is a normal random variable with mean a and variance b 2 . Let

√ n[ µ MC n -µ]/ψ ⇒ N (0, 1) as n → ∞ (
ψ 2 n = ∑ i=1 [h(U i ) -µ MC n ] 2 /(n -1) be the sample variance of the i.i.d. h(U i ).
For a desired confidence level 0 < γ < 1, we can exploit the CLT to construct an approximate γ-level CI for µ as

I MC n,γ ≡ [ µ MC n ± z γ ψ n / √ n],
where the critical point z γ satisfies Φ(z γ ) = 1 -(1γ)/2 and Φ is the N (0, 1) cumulative distribution function (CDF). Providing a probabilistic measure of the MC estimator's error, I MC n,γ is an AVCI in the sense that lim n→∞ P(µ ∈ I MC n,γ ) = γ (Asmussen and Glynn 2007, p. 71).

Quasi-Monte Carlo

QMC replaces MC's i.i.d. uniforms with carefully placed deterministic points from a low-discrepancy sequence Ξ = (ξ i ) i≥1 , such as a digital net (e.g., a Sobol' sequence) or lattice; see, e.g., Chapters 3-5 of [START_REF] Niederreiter | Random Number Generation and Quasi-Monte Carlo Methods[END_REF]. Using the first n points from Ξ leads to QMC approximating µ by µ Q n = ∑ n i=1 h(ξ i )/n. We can bound the error | µ Q n -µ| via the Koksma-Hlawka inequality (Niederreiter 1992, Section 2.2):

| µ Q n -µ| ≤ V HK (h)D * n (Ξ) (2)
for all n > 1, where D * n (Ξ) is the star-discrepancy of the first n points in Ξ, and V HK (h) is the Hardy-Krause variation of the integrand h. In (2), V HK (h) ≥ 0 quantifies the "roughness" of h, and D * n (Ξ) ∈ [0, 1] measures the "nonuniformity" of Ξ. Low-discrepancy sequences often have

D * n (Ξ) = O(n -1 (ln n) s ), as n → ∞, (3)
where

f (n) = O(g(n)) (resp., f (n) = Θ(g(n)))
as n → ∞ for functions f and g means that there exist positive constants a 0 , a 1 , and

n 0 such that | f (n)| ≤ a 1 |g(n)| (resp., a 0 |g(n)| ≤ | f (n)| ≤ a 1 |g(n)|) for all n ≥ n 0 . Thus, if V HK (h) < ∞,
(2) and (3) imply that the QMC error shrinks as O(n -1 (ln n) s ) as n → ∞, better than the Θ(n -1/2 ) rate at which MC's RMSE decreases. While theoretically useful, the bound in (2) has limited practical value as it is not easily computed and is often quite loose. There are other related error bounds (e.g., [START_REF] Hickernell | A Generalized Discrepancy and Quadrature Error Bound[END_REF][START_REF] Hickernell | Error Analysis for Quasi-Monte Carlo Methods[END_REF][START_REF] Lemieux | Quasi-Random Number Techniques[END_REF][START_REF] Niederreiter | Random Number Generation and Quasi-Monte Carlo Methods[END_REF]), but all suffer from the same issues. RQMC employs such a randomization r ≥ 1 i.i.d. times, and for each j = 1, . . . , r, let Ξ j ≡ (U i, j ) i≥1 be the jth randomized sequence. Given a computation budget of n integrand evaluations, we specify the number m of points to use from each Ξ j so that mr ≈ n, leading to the RQMC estimator of µ as

Randomized

µ RQ m,r = 1 r r ∑ j=1 X j , where X j = 1 m m ∑ i=1 h(U i, j ). (4)
The X j , j = 1, . . . , r, are i.i.d., and let σ 2 m,r = ∑ r j=1 (X jµ RQ m,r ) 2 /(r -1) be their sample variance when r ≥ 2. We then arrive at a possible γ-level CI I RQ m,r,γ ≡ [ µ RQ m,r ± z γ σ m,r ] for µ. The literature includes several methods to construct Ξ , including scrambled digital nets [START_REF] Owen | Randomly Permuted (t,m,s)-Nets and (t,s)-Sequences[END_REF][START_REF] Owen | Monte Carlo Variance of Scrambled Net Quadrature[END_REF]) and digital shifts (L'Ecuyer 2018). To simplify the discussion, we describe only one approach: random shifts [START_REF] Cranley | Randomization of Number Theoretic Methods for Multiple Integration[END_REF]. Here, randomization j generates a single U j ∼ U [0, 1] s and adds it (modulo 1) to each point in Ξ, so the ith point in the jth randomized sequence Ξ j is U i, j = U j + ξ i , where x is the modulo-1 operator applied to each coordinate of x ∈ ℜ s . The U j across randomizations j = 1, 2, . . . , r, are independent. It is easy to show that each U i, j ∼ U [0, 1] s , so µ RQ m,r and each X j are unbiased estimators of µ. But for each randomization j, the sequence Ξ j has dependent points because they all share the same uniform U j .

With random shifts, each randomized sequence Ξ j satisfies (Tuffin 1997, Theorem 2)

D * m (Ξ j ) ≤ 4 s D * m (Ξ). (5) Thus, if V HK (h) < ∞, the estimator X j in (4) from a single randomization of m points satisfies RMSE[X j ] = O(m -1 (log m) s ) as m → ∞, an improvement over the Θ(m -1/2
) rate in (1) for MC using the same number m of integrand evaluations. Even faster convergence rates can be achieved for special classes of functions and specific sequences Ξ called lattice rules [START_REF] Tuffin | Variance Reduction Order Using Good Lattice Points in Monte Carlo Methods[END_REF][START_REF] L'ecuyer | Variance Reduction via Lattice Rules[END_REF].

Although intuitively appealing, the CI I RQ m,r,γ in general lacks theoretical justification, as it implicitly relies on µ RQ m,r obeying a Gaussian CLT. For m → ∞ with r ≥ 1 fixed, [START_REF] Loh | On the Asymptotic Distribution of Scrambled Net Quadrature[END_REF] establishes an RQMC CLT that covers solely the case of fully nested scrambling of a digital net, which is computationally prohibitive, limiting its adoption by practitioners. For random shifts of a lattice, the RQMC estimator µ RQ m,r may not obey a Gaussian CLT as m → ∞ for fixed r ≥ 1, as shown by L' Ecuyer et al. (2010). Indeed, they prove that for r = 1, the limiting error distribution has simple non-Gaussian forms for dimension s = 1, and s > 1 generally leads to non-Gaussian limits with no such easy characterizations, so the same holds for any fixed r ≥ 1. Thus, we see the need for general Gaussian CLTs for RQMC, which is our aim.

Assumptions and Preliminary Results

We want to study the asymptotic behavior of the RQMC estimator in (4) as the computation budget n for the number of integrand evaluations grows large. To do this, we take the number m ≡ m n ≥ 1 of points from the randomized sequence and the number r ≡ r n ≥ 1 of randomizations to be functions of n satisfying Assumption 1.A m n r n ≤ n for each n ≥ 1, with m n → ∞, r n → ∞, and m n r n /n → 1 as n → ∞.

Under Assumption 1.A, the RQMC estimator in (4) becomes

µ RQ m n ,r n = 1 r n r n ∑ j=1 X n, j , where X n, j = 1 m n m n ∑ i=1 h(U i, j ), (6)
so X n, j averages h on the first m n points of the jth randomized sequence. Our goal is to provide conditions on h and (m n , r n ) that yield (as n → ∞) a Gaussian CLT (Section 3) or AVCI (Section 4). Other papers (e.g., [START_REF] Glynn | Limit Theorems for the Method of Replication[END_REF][START_REF] Damerdji | Strong Consistency of the Variance Estimator in Steady-State Simulation Output Analysis[END_REF]) adopt frameworks akin to Assumption 1.A to study MC methods for analyzing steady-state behavior via multiple replications or batching. Assumption 1.A requires r n → ∞ because otherwise, the limiting error distribution may not be Gaussian, as noted at the end of Section 2.3. We simplify the discussion by further having m n → ∞ in Assumption 1.A, but this is not necessary; [START_REF] Nakayama | Sufficient Conditions for Central Limit Theorems and Confidence Intervals for Randomized Quasi-Monte Carlo Methods[END_REF] also analyze the special case that m n ≡ m 0 for a fixed m 0 ≥ 1. Section 5 will adopt the following specialization of Assumption 1.A. Assumption 1.B m n = n c and r n = n 1-c with c ∈ (0, 1).

We should define, e.g., m n = n c and r n = n 1-c ( • is the floor function) so that m n and r n are integers, but for simplicity, we ignore this technicality. Section 5 will determine constraints on h and c ∈ (0, 1) that secure a CLT or AVCI, and in each case, the optimal such c that minimizes the rate at which RMSE[ µ RQ m n ,r n ] shrinks as n → ∞. Also, we will examine the tradeoffs in the conditions on h and c. For a randomized sequence Ξ constructed from scrambling or a digital shift of a digital net, or for a randomly shifted lattice rule, the randomization preserves the partitioning structure of the original sequence Ξ: a randomly shifted lattice is still a lattice, and scrambling or digitally shifting a digital net retains the original sequence's finer-grain properties [START_REF] Owen | Randomly Permuted (t,m,s)-Nets and (t,s)-Sequences[END_REF][START_REF] Owen | Monte Carlo Variance of Scrambled Net Quadrature[END_REF][START_REF] L'ecuyer | Randomized Quasi-Monte Carlo: An Introduction for Practitioners[END_REF]. Moreover, these Ξ obey similar discrepancy bounds as Ξ. Specifically, consider any low-discrepancy sequence Ξ for which (3) holds, so there exists some constant 0 < w 0 < ∞ such that D * m (Ξ) ≤ w 0 m -1 (ln m) s for all m > 1. Then its random shift Ξ satisfies D * m (Ξ ) ≤ w 0 m -1 (ln m) s with w 0 = 4 s w 0 by (5), and scrambling or digital shifting digital nets yields analogous bounds. Thus, all of these randomizations fulfill the following assumption, which we use to analyze RQMC estimators when V HK (h) < ∞. Assumption 2 For the RQMC method used, there exists a constant 0 < w 0 < ∞ such that each randomized sequence Ξ satisfies D * m (Ξ ) ≤ w 0 m -1 (ln m) s for all m > 1, where w 0 depends on the RQMC method but not on the randomization's realization (e.g., of U ∼ U [0, 1] s in a random shift).

We often will further impose one of the following conditions on the integrand h. The conditions are presented in order of decreasing strength (see Proposition 1 below), and Section 5 will show that this leads to corresponding tradeoffs in our conditions on (m n , r n ) to ensure a CLT or AVCI. Assumption 3.A The integrand h is of bounded Hardy-Krause variation, i.e., V HK

(h) < ∞. Assumption 3.B The integrand h is bounded; i.e., |h(u)| ≤ t 0 for all u ∈ [0, 1] s for some constant t 0 < ∞. Assumption 3.C There exists b > 0 such that E |h(U) -µ| 2+b < ∞, where U ∼ U [0, 1] s .
Limiting the roughness of h over [0, 1] s , Assumption 3.A imposes substantial restrictions; it does not hold, e.g., in dimension s ≥ 2 when h is an indicator function (so µ is a probability) with discontinuities not lining up with the coordinate axes [START_REF] Owen | A Strong Law of Large Numbers for Scrambled Net Integration[END_REF]. In contrast, Assumption 3.C constrains the heaviness of the tails of the distribution of h(U). Proposition 1 Assumption 3.A is strictly stronger than Assumption 3.B, itself strictly stronger than Assumption 3.C.

Using different conditions on h, we next derive two bounds on absolute central moments of the estimator X n,1 in (6) from a single randomization. The first lemma, for V HK (h) < ∞ (Assumption 3.A), follows from Theorem 2 of [START_REF] Tuffin | Variance Reduction Applied to Product-Form Multi-Class Queuing Networks[END_REF]; the second applies Minkowski's inequality (Billingsley 1995, eq. (5.40)) when Assumption 3.C holds for 2 + b replaced by q ≥ 1. Lemma 1 Under Assumptions 1.A, 2, and 3.A, for any q > 0 and for all n such that m n > 1,

η n,q ≡ E [|X n,1 -µ| q ] ≤ E (V HK (h)D * m n (Ξ )) q ≤ w 0 V HK (h)(ln m n ) s m n q < ∞. ( 7
)
Lemma 2 Under Assumption 1.A, for any

q ≥ 1, if E[|h(U) -µ| q ] < ∞ for U ∼ U [0, 1] s , then η n,q ≤ E[|h(U) -µ| q ] for every n.
For a single randomization of m points, RQMC typically has

σ m ≡ (Var[∑ m i=1 h(U i,1 )/m]) 1/2 = O(m -α ) as m → ∞ with α > 1/2 (e.g., see (7) when V HK (h) < ∞). This improves on MC's RMSE convergence rate, which satisfies RMSE[ µ MC m ] = Var[ µ MC m ] = ψm -1/2
by (1). Assume the following limit exists:

α * = -lim m→∞ ln(σ m ) ln(m) , (8)
so α * is the constant such that σ m decreases, as m → ∞, at a rate (ignoring leading coefficients and lower-order terms) strictly faster than m -α * +ε and strictly slower than m -α * -ε for every ε > 0; i.e., σ m = o(m -α * +ε ) and

σ m = ω(m -α * -ε ) as m → ∞ for any ε > 0, where f (m) = o(g(m)) as m → ∞ means that f (m)/g(m) → 0 as m → ∞, and f (m) = ω(g(m)) as m → ∞ means that f (m)/g(m) → ∞ as n → ∞. By (7), α * ≥ 1 when V HK (h) < ∞, (9)
as in Assumption 3.A, and more generally, as is typical of RQMC, we assume that

α * > 1 2 . ( 10
)
The value of α * depends on the particular integrand h and the RQMC method applied, but not on how (m n , r n ) or c are specified in Assumptions 1.A and 1.B. For any randomized sequence

Ξ = (U i ) i≥1 , let E m (Ξ ) = 1 m ∑ m i=1 h(U i )
µ be the error of the estimator based on the first m points. Then define its exponential rate as

α(Ξ ) = -lim m→∞ ln(|E m (Ξ )|) ln(m) , (11)
assuming the limit always exists, and define the worst-case rate among all randomizations as

α = inf Ξ α(Ξ ). (12)
Thus, for any randomized sequence Ξ , we see that

|E m (Ξ )| = o(m -α+ε ) and |E m (Ξ )| = ω(m -α-ε ) as m → ∞
for each ε > 0, and α ≤ α * always holds. Assumption 4 The convergence rate exponent α of the worst-case error among all randomizations is the same as the standard deviation rate exponent α * ; i.e., α = α * . Corollary 1 in Section 5.1 will later show that when Assumption 4 holds, the RQMC estimator µ RQ m n ,r n will obey a CLT for (m n , r n ) = (n c , n 1-c ) for any c ∈ (0, 1) (Assumption 1.B). However, establishing Assumption 4 in practice may be difficult, so much of our paper focuses on providing other more verifiable sufficient conditions that secure a CLT.

For a given total number n of integrand evaluations, RQMC papers often suggest choosing m n as large as possible to gain from the fast convergence rate of QMC. But as noted at the end of Section 2.3, we still should specify r n big enough to ensure a Gaussian CLT. Under Assumptions 1.B, 2 and 3.A, (7) implies

RMSE[ µ RQ m n ,r n ] ≤ [w 0 V HK (h)(ln m n ) s /m n ] √ r n = Θ (c ln n) s n (1+c)/2
, so larger c leads to faster convergence. While c = 1 minimizes the RMSE bound, a Gaussian CLT may not be guaranteed, as noted earlier at the end of Section 2.3.

GENERAL CONDITIONS FOR A CENTRAL LIMIT THEOREM

We next study limiting properties of µ RQ m n ,r n in (6) as n → ∞. The estimator averages X n,1 , X n,2 , . . . , X n,r n , but their distribution changes with n, complicating the asymptotic analysis. To develop a theoretical framework for handling this under Assumption 1.A, note that (X n, j ) n=1,2,...; j=1,2,...,r n forms a triangular array (Billingsley 1995, p. 359), also called a double array. In a triangular array, the r n variables within a row n are independent, but there may be dependence across rows. While the general formulation allows for the CDFs of the r n variables within a row n to differ, RQMC actually has X n,1 , X n,2 , . . . , X n,r n are i.i.d., each with some distribution F n ,

where F n may change with n, as occurs in (6). To preclude trivialities, we impose another assumption, without which the exact result is eventually always returned by the RQMC estimator. Assumption 5 σ 2 m n ≡ Var[X n,1 ] > 0 for all sufficiently large n. The Lindeberg and Lyapounov CLTs (Billingsley 1995, Theorems 27.2 and 27.3) apply for the RQMC structure in (13). To set them up, write the variance of the sum of the r n random variables in (13) as s 2 n ≡ r n σ 2 m n . Denote the CDF of X n, jµ by G n , which does not depend on j by (13). Note that σ 2 m n = y∈ℜ y 2 dG n (y), and let τ 2 n (t) = |y|>ts n y 2 dG n (y) for t > 0. Then the RQMC estimator in (6) satisfies the following. Theorem 1 If Assumptions 1.A and 5 hold and also the Lindeberg condition

τ 2 n (t) σ 2 m n → 0, as n → ∞, ∀t > 0, (14)
then the RQMC estimator in (6) satisfies the CLT

µ RQ m n ,r n -µ σ m n / √ r n ⇒ N (0, 1), as n → ∞. (15)
Also, ( 14) holds if, for some b > 0, E |X n,1 -µ| 2+b < ∞ for each n and the Lyapounov condition holds:

E |X n,1 -µ| 2+b r b/2 n σ 2+b m n → 0, as n → ∞. ( 16
)
The Lindeberg condition (14) constrains the tail behavior of G n . In the general setting of independent but not identically distributed X n, j , 1 ≤ j ≤ r n , the analogous version of (14) (Billingsley 1995, eq. (27.8)) ensures that the contribution of any single X n, j to their sum's variance s 2 n is negligible for large n. We can show (Billingsley 1995, p. 361) that ( 14) is even necessary for the CLT (15) when ( 13) holds. Imposing restrictions on moments rather than tail properties, the Lyapounov condition ( 16) can sometimes be easier to apply than the Lindeberg condition ( 14). Section 5.1 will obtain sufficient conditions for ( 14) or ( 16) under Assumption 1.B to secure Theorem 1 for each of our Assumptions 3.A-3.C on the integrand h.

To build a CI from the CLT (15), suppose that r n ≥ 2, which Assumption 1.A ensures for all n large enough. As X n, j , j = 1, 2, . . . , r n , are i.i.d. by ( 13), their sample variance

σ 2 m n ,r n = ∑ r n j=1 X n, j -µ RQ m n ,r n 2 /(r n -1)
provides an unbiased estimator of σ 2 m n = Var[X n,1 ]. For a given desired confidence level γ ∈ (0, 1), we get

I RQ m n ,r n ,γ ≡ µ RQ m n ,r n ± z γ σ m n ,r n / √ r n ( 17
)
as the RQMC CI for µ. The next result imposes conditions guaranteeing that I RQ m n ,r n ,γ is an AVCI in the sense that (20) below holds. Theorem 2 Suppose that Assumptions 1.A and 5 hold. Also, suppose that E (X n,1µ) 4 < ∞ and that

E (X n,1 -µ) 4 r n σ 4 m n → 0, as n → ∞. ( 18
) Then µ RQ m n ,r n -µ σ m n ,r n / √ r n ⇒ N (0, 1), as n → ∞, (19)
and

lim n→∞ P(µ ∈ I RQ m n ,r n ,γ ) = γ. (20)
As ( 18) is the same as ( 16) for b = 2, Theorem 1 implies CLT ( 15), which is expressed in terms of the exact σ m n . But the left side of ( 19) instead uses the estimator σ m n ,r n . Theorem 2's conditions further ensure σ m n ,r n /σ m n ⇒ 1 as n → ∞, so Slutsky's theorem (Billingsley 1995, p. 340) verifies (19), securing AVCI (20). Section 5.1 will provide sufficient conditions under Assumption 1.B that yield Theorem 2 for two of our conditions on the integrand h (Assumptions 3.A and 3.C).

ANALYSIS WHEN

(m n , r n ) = (n c , n 1-c ) (ASSUMPTION 1.B)
Assumption 1.B specializes Assumption 1.A by taking (m n , r n ) = (n c , n 1-c ) for some c ∈ (0, 1). We next will determine the values of c that imply CLT (15) through Theorem 1 or that guarantee AVCI (20) via Theorem 2. For those c, we then find the ones leading to RMSE[ µ RQ m n ,r n ] shrinking fastest as n → ∞. Under Assumption 1.B, we have that m n = n c with c ∈ (0, 1), so (8) implies that

σ m n = ω n -cα * -ε and σ m n = o n -cα * +ε as n → ∞, for all ε > 0. ( 21
)
Taking ε > 0 arbitrarily small in ( 21) leads to σ m n ≈ Θ(n -cα * ) as n → ∞. Thus, a combination of r n = n 1-c with ( 6) and ( 13) yields, for any c ∈ (0, 1),

RMSE µ RQ m n ,r n = σ m n √ r n ≈ Θ n -v(α * ,c) as n → ∞, where v(α * , c) ≡ c α * - 1 2 + 1 2 . ( 22
)
Our assumption (10) guarantees that v(α * , c) > 1/2 for c ∈ (0, 1). Hence, the convergence rate of RQMC's RMSE for any c in Assumption 1.B is better than RMSE

[ µ MC n ] = Θ(n -v MC ) as n → ∞ for MC, where v MC ≡ 1 2 (23)
by ( 1). For any fixed α * satisfying (10), v(α * , c) strictly increases in c by ( 22), so RQMC's RMSE shrinks faster for larger c. We thus want to determine how large c can be and still ensure CLT (15) or AVCI (20).

Section 5.1 will provide various corollaries of the CLT and AVCI theorem in Sections 3 and 4. Each such Corollary k will result in restricting c as

c < c k (α * ) (24)
for some 0 < c k (α * ) ≤ 1 depending on the particular Corollary k. As we will see, most of the c k (α * ) are strictly decreasing in α * . Thus, as α * increases, ( 24) further restricts the choices of c, thereby reducing the maximum allowable number of QMC points and increasing the minimum number of randomizations because (m n , r n ) = (n c , n 1-c ). But by ( 8), larger α * corresponds to a better convergence rate for the estimator based on a single randomization, so in some sense, securing a CLT or AVCI often entails hampering better QMC performance. Because ( 22) implies that larger c leads to RMSE shrinking at a faster rate, the "optimal" c that maximizes the rate subject to the constraint ( 24) is c = c k (α * )δ for infinitesimally small δ > 0. Accordingly, an analysis akin to the arguments applied to achieve (22) arrives at the optimal approximate rate:

RMSE[ µ RQ m n ,r n ] ≈ Θ n -v k (α * ) as n → ∞, (25)
where, for each Corollary k (and k ) in Section 5.1, the exponent v k (α * ) appears below in (26).

Proposition 2 Under Assumption 1.B and (10), the optimal approximate RMSE rate exponent in ( 25) is

v k (α * ) ≡ c k (α * ) α * - 1 2 + 1 2 > v MC (26) for v MC in (23), so RQMC outdoes MC. If c k (α * ) = 1 in (24), then v k (α * ) = α * . Also, for any k and k , v k (α * ) > v k (α * ) if and only if c k (α * ) > c k (α * ). (27)
When c k (α * ) = 1, (24) becomes the weakest possible constraint satisfying Assumption 1.B. In this case, Proposition 2 implies that v k (α * ) = α * , so the RMSE of the multiple-randomization RQMC estimator µ RQ m n ,r n decreases at about the same rate as for a single randomization with full length m = n.

The next subsection will specialize c k (α * ) in ( 24) and v k (α * ) in ( 26) for various corollaries. Section 5.2 will compare the resulting values graphically.

Corollaries of Theorems 1 and 2

We first provide a corollary of Theorem 1 based on Assumption 4, which imposes constraints on both the integrand and RQMC sequence. Corollary 1 Suppose that Assumptions 1.B, 4, and 5 hold. If c < 1 ≡ c 1 (α * ), then the Lindeberg condition (14) and CLT (15) hold. Moreover, ( 25) and ( 26

) have v k (α * ) = v 1 (α * ) ≡ α * .
Under Assumption 4 (the worst-case error and the standard deviation decrease at the same exponential rate), Corollary 1 secures a CLT for (m n , r n ) = (n 1-ε , n ε ) with ε > 0 as small as we wish. Thus, although r n → ∞ is needed, choosing ε > 0 small allows taking a large number m n = n 1-ε of points from the low-discrepancy sequence, which enables exploiting QMC's superior convergence rates.

As establishing Assumption 4 may be difficult in practice, we next provide other conditions that are more readily verifiable to ensure CLT (15). We give corollaries corresponding to each of our restrictions on the integrand h in Assumptions 3.A-3.C, which are in decreasing order of strength (Proposition 1). We first specialize (16) of Theorem 1 to establish a CLT when V HK (h) < ∞, which enables using Lemma 1. Corollary 2 Suppose that Assumptions 1.B, 2, 3.A (V HK (h) < ∞), and 5 hold. Also, suppose that for some constant λ ∈ (0, 1),

c < 1 -λ 2α * -1 -λ ≡ c 2 (α * ), (28)
where λ may be chosen arbitrarily small. Then the Lyapounov condition ( 16) and CLT (15) hold. Moreover, for each α * ≥ 1, as in (9), c 2 (α * ) satisfies 0 < c 2 (α * ) ≤ 1, and ( 25) and ( 26)

have v k (α * ) = v 2 (α * ) with v 2 (α * ) ≡ 2α * -1 -λ α * 2α * -1 -λ , where 1 2 < v 2 (α * ) ≤ 1. ( 29
)
Corollary 2 allows taking λ ∈ (0, 1) as arbitrarily small in ( 28) and ( 29), so c 2 (α * ) ≈ 1/(2α * -1) and v 2 (α * ) ≈ 1. The next corollary of Theorem 1 exploits ( 14) to yield a CLT when the integrand h is bounded. Corollary 3 Suppose that Assumptions 1.B, 3.B (h is bounded), and 5 hold. If

c < 1 2α * + 1 ≡ c 3 (α * ),
then the Lindeberg condition ( 14) and CLT (15) hold. Moreover, for each α * > 1/2, as in (10), c 3 (α * ) satisfies 0 < c 3 (α * ) < 1/2, and ( 25) and ( 26)

have v k (α * ) = v 3 (α * ) with v 3 (α * ) ≡ 2α * 2α * + 1 , where 1 2 < v 3 (α * ) < 1.
The following corollary of Theorem 1 imposes a moment condition on h(U) (Assumption 3.C) to apply Lemma 2 to ( 16) to obtain a CLT, in contrast to requiring V HK (h) < ∞, as in Corollary 2. Corollary 4 Suppose that Assumptions 1.B, 3.C (finite absolute central moment of order 2 + b for some b > 0), and 5 hold. If

c < 1 2α * (1 + 2 b ) + 1 ≡ c 4 (α * , b),
then the Lyapounov condition ( 16) and CLT (15) hold. Moreover, for each b > 0 and α * > 1/2, as in (10), c 4 (α * , b) satisfies 0 < c 4 (α * , b) < 1/2, and ( 25) and ( 26)

have v k (α * ) = v 4 (α * , b) with v 4 (α * , b) ≡ 2α * (1 + 1 b ) 2α * (1 + 2 b ) + 1
, where 1 2 < v 4 (α * , b) < 1.

For I RQ m n ,r n ,γ in (17) to be AVCI (20), Theorem 2 assumes that (18) holds, which yields the CLTs in ( 15) and ( 19). We next consider conditions that enable verifying AVCI. Corollary 5 Suppose that Assumptions 1.B, 2, 3.A (V HK (h) < ∞), and 5 hold. If

c < 1 4α * -3 ≡ c 5 (α * ),
then the CLT (19) and AVCI (20) hold. Moreover, for each α * ≥ 1, as in (9), c 5 (α * ) satisfies 0 < c 5 (α * ) ≤ 1, and ( 25) and ( 26)

have v k (α * ) = v 5 (α * ) with v 5 (α * ) ≡ 3α * -2 4α * -3 , where 3 4 < v 5 (α * ) ≤ 1.
While Corollary 5 requires V HK (h) < ∞, we next ensure AVCI instead through a moment condition. Corollary 6 Suppose that Assumptions 1.B and 5 hold, as well as Assumption 3.C (finite absolute central moment of order 2 + b) for b = 2. If

c < 1 4α * + 1 ≡ c 6 (α * ),
then the CLT (19) and AVCI (20) hold. Moreover, for each α * > 1/2, as in (10), c 6 (α * ) satisfies 0 < c 6 (α * ) < 1/2, and ( 25) and ( 26)

have v k (α * ) = v 6 (α * ) with v 6 (α * ) ≡ c 6 (α * ) α * - 1 2 + 1 2 = 3α * 4α * + 1
, where 1 2 < v 6 (α * ) < 3 4 . In the left panel of Figure 2, the upper bounds c k (α * ) on c decrease as α * grows, so ensuring CLT (15) or AVCI (20) for larger α * requires putting more effort on the MC part (i.e., r n = n 1-c grows as c decreases) and correspondingly less on the QMC (i.e., m n = n c shrinks as c gets smaller). By (26), the tradeoff could potentially harm the rate exponent v k (α * ) governing how quickly the RQMC estimator's optimal RMSE decreases, but this does not occur for most k. The one exception is v 5 (α * ) for the AVCI Corollary 5 when V HK (h) < ∞, which we explain by examining the corresponding c 5 (α * ) in the left panel of Figure 2. While c 5 (α * ) starts off at α * = 1 very high, it quickly drops off, so m n must decrease rapidly as α * grows to secure AVCI when V HK (h) < ∞, leading to less benefit from the QMC. Even so, we have that v 5 (α * ) > v 6 (α * ) for all α * , so the optimal rate exponent when establishing AVCI is better for V HK (h) < ∞ than through the moment condition of Corollary 6.

(α * = 1.25) c * = 1 c 2 (α * ) c 3 (α * ) c 4 (α * , b) c 5 (α * ) c 6 (α * )

CONCLUDING REMARKS

We presented conditions that ensure the RQMC estimator of a mean µ obeys a Gaussian CLT or guarantee AVCI. We also examined the tradeoffs in the restrictions. While our paper gave sufficient conditions, we are currently looking into relaxing the requirements. Other current work includes devising procedures to estimate the upper bounds c k (α * ) in (24) and Section 5.1, which will allow practitioners to apply our theoretical results. We are further investigating analogous theory for biased estimators, as for quantiles.

Figure 1 :

 1 Figure 1: Plots of the upper bounds c k (α * ) in (24) of c in Assumption 1.B for different Corollaries k from Section 5.1. The plots display the c k (α * ) as functions of b for different fixed values of α * . The upper left panel does not include c 2 (α * ) and c 5 (α * ) as these require V HK (h) < ∞, which then implies α * ≥ 1 by (9). The plots show that stronger conditions on h correspond to loosening constraints on c.

  Quasi-Monte CarloRQMC applies i.i.d. randomizations of the QMC sequence Ξ to produce i.i.d. estimators of µ, and builds an approximate CI via their sample variance. A randomization creates from Ξ another sequence Ξ ≡ (U i ) i≥1 that retains the low-discrepancy properties of Ξ. Each U i ∼ U [0, 1] s , but the points in Ξ are dependent.

ACKNOWLEDGMENTS

The authors would like to thank Art Owen for pointing out references, used in the proof of Proposition 1, showing that if the integrand h has bounded Hardy-Krause variation, then h is bounded.