Spectral quantification of nonlinear behaviour of the nearshore seabed and correlations with potential forcings at Duck, N.C., U.S.A
SSpectral quantification of nonlinear behaviour of thenearshore seabed and correlations with potentialforcings at Duck, N.C., U.S.A.
Magar, V. a, ∗ , Lefranc, M. b , Hoyle, R. B. c , Reeve, D. E. a a Coastal Engineering Research Group, Marine Institute, School of Marine Science andEngineering, University of Plymouth, Drake Circus, Plymouth PL48AA. Tel: +44(0)1752 586137, Fax: +44 (0)1752 232638 b Laboratoire de Physique des Lasers, Atomes, Mol´ecules (PhLAM) and Centre d’ ´Etudeset de Recherches Lasers et Applications (CERLA), UFR de Physique, Bat. P5,Universit´e des Sciences et Technologies de Lille, F-59655 Villeneuve d’Ascq CEDEX,France c Department of Mathematics, University of Surrey, Guildford, Surrey, GU2 7XH, UK
Abstract
Local bathymetric quasi-periodic patterns of oscillation are identified frommonthly profile surveys taken at two shore-perpendicular transects at theUSACE field research facility in Duck, North Carolina, USA, spanning 24.5years and covering the swash and surf zones. The chosen transects are thetwo furthest (north and south) from the pier located at the study site. Re-search at Duck has traditionally focused on one or more of these transects asthe effects of the pier are least at these locations. The patterns are identifiedusing singular spectrum analysis (SSA). Possible correlations with potentialforcing mechanisms are discussed by 1) doing an SSA with same parametersettings to independently identify the quasi-periodic cycles embedded withinthree potentially linked sequences: monthly wave heights (MWH), monthlymean water levels (MWL) and the large scale atmospheric index known asthe North Atlantic Oscillation (NAO) and 2) comparing the patterns withinMWH, MWL and NAO to the local bathymetric patterns. The results agreewell with previous patterns identified using wavelets and confirm the highly ∗ Corresponding author
Email addresses: [email protected] (Magar, V.), [email protected] (Lefranc, M.),
[email protected] (Hoyle, R. B.), [email protected] (Reeve, D. E.) a r X i v : . [ n li n . C D ] J u l onstationary behaviour of beach levels at Duck; the discussion of potentialcorrelations with hydrodynamic and atmospheric phenomena is a new con-tribution. The study is then extended to all measured bathymetric profiles,covering an area of 1100m (alongshore) by 440m (cross-shore), to 1) anal-yse linear correlations between the bathymetry and the potential forcingsusing multivariate empirical orthogonal functions (MEOF) and linear cor-relation analysis and 2) identify which collective quasi-periodic bathymetricpatterns are correlated with those within MWH, MWL or NAO, based ona (nonlinear) multichannel singular spectrum analysis (MSSA). MEOF andlinear correlation analysis are equivalent, with spatial distributions of thelinear correlation coefficient between bathymetry and the potential forcingbeing the same as those of the relevant EOF. Regions of high correlations aresupported by what would be expected from knowledge about local physicalprocesses. At the two transects furthest from the pier, regions with highcorrelation and anti-correlation broadly agree with the SSA findings. TheMEOF/correlation analysis was used as a preliminary step towards the iden-tification of collective spatio-temporal patterns between all bathymetric pro-files and the potential forcings. This final analysis, based on MSSA, showedthat no collective interannual patterns of oscillations are present throughoutthe bathymetry and the three potential forcings - which was to be expectedas SSA showed the highly localised nature of the interannual bathymetricpatterns. Also, annual and semi-annual cycles within the bathymetry arestrongly correlated to the monthly wave height, in agreement with the SSAfindings. Other collective intra-annual cycles besides the semi-annual wereidentified; they were all correlated to the NAO. Key words:
Sandy beaches, Beach features, Coastal morphology, Spectralanalysis, Oscillations, Duck
1. Introduction
Understanding the evolution of beaches in the long term (at yearly todecadal time scales) poses an important challenge to researchers. This isbecause it is unclear how the mechanisms forcing the morphology, such ashydrodynamic and sediment transport processes, affect beaches at such timescales (Southgate et al., 2003; Dodd et al., 2003), and also because datasetsfrom which yearly to decadal patterns of behaviour may be extracted areat present scarce. Because of the threat of climate change, improved un-2erstanding and better predictive tools are crucial, both for researchers andend-users alike.The term ’coastal morphodynamics’ has been coined to describe the in-vestigation of changes in the shape (or morphology) of coastal features suchas beaches, sandbanks and so on. There are several modelling approachesadopted in coastal morphodynamics, for instance sequential models, equilib-rium models, equation-based or data-based models. In long-term morpho-dynamical modelling data-based approaches tend to be preferred becauseof the difficulties in assessing how the different processes affect the dynam-ics in the long-term, as was mentioned above. Equilibrium-based models ofsandbar systems assume sandbar migration towards a wave height dependentequilibrium location (see, e.g. Plant et al., 1999). Sequential models eitheranalyse single cross-shore profiles or incorporate longshore variability; in theformer the aim is an analysis of cross-shore sediment transport, in particularthe characterisation of cyclic patterns, while in the latter the aim consists ofidentifying morphological beach states that would characterise morphologicalsequences and thus help understand and classify nearshore bars in terms oftheir shape and dynamics (Wright and Short, 1984; Lippmann and Holman,1990). An important aspect of sandbar dynamics is the observed interactionsbetween sandbars when several bars are present. At Duck, such interactionshave been studied by Lippmann et al. (1993), who observed significant non-linear behaviour of inner bars when outer bars were present.There are two possible signal selection criteria used to analyse nearshoremorphodynamics. One is to analyse spatial patterns, for instance the evolu-tion of the sandbar crests and troughs or the evolution of the shoreline, andthe other is to find patterns optimizing a statistical measure, for instancepatterns identified with empirical orthogonal functions (EOFs), CanonicalCorrelation Analysis (CCA) or wavelet analysis. The spatial pattern ap-proach has been extremely popular and has contributed significantly to theunderstanding of sandbar dynamics (see Lippmann and Holman, 1990; Lipp-mann et al., 1993; Wijnberg and Terwindt, 1995; Plant et al., 1999; Wijnbergand Kroon, 2002; R´o˙zy´nski, 2003; Kuriyama et al., 2008; Pape and Ruessink,2008, and references therein). These studies focus on descriptions of bar mi-gration, bar amplitude evolution and development of bar asymmetry. Us-ing an optimal statistical measure as a pattern identification criterion hasalso been very popular; these methods are particularly useful for temporalor spatio-temporal pattern analysis. For instance, EOFs identify the linearmodes of maximal variance, while CCA leads to linear modes that maximize3he correlations between two time series. CCA has been used recently at Duckto analyse correlations between wave climate and bathymetric evolution atDuck, with reasonable success (Larson and Kraus, 1994; Horrillo-Caraballoand Reeve, 2008). Reeve et al. (2008), on the other hand, used wavelets toanalyse dominant temporal patterns of behaviour at a cross-shore transectin Duck. In some cases, it has been possible to link the spatial patternsto the patterns optimizing a statistical measure (see, e.g., Lippmann et al.,1993; Wijnberg and Terwindt, 1995); in others the interpretation of the pat-terns optimizing a statistical measure may be less obvious, for instance forthe 1-2 yearly temporal patterns embedded in the seabed dynamics at Duck,identified by Reeve et al. (2008), in which case additional analyses may benecessary to assess their origin.We postulate that different bathymetric regions have different nonlinearbehaviour depending on their cross-shore and longshore location, and thatseveral of the observed patterns are well correlated with patterns embed-ded within potential forcing mechanisms such as the monthly wave patterns,monthly sea level patterns or the NAO. Hence, nonlinear patterns of be-haviour of bed levels at Duck will be quantified using singular spectrum anal-ysis (SSA) and spectral density analysis SDA to characterise quasi-periodicoscillations. Possible correlations between these patterns and those embed-ded within large-scale phenomena that might potentially be at the originof the observed behaviour will be dicussed, and correlations between themwill be confirmed by further multivariate EOF (MEOF) and multichannelsingular spectrum analysis (MSSA) studies.SSA and MSSA are two data-based, nonlinear methods that optimizemodes of maximal variance. The methods are applied to time-lagged copiesof bathymetric (and hydrodynamic or atmospheric index) time series andwere chosen because they are good at identifying underlying temporal andspatio-temporal patterns of oscillation even when the data is short and noisy(Vautard and Ghil, 1989), as is the case for most coastal morphodynamicsdatasets. Vautard and Ghil (1989) have shown, for instance, that SSA iscapable of describing the dynamics of paleoclimatic oscillations obtained frommarine sediment core measurements that had as few as 184 data points, aswell as the statistically significant number of degrees of freedom. There areother measures of dimension aside from the number of degrees of freedom ofa dynamical system, such as the correlation dimension, which is a measureof the dimensionality of the space occupied by a set of points. However,Grassberger (1986) and Vautard and Ghil (1989) have showed that reliable4orrelation dimension computations for the same core measurements were notpossible due to the data’s shortness and noisiness. Indeed, it is well knownthat in general the number of data points needed for accurate computationof the correlation dimension would be much larger, more likely in the orderof a thousand or larger (Grassberger and Procaccia, 1983; Ruelle, 1990).When analysing temporal patterns, however, it is important as well that thefrequency of sampling is significantly larger than the frequency of the studiedpattern, regardless of the size of the dataset (Gilmore and Lefranc, 2002). Forour purposes SSA (and by extension MSSA) is a perfectly suitable techniquesince our interest focuses on the quasi-oscillations underlying the dynamicsand these oscillations are well captured by such techniques.SSA has some similarities with empirical orthogonal functions (EOF),first applied to analysing beach variation by Winant et al. (1975). Both SSAand EOF, as well as MSSA and Complex EOF (CEOF), have very recentlybeen exploited to analyse long-term sea level changes (Papadopoulos andTsimplis, 2006), as well as behavioural patterns in coastal features such asshorelines (R´o˙zy´nski et al., 2001; Southgate et al., 2003; R´o˙zy´nski, 2005;Miller and Dean, 2006) and nearshore sandbars (Wijnberg and Terwindt,1995; Rattan et al., 2005). SSA may be used together with spectral densityanalysis (SDA) to characterise the frequencies identified with SSA. However,one major benefit of SSA over spectral analysis is that, just as with waveletanalysis, SSA is capable of capturing nonstationarity in patterns (Li et al.,2005; Reeve et al., 2007).Southgate et al. (2003) used SSA to identify the long-term trend of theshoreline at Duck, between 1980 and 1993. The authors analysed the compo-nent of the signal resolved by the first three eigenvalues. Their reconstructionshowed no evidence of quasi-periodic oscillations of the shoreline. Rather, itled to the identification of an accretional trend from 1981 to 1990, followedby significant erosion from 1990 onwards. This contrasted with their findingsat Ogata beach, Japan, where a clear linear trend of erosion, with a 5-yearlycycle superimposed, could be identified.As R´o˙zy´nski et al. (2001) indicate, the specification of what is a signaland what is noise is crucial in SSA. These authors have used this techniqueto analyse shoreline evolution at Lubiatowo, Poland, between 1983 and 1999.Three long-term patterns were identified, with periods of 9, 16 or 32 years,at different positions along the coastline. Interannual cycles of 2-3 and 3-4 years were also found. Significant variability occurred at some positions,however. This study also pointed to the possible existence of chaotic be-5aviour as resolved by the second eigenvalue, which the authors linked toa response to extreme events. An association between the rhythmicity ofthe bed level variations and sand bar propagation was also postulated, sup-ported by the existence of high correlations between shoreline positions andthe inner bar crests (Pruszak et al., 1999). The authors also explain somefindings in terms of self-organisation. While the implementation of EOFis relatively straightforward, with SSA intermediate computations and testsare necessary to obtain a set of uncorrelated reconstructed components andavoid misinterpretations. However, SSA is superior to EOF because it canidentify nonlinear, temporal patterns underlying the dynamics, while EOF isa linear method which cannot identify nonlinear patterns at all. MSSA fol-lows the same principles as SSA, but is performed on a vector that containsvariations at different positions rather that at a single location. Thus, MSSAis more likely to identify collective patterns of behaviour. MSSA is some-times called multivariate extended empirical orthogonal function (MEEOF)analysis (Mote et al., 2000).MSSA has been applied to Lubiatowo, to identify the collective patternsof behaviour of the shoreline (R´o˙zy´nski, 2005). However, such patterns willnecessarily appear as well in the SSA at individual positions, but possibly noteverywhere as the dominant patterns of behaviour. At Lubiatowo, three col-lective quasi-periodic oscillations were identified, with periods of 7 to 8 years,20 years and several decades. The authors commented that the 7 to 8-yearlycycle might be forced by the North Atlantic Oscillation (NAO), which hasan important influence in the Baltic Sea. Very recently, (R´o˙zy´nski, 2010) ex-tended his analysis and found supporting evidence that the NAO may affectthe bathymetry through a gentle coupling between the NAO, the hydrody-namics and the seabed.As indicated by Falques et al. (2000), the spatial regularity of nearshorepatterns implies that nearshore dynamics may be understood and, most im-portantly, described using simple physical mechanisms. That we need to lookfor mechanisms acting in the long term is particularly important. Within thesurf zone and at the shore, the occurrence of spatially rhythmic patterns hasbeen attributed to the presence of low-frequency waves, such as infragrav-ity waves (Bowen and Guza, 1978; Bowen, 1980; Holman and Bowen, 1982).However, infragravity waves are responsible for the short term response ofsome beaches, not the interannual behaviour. At such large time scales thewave-related parameter that is of importance is the monthly averaged waveheight. Note that we are not implying there ought to be a simple relationship6etween the waves and the bathymetry, whose response to forcings has beenshown to be nonlinear in other studies (Falques et al., 2000, and referencestherein). This nonlinearity could be reflected in the difference between am-plitudes and phases, as well as in the characteristics of the regime changesobserved in the corresponding time series. In the case of longshore paral-lel bars, the bar crests and troughs may cause redistribution of wave energy,variations in the wave breaking point and wave transformation processes (i.e.wave refraction, reflection and diffraction) which in turn, produce a radiationstress which is no longer in equilibrium with the set-up and set-down, hencecreating a steady circulation (Falques et al., 2000). Cellular automata mod-els have been developed to analyse free evolution of the seabed (Coco et al.,1999), while nonlinear morphodynamic evolution models have permitted theanalysis of difference in bed response depending on an initial perturbation(see, e.g., Klein and Schuttelaars, 2006). Aside from the monthly averagedwave heights and atmospheric phenomena such as NAO, it is likely that thetides and storm surges also have an effect on long term morphodynamics.Tides are fully predictable shallow water waves generated by gravitationalforces between the Earth, the Sun and the Moon. Storm surges changethe water depth through pressure forces, wind stresses and wave set-up andset-down, and this in turn changes the water depth at the bar crests andtroughs, and hence the radiation stresses mentioned above. Given that tidesand storm surges may be combined in one single parameter, namely the meanwater level, some of the patterns observed in the bathymetric evolution maypotentially be correlated with patterns found within a water level time series.In Sect. 2 below the case study data is presented and the SSA and SDAtechniques are explained in more detail, together with the MEOF and theMSSA. The results are presented in Sect. 3 and they are then discussedin Sect 4. Sensitivity to window length changes is discussed in Sect. 4.1.Comparisons with previous investigations are presented in Sect. 4.2. InSect. 4.3, SSA and SDA of the monthly averaged wave heights at a nearshorewave gauge are performed and the patterns are compared to those of thebathymetry, together with similar analyses for the North Atlantic Oscillation(taken as an average over the whole area of study) to illustrate the possiblecorrelations with large-scale phenomena. Sect. 4.4 discusses linear correla-tions between the full set of bathymetric measurements and the monthlymean wave heights, monthly mean sea level and monthly NAO. Coherentspatio-temporal correlations identified via MSSA are also discussed. Finally,Sect. 5 contains a summary and some concluding remarks. Websites are re-7erred to as data or software sources, when appropriate, in accordance withcopyright requirements.
2. Summary of data and techniques y and the cross-shore coordinate is x . The baseline is a shore parallelline with its origin at the South-East corner of the FRF property. On-goingquality control and early publication of the surveys for those transects led toa significant number of studies based on them (see introduction). Moreover,the bathymetric evolution at these profiles is the least affected by the pres-ence of the pier(Nicholls et al., 1998). However, important differences in thebehaviour north and south of the pier, in particular regarding the sandbarsystem, have been observed (Shand et al., 1999).Duck is an example of a net offshore bar migration (NOM) system, ora system with interannual cyclic sandbar behaviour (Shand et al., 1999;Ruessink and Terwindt, 2000). NOM is a cyclic phenomenon consisting ofthree stages (Ruessink and Kroon, 1994): bar generation near the shoreline(stage 1); systematic offshore migration of the bar across the surf zone dur-ing high energy periods (stage 2) and; disappearance of the bar in the outersurf zone (stage 3). Sandbar formation and migration may also be triggered8y previous profile geometries, with the degeneration of an outer bar oftenleading to formation of a new bar near the shoreline (Ruessink and Terwindt,2000). NOM systems are often alongshore coherent, but at Duck this coher-ence is somewhat broken by the presence of the pier. Profiles at Duck varyfrom unbarred to triple barred, with the most common profile configurationbeing a double bar system consisting of a narrow inner bar and a broad outerbar (Howd and Birkemeier, 1987; Lee and Birkemeier, 1993; Lee et al., 1998).Three parameters may characterise a NOM cycle at each stage: the av-erage cross-shore distance over which bars migrate; the average duration ofthe bar migration and; the average return period of migration cycles. Theseparameters are site specific and in Duck they vary North and South of thepier. Shand et al. (1999, Fig. 7) summarise site parameter characteristicsfor a number of NOM systems including Duck North (DN) and Duck South(DS). While the total distance and the total duration of bar migration arerelatively similar at DN and DS, the average return period is significantly dif-ferent. The total distance travelled by the bars is 289.6 and 289.1 metres andtheir total duration (adding all three stages) is of 4.4 years and 4.1 years atDS and DN, respectively. However, for stage 2 for instance, the most activestage, the return period at DS is 3.2 years while at DN it is 6.8 years. Thesand is transported back to the shoreline during fairweather periods througha variety of mechanisms: in the outer surf zone transport is driven by waveasymmetry while in the inner surf zone it also occurs through migrating bod-ies such as transverse bars or bar bifurcates generated during bar splittingprocesses (Holman and Sallenger, 1986; Holman and Lippmann, 1987; Leeet al., 1998; Shand, 2007). Since the wave climate has a very strong seasonalpattern with storms occurring mostly during the winter, annual patterns areexpected to appear within the bathymetry.In this study, the two outermost transects, i.e. transects 58 (Duck North)and 190 (Duck South), were first chosen to analyse the local quasi-oscillationsembedded in the most active part of the surf zone. This reduced the cross-shore range to 390 m, between x = 100 and x = 490; this range coincidedwith that analysed by Southgate and M¨oller (2000) and includes most of theNOM width, that is, the region of sandbar migration (see Lee et al., 1998,Fig. 3). However the study area was extended to x = 80 m (to includethe shoreline processes in more detail) and to x = 520 m (to check the typeand strength of the correlations further offshore) for the linear correlationand MSSA analyses.] A spatial and temporal Akima interpolation (Akima,1970) was performed on the data to obtain depth values at constant 10 m9ntervals and at even timesteps of 30 days. These time and space interpolationintervals were also chosen by Southgate and M¨oller (2000) and prevent over-interpolation of the data. This procedure led to a 299-point time seriesat every bathymetric position, spanning from July 1981 to January 2006.Cross-shore bathymetric anomalies (with respect to the time-average at eachposition along the transect) are shown in Figs (a,i) and (b,i) within Fig. 1for transects 58 and 190, respectively. Two complementary methods were used to identify local patterns of be-haviour: Singular Spectrum Analysis (SSA) and Spectral Density Analysis(SDA). Both methods are described in detail in Appendix A.1 so here weonly present a brief summary of SSA. Essentially, SSA is a nonlinear eigen-value problem for a M-dimensional sequence Z n = ( z n , z n +1 , . . . , z n + M − ), for n = 1 , , . . . , N − M + 1, of time-lagged copies of the bed level z ( t ), at times t = iτ s , i = 1 , , . . . , N , N being the length of the time series, τ s the samplinginterval and M < N is the window length. Time scales of the dynamics thatcan be reconstructed from this time series are between τ s and the windowtime span. Reconstructed components (RCs), linked to the eigenvalues of Z n , have the same time span as the original time series and isolate one ormore or the quasi-oscillatory patterns (depending on the number of eigen-values used for the reconstruction). In this way, the signal was separatedinto a trend and a detrended signal using a window length, M , of 4 years.The trend corresponds to linear variations as well as patterns with periodsabove 3 to 5 years and is the RC of the first few eigenvalues of the SSAwith M = 4. The trend can be analysed further with an SSA with windowlength of, say, 9 years, to identify long-period oscillations (LPOs) embeddedwithin the trend signal. The detrended signal, on the other hand, containsshort-period oscillations (SPOs), of periodicities of less than 5 years. Again,the periodicities of the SPOs can be determined with a further SSA on thedetrended part of the signal. Note the overlap between LPOs and SPOs;this is due to the chosen window length of 4 years for the initial SSA. Therobustness of the observations to small changes in the window length will bediscussed in more detail further on.Patterns embedded within monthly mean wave heights (MWH), monthlymean water levels (MWL) and the monthly North Atlantic Oscillation (NAO)were also extracted following the same procedure. The details of the hydro-dynamic and atmospheric data are discussed in Appendix A.1. The patterns10btained were compared to those within the bathymetry to identify locationswhere similar patterns are observed. This led to an evaluation of correlationsbetween bathymetric and hydrodynamic/atmospheric patterns. Note thatSSA would give a nonlinear measure of the correlations. Linear correlationswere found as a preliminary step for in the analysis of spatio-temporal pat-terns of behaviour; the correlation analysis and the methods involved in thespatio-temporal study are introduced in the next section. Spatio-temporal collective bedlevel patterns were then analysed, not onlyon Transects 58 and 190 as before but also for all bathymetric surveys be-tween them. A linear correlation analysis was first performed, by computingthe correlation coefficient matrix, r , between the bedlevel time series andthat of each of the potential forcings. This led to three different correla-tion coefficient matrices, one for the correlations of the bathymetry withthe NAO, one for correlations with the MWH and one for correlations withthe MWL, respectively. The values in r can be used to produce correlationcoefficient contours to identify locations where the correlation between thebedlevels and the forcing is strongest. This analysis is intrinsically linkedwith the MEOF analysis carried out as a preliminary step for MSSA for allthe bathymetry transects, because the correlation coefficient is equivalent tothe covariance matrix, C , normalised by the square root of the product ofthe covariance at the associated diagonal elements: r ( i, j ) = C ( i, j ) (cid:112) C ( i, i ) C ( j, j ) , and MEOF is based on the covariance of the multivariate matrix. The MEOFmethod is described in detail in Sect. A.2 within the Appendix. MEOFpermits reduction in the number of channels from 1128 (one channel foreach bathymetric location and one channel for each of the three potentialforcings evaluated at a single location - see the Appendix for details) to threechannels corresponding to the MEOF PCs in directions of maximal varianceof the potential forcings (see Table A.3), with EOF1 corresponding to NAO,EOF2 to the mean wave heights (MWH) and EOF3 to the mean water levels(MWL). This reduction from 1128 to 3 channels is called a prefiltering. Thesystem of 3 channels is then analysed with MSSA to identify which collective,11edlevel spatio-temporal patterns may be correlated with NAO, MWH orMWL. The MSSA technique is described in Sec A.3 in the Appendix.
3. Results
The trends, shown in Figs. 1(a,ii) and 1(b,ii) for transects 58 and 190,respectively, capture well the dominant patterns observed in Figs. 1(a,i)and 1(b,i), which seem to correspond to the extrema in the bathymetry andtherefore are likely to correspond to the sandbars extrema. This is supportedby the fact that some of the travelling characteristics of the sandbars havealso been captured, because the extrema as well as surrounding bathymetriccontours seem to travel offshore as time progresses.For Transect 58, the patterns extracted resolve approximately 50% of thevariance in regions where the extrema are well resolved. At other locations,such as between 130 and 220 m offshore, the variance resolved by the trendis slightly smaller than 50%. However, these results indicate that the trendand the detrended signal are equally important because they resolve an equalproportion of the total variance, and hence both should be analysed further.For Transect 190, the results are very similar, i.e. the variance resolvedby the trend and the detrended signals is about the same, though with thetrend accounting for slightly less than 50% of the variance in this case, and soboth signals need to be analysed. As with transect 58, the detrended signal,which contains the higher frequency oscillations (with periods generally below4 years), seems to resolve slightly more of the variance between 130 and 220 moffshore.The periods of the underlying cycle were computed all along transects 58and 190; Table 2 shows the results at selected locations along the profiles.The uncertainties in the periods were estimated from figures of the powerspectrum as well as from the uncertainty of the power values at each fre-quency. An example of a power spectrum is shown in Fig. 3, correspondingto position x = 160 m at transect 58. If the uncertainty interval of the powervalue extends to or beyond the peak power value, then the uncertainty of theposition of the peak also extends to the frequency at that power value. Long period oscillations (LPOs) and short period oscillations (SPOs) -with periods above and below 4 years, respectively - were then extracted12rom the trend and the detrended time series. The RCs obtained whencombining all the LPOs are shown in Figs. 2(a) and 2(c) for transects 58 and190, respectively while the SPOs are shown in Figs. 2(b) and 2(d).In relation to the LPOs, it is clear that these patterns have an importantcontribution below 350 to 400 m offshore, because the variance resolved bythese patterns - shown in the upper plot of Figs. 2(a) and 2(c) - typicallytakes large values there. At both transects the signal contains roughly threeregions where quasi-oscillations are important. At transect 58 these regionsextend between x = 130 and 140 m, x = 210 and 260 m and x = 340and 350 m; at transect 190, they extend between x = 130 and 180 m, x =230 and 240 m and x = 290 to 370 m. Hence there is agreement in thenumber of regions but their extent differs. The main periods of the quasi-oscillatory patterns are summarised in Table 2. The table shows there is alot of variability in the periods of the patterns, however, between 220 and370 m offshore the periods of the patterns are close to those found withinthe NAO.The SPOs, in contrast, may have an important contribution throughoutthe region under consideration, as the variance they resolve - see upper plotsin Figs. 2(b) and 2(d) - is usually around 50%. Based on the information inFig. 2(b) and Table 2, at transect 58 the signal contains a very clear yearlycycle in the region between x = 100 and 300 m. A 1 to 2-yearly cycle is alsopresent between 200 and 300 m offshore. From x = 300 m onwards, the SPOshave periods between 1 . . x = 310 and between 2 . . . The questions we want to address in this section are whether there arecorrelations between the bathymetry and the NAO, the MWL and the MWHand whether these phenomena have a local influence or act throughout thewhole bathymetry. These phenomena are evaluated at a single location, asexplained in Sec. 2, while there are measurements of the bathymetry at sev-eral locations along each transect used in our analyses. The correlationsbetween the bathymetry and the NAO, the MWH and the MWL are com-puted using a simple correlation analysis. The correlations obtained at allbathymetric locations are shown in Fig. 5. The linear correlation analysisperformed at all transects seems to corroborate the local correlations be-tween embedded quasi-oscillations that were speculated about at transects58 and 190 when comparing the LPOs and SPOs (shown in Table 2) to theatmospheric/hydrodynamic patterns deduced with SSA (see Sect.4.3). How-ever, it is important to point out that a linear correlation analysis wouldpermit neither to characterise the bedlevel patterns nor to suggest nonlinearcorrelations as was achieved with SSA.14he strength of the correlation is given by the correlation coefficient, r ,shown in contour form in Figs 5. The maxima and minima of r show loca-tions where the correlations are strongest, independently of the magnitudeof r . In order to assess the overall strength of the (linear) correlation foreach of the variables the following criteria may be useful: a) strong corre-lation when 0 . < | r | (cid:54)
1; b) moderate correlation when 0 . < | r | (cid:54) . . < | r | (cid:54) . < | r | (cid:54) .
1. With these criteria, the percentage correlationbetween bedlevels and potential forcings may be computed for all the bathy-metric profiles. If all the bathymetric surveys are considered, the overallstrength of the linear correlations are as follows: 65.6% of the bathymetryis very weakly correlated and 34.4% is weakly correlated with NAO; 67.6%of the bathymetry is very weakly correlated and 32.4% is weakly correlatedwith MWH and; 28.3%of the bathymetry is very weakly correlated, 50.5%is weakly correlated and 21.2% is moderately correlated with MWL. Notethe values of these correlation coefficients are very low, but they only givea measure of the linear correlation between the bathymetric evolution andthe potential forcings, while SSA and MSSA give a nonlinear measure ofcorrelations.Table 3 shows for the first 4 EOFs from the MEOF analysis, the rowscorresponding to each of the potential forcings, for transects 58 and 190 andthen for the full set of bathymetric measurements at locations between andincluding these two transects. The potential forcings have been organisedso that the one at the top has the largest (absolute value) contribution forEOF1, the second for EOF2, and so on. It is clear that each potential forcingmay be unambiguously associated with a particular EOF. EOF4 has beenincluded to show that for this eigenvector all potential forcings generallyhave a contribution at least an order of magnitude smaller than for eigen-vectors EOF1 to EOF3. The eigenvector EOF1, corresponding to the largesteigenvalue, is correlated with NAO, regardless of whether we are consideringtransects 58 and 190, or all measured bathymetric profiles. Similarly, EOFs2 and 3 are correlated with MWH and MWL respectively in both cases.Also, the magnitudes of the components corresponding to potential forcingsin EOFs 1 to 3 are very close (in fact, they are equal for EOF1 and EOF2 atthe level of accuracy considered here) whether we consider the two transectsor all the bathymetric measurements together. Hence, the potential forcingsare correlated in the same way with transects 58 and 190 as they are withthe full set of measurements of the bathymetry and so we can use the latter15or subsequent analyses.The MEOF analysis indicates that the first 3 EOFs are strongly linkedto the three forcings considered and together they resolve 98 .
25% of thevariance. Note how large the resolved variance is compared to the values ofthe correlation coefficients. This shows that even though the correlations aresmall, correlations with any other phenomenon would be significantly smallerand thus we can confidently say that NAO, MWH and MWL are the onlyrelevant potential forcings. Also, since these three EOFs resolve so much ofthe variance, their associated PCs are used in the MSSA as input channelsinstead of the original dataset. This simplifies the MSSA considerably as thenumber of input channels is reduced from 1128 to 3.The first 14 MSSA components of the filtered system – that is, the systemconsisting of the first 3 PCs as input channels, were analysed. Each of thesecomponents has its corresponding eigenvalue, spatiotemporal PC (ST-PC)and spatiotemporal EOF (ST-EOF). The criteria presented by Plaut andVautard (1994) were used to identify coherent patterns within the ST-PCs.The method consists of finding pairs of consecutive eigenvalues which pass anumber of tests; the pairs passing such tests correspond to a quasi-oscillatorypattern. The technique is described in more detail in the Appendix. It wasfound that all 7 pairs considered passed the test. No more pairs were analysedas these pairs resolve most of the variance.Once pairs of eigenvalues corresponding to embedded quasi-oscillationsare identified, the periods of these quasi-oscillations can be assessed. Thefirst two pairs correspond to a yearly and a semi-yearly cycle (see Fig. 7;the time span is 120 months, or 10 years), and MSSA confirms that thesecycles are correlated with the monthly wave heights as their amplitude islargest at channel 2, which corresponds to MWH. The other pairs have largestamplitude at channel 1, hence they are all correlated with the NAO. The first4 of these pairs are linked to cycles of period of 3.92 mths, 3.14 mths, 8.62mths and 2.96 mths. The last pair is associated with two cycles, of 5.4 mthsand 2.98 mths. It is noteworthy that no collective interannual patterns seemto exist, even if such patterns have an important influence locally.16 . Discussion
At this transect a window length of 4 years chosen for the detrendinganalysis. However, the robustness of the results to small changes in windowlength was tested by varying the window between 3.3 and 6.6 years (so thewindow length span includes the sandbar life cycles observed north and southof the pier). It was found that the extrema had approximately the samevalues for all the window lengths considered, with some small deviations oftheir location along the transect depending on the window length considered.. The main difference observed occurred in the region of small variance, i.e.between 180 and 220 m offshore, where the trend signal did not capture anyof the extrema when M = 3 . . .
12 years. Again, the dominant patterns were found tobe robust to changes in window length, but the window length of 9 yearscaptured some of the fluctuations better.
Since the window length of 4 years was best at transect 58 detrendingwas performed initially at this window length. The trends shown in Fig. 1are at this window length for most of the transect except between 130 and170 m offshore. In this region the 4 year window did not resolve the travellingpatterns adequately, so it was necessary to experiment with slight changes towindow length to obtain a coherent picture. The amount of variance resolvedby the patterns, however, is the same as that with a window of 4 yrs, andthe characteristics of the travelling patterns are therefore robust (becausethese slight changes do not introduce more information). It seems that thesensitivity in window length in this region is due to the change in type ofbehaviour between the region closer to the shore and the patterns travellingseawards between 140 and 270 m (see the patterns developing from year 15onwards in Plot b,ii) of Fig. 1, for instance).
Our results compare well with the wavelet analyses presented by Reeveet al. (2007), which focused on temporal scales of variability at Tr. 62. We17hall use our results at Tr. 58 for the comparison, since Tr. 62 lies onlyaround 100 m south of Tr. 58 and hence sufficiently nearby for variations inunderlying oscillations to be small.Reeve et al. (2007) found that between 100 and 190 m offshore mostof the variance corresponded to periodicities between 8 and 12 . x = 260 m, 25 .
95% of the variance correspondedto periodicities between 1 . . x = 410 mthe pattern with a period of 6 − t = 0 to t = 3 yrs, while Lippmann et al.’s study took place between October 1986 andSeptember 1991, i.e. from t = 5 .
35 to 10 .
17 yrs. In both studies the shorelinemoved on and offshore, always lying between x = 100m and x = 140 m. Thisagrees well with the first spatial region identified with SSA at both windowlengths, which implies that these quasi-oscillatory patterns are linked to thoseof the shoreline.There is also good agreement between Birkemeier and Lippmann et al.’sstudies and this one in relation to the region identified as most dynamic. Inthis study it is found, for instance, that between 130 m and 350 m severalquasi-periodic cycles may dominate the dynamics (see Plots (a,ii) and (a,iii)as well as Plots (b,ii) and (b,iii) in Fig. 1); Birkemeier and Lippmann etal. also found that this was the most dynamic region in terms of numberof bars and in terms of bar motion. In their study, the inner bars movedon and offshore between 130 m and 330 m, and the outer bars generallybetween 200 and 450 m. Moreover, Lippmann et al. (1993) showed thatthe inner bar oscillates seasonally around a fixed position of around 200 m,while Birkemeier (1984, as cited by Lippmann et al. (1993)) showed that attransects 62 and 188 the inner bars had significant intra-annual variationsand remained generally shorewards of 200-210 m offshore. In this study it isalso found that several intra-annual patterns exist shorewards of 200-210 m,as shown in Plots (a,iii) and (b,iii) of Fig. 1. It is expected that a yearly cycle will be embedded within the wave dataand hence may potentially be correlated to the bathymetric yearly cycle. Toinvestigate this, we used the wave data from a wave gauge array, consistingof 15 pressure gauges mounted 0 . o and 108 o North, with 68 .
5% of them in the 73 o -90 o range. More-over, at this location the mean wave period is around 9 s, with a standarddeviation of 0 .
76 s; thus, there is little variation in the monthly averaged19ave period, T p . . . . . As shown in Sect. 3.3, the MEOF has permitted the isolation of the effectsof the different potential forcings on the bathymetry because each potentialforcing is most strongly correlated with one of the first 3 EOFs: the NAOwith EOF1, the monthly wave heights (MWH) with EOF2 and the monthlymean water level (MWL) with EOF3.Figs. 5 permit the identification of regions where each of the potentialforcings is more strongly correlated with seabed evolution, simply by locatingthe regions where the correlations are highly positive or highly negative. Notethat the region where the analysis was performed extends from 80 to 520 moffshore, rather than between 100 and 490 m, as was the case in the previoussections. As explained before, this was done to capture the behaviour aroundthe shoreline and the behaviour offshore more fully.Since the first EOF captures the largest proportion of the variance andit is linked to the NAO, it appears that the NAO is more strongly corre-lated with the seabed than are the waves at the time and space scales underconsideration. Plot 5(a) shows the distribution of correlations between theseabed and NAO, and highlights the complex interactions between them; theregions of positive and negative correlation are clearly visible in the plot. Itshows the correlations are either very weak (with | r | < .
1) or weak (with0 . < | r | (cid:54) . ± . . . Summary and Conclusions In summary, local and collective bathymetric quasi-periodic patterns ofoscillation were identified from monthly profile surveys spanning 26 years,from July 1981 to January 2006, at the USACE field research facility inDuck, North Carolina. The local analyses were based mostly on the outmosttransects of the surveyed domain, but the collective patterns were computedusing all surveyed bathymetric profiles, including those in the vicinity of thepier. Correlations with three potential forcings, namely the monthly waveheights (MWH), monthly mean water levels (MWL) - which include thestorm surges as well as the tide - and the large scale atmospheric index theNorth Atlantic Oscillation (NAO), were discussed.The local seabed patterns and their spatial correlations with NAO, MWHand MWL were first analysed using SSA, which showed the highly nonsta-tionary behaviour of the bedlevels along transects 58 and 190. It was shown1) that a quasi-yearly cycle was embedded at most bathymetric locationsand this is likely to be correlated with the quasi-yearly pattern within themonthly averaged wave climate; 2) the local behaviour is highly nonstation-ary as observed by previous studies; 3) patterns with periodicities in therange between 1 . .
6. Acknowledgements
This work was supported by an RCUK Academic Fellowship from the EP-SRC (Grant number EP/C508750/1) and a Research and Innovation Fellow-ship award from the University of Plymouth, UK. Thanks to Tom Lippmannand David Huntley for useful discussions. The editors and two anonymousreferees are thanked for valuable comments.
References
Akima, H., Oct 1970. A new method of interpolation and smooth curve fittingbased on local procedures. Journal of the ACM 17 (4), 589–602.Allen, M. R., Robertson, A. W., 1996. Distinguishing modulated oscillationsfrom coloured noise in multivariate datasets. Clim. Dyn. 12 (11), 775–784.Birkemeier, W. A., 1984. Time scales of nearshore profile changes. In: Pro-ceeding of the 19 th Conference on Coastal Engineering. New York (ASCE),pp. 1507–1521.Blackman, R. B., Tukey, J. W., 1958. The measurement of power spectrafrom the point of view of communication engineering. Dover, New York.Bowen, A. J., 1980. Simple models of nearshore sedimentation: beach profilesand longshore bars. In: McCann, S. B. (Ed.), The Coastline of Canada.Geological Survey of Canada, Ottawa, pp. 1–11.Bowen, A. J., Guza, R. T., 1978. Edge waves and surf beat. J. Geophys. Res.83, 1913–1920. 26owen, A. J., Inman, D. L., 1971. Edge waves and crescentic bars. J. Geo-phys. Res. 76.Broomhead, D. S., King, G. P., 1986. Extracting qualitative dynamics fromexperimental data. Physica D 20, 217–236.Cayocca, F., 2001. Long-term morphological modeling of a tidal inlet: theArcachon basin, France. Coast. Eng. 42, 115–142.Coco, G., O’Hare, T. J., Huntley, D. A., 1999. Beachcusps: A comparisondata and theoriesfor their formation. J. of Coastal Res. 15, 741–749.Dodd, N., Blondeaux, P., Calvete, D., de Swart, H. E., Falqu´es, A., Hulscher,S., R´o˙zy´nski, G., Vittori, G., 2003. Understanding coastal morphodynam-ics using stability methods. Jour. Coast. Res. 19 (4), 849–865.Evans, R. L., November 2004. Rising sea levels and moving shorelines.
Oceanus
WHOI Magazine.Falques, A., Coco, G., Huntley, D. A., 2000. A mechanism for the generationof wave-driven rhythmic patterns in the surf zone. Jour. Geophys. Res. -Oceans 105, 24071–24087.Gallagher, E. L., Elgar, S., Guza, R. T., 1998. Observations of sand barevolution on a natural beach. Jour. Geophys. Res. 103 (C2), 3203–3215.Garnier, R., Calvete, D., Falques, A., Caballeria, M., 2006. Generation andnonlinear evolution of shore-oblique/transverse sand bars. J. Fluid Mech.567, 327–360.Gehrels, W. R., Long, A. J., 2007. Sea level is not level: the case for a newapproach to predicting UK sea-level rise. Geography 93 (1), 11–16.Ghil, M., Allen, M. R., Dettinger, M. D., Ide, K., Kondrashov, D., Mann,M. E., Robertson, A. W., Saunders, A., Tian, Y., Varadi, F., Yiou, P.,2002. Advanced spectral methods for climatic time series. Rev. Geophys.40 (1), –1– –41.Ghil, M., Mo, K. C., 1991. Intraseasonal oscillations in the global atmosphere.part I: Northern Hemisphere and tropics. J. Atmos. Sci. 48, 752–779.27ilmore, R., Lefranc, M., 2002. The topology of Chaos; Alice in Stretch andSqueezeland, First Edition, John Wiley and Sons, Inc, New York, 518 pp.Grassberger, P., 1986. Do climate attractors exist? Nature 323, 609–612.Grassberger, P., Procaccia, I., 1983. Measuring the strangeness of strangeattractors. Phys. D 9, 189–208.Grinsted, A., Moore, J. C., Jevrejeva, S., 2004. Application of the crosswavelet transform and wavelet coherence to geophysical time series. Nonlin.Proc. Geophys. 11, 561–566Holman, R. A., Bowen, A. J., 1982. Bars, bumps, and holes: Models for thegeneration of complex beach topography. J. Geophys. Res. 87 (C1), 457–468.Holman, R. A., Lippmann, T. C., 1987. Remote sensing of nearshore barsystems – making morphology visible. Proceedings of Coastal Sediments’87 ASCE, 927–944Holman, R. A., Sallenger, A. H., 1986. High energy nearshore processes. EoSTrans. AGU 67, 1369–1371.Horrillo-Caraballo, J. M., Reeve, D. E., 2008. An investigation of the linkbetween beach morphology and wave climate at Duck, N.C., USA. Jour.Flood Risk Management 1 (2), 110–122.Howd, P. A., Birkemeier, W. A., 1987. Beach and nearshore survey data:1981-1984. CERC Field Research Facility. Tech. rep., U. S. Army Eng.Waterways Expt. Stn., Coast. Eng. Res. Ctr., Vicksburg, MS.Kirby, J. R., Kirby, R., 2008. Medium timescale stability of tidal mudflatsin Bridgwater Bay, Bristol Channel, UK: Influence of tides, waves andclimate. Cont. Shelf Res. 28, 2615–2629Klein, M. D., Schuttelaars, H. M., 2006. Morphodynamic evolution of double-barred beaches. J. Geophys. Res. 111 (C06017).Komar, P. D., 1971. Nearshore circulation and the formation of giant cusps.Geol. Soc. Am. Bull. 82, 2643–2650.28uriyama, Y., Ito, Y., Yanagishima, S., 2008. Medium-term variations of barproperties and their linkages with environmental factors at Hasaki, Japan.Mar. Geol. 248 (1-2), 1–126.Larson, M., Kraus, N. C., 1994. Temporal and spatial scales of beach profilechange, duck, north carolina. Mar. Geo. 117, 75–94.Lee, G., Nicholls, R. J., Birkemeier, W. A., 1998. Storm-driven variability ofthe beach-nearshore profile at duck, north carolina, usa, 1981-1991. Mar.Geol. 148, 163–177.Lee, G. H., Birkemeier, W. A., 1993. Beach and nearshore survey data:1985-1991. CERC Field Research Facility. Tech. rep., Defense TechnicalInformation Center OAI-PMH Repository [http://stinet.dtic.mil/oai/oai](United States).URL http://handle.dtic.mil/100.2/ADA266371 Li, Y., Lark, M., Reeve, D., 2005. Multi-scale variability of beach profiles atDuck: A wavelet analysis. Coast. Eng. 52 (12), 1133–1153.Lippmann, T. C., Holman, R. A., 1990. The spatial and temporal variabilityof sand bar morphology. Jour. Geophys. Res. 95 (C5), 11575–11590.Lippmann, T. C., Holman, R. A., Hathaway, K. K., 1993. Episodic, nonsta-tionary behavior of a double bar system at Duck, North Carolina, U. S. A.Jour. Coastal Res. 15, 49–75.Miller, J. K., Dean, R. G., 2006. Shoreline variability via empirical orthogonalfunction analysis: Part I. Temporal and spatial characteristics. Coast. Eng.54 (2), 111–131.Moron, V., Vautard, R., Ghil, M., 1998. Trends, interdecadal and interannualoscillations in global sea-surface temperatures. Climate Dynamics 14, 545–569.Mote, P. W., Clarke, H. L., Dunkerton, T. J., Harwood, R. S., Pumphrey,H. C., 2000. Intraseasonal variations of water vapor in the tropical uppertroposphere and tropopause region. Jour. Geophys. Res. 105 (D13), 17457–17470. 29icholls, R. J., Birkemeier, W. A., Lee, G., 1998. Evaluation of depth ofclosure using data from Duck, NC, USA. Mar. Geo. 148, 179–201.Niedoroda, A. W., Tanner, W. F., 1970. Preliminary study on transversebars. Mar. Geol. 9, 41–62.North, G. R., Bell, T. L., Cahalan, R. F., Moeng, F. J., 1982. Sampling errorsin the estimation of empirical orthogonal functions. Mon. Wea. Rev. 110,699–706.Papadopoulos, A., Tsimplis, M. N., 2006. Coherent coastal sea-level variabil-ity at interdecadal and interannual scales from tide gauges. J. Coast. Res.22 (3), 625–639.Pape, L., Ruessink, B. G., 2008. Multivariate analysis of nonlinearity insandbar behavior. Nonlinear Processes in Geophysics 15 (1), 145–158.URL
Plant, N. G., Holman, R. A., Freilich, M. H., Birkemeier, W. A., July 1999. Asimple model of interannual sandbar behaviour. Jour. Geophys. Res. 104,15755–15776.Plaut, G., Vautard, R., 1994. Spells of low-frequency oscillations and weatherregimes in the Northern Hemisphere. J. Atmos. Sci 51 (2), 210–236.Pruszak, Z., R´o˙zy´nski, G., Szmytkiewicz, M., Skaja, M., 1999. Quasi-seasonal morphological shore evolution response to variable wave climate.In: Coastal Sediments Conference Proceedings, ASCE. pp. 1081–1093.Ranasinghe, R., McLoughlin, R., Short, A., Symonds, G., 2004. The South-ern Oscillation Index, wave climate and beach rotation. Mar. Geol. 204,273–287.Rattan, S. S. P., Ruessink, B. G., Hsieh, W. W., 2005. Non-linear complexprincipal component analysis of nearshore bathymetry. Nonlin. ProcessesGeophys. 12 (5), 661–670.Reeve, D., Li, Y., Lark, M., Simmonds, D., May 2007. An investigation ofthe multi-scale temporal variability of beach profiles at Duck using waveletpacket transforms. Coast. Eng. 54 (5), 401–415.30eeve, D. E., Horrillo-Caraballo, J. M., Magar, V., 2008. Statistical analysisand forecasts of long-term sandbank evolution at Great Yarmouth, UK.Estuar. Coast. Shelf Sci. 79 (3), 387–399.Robertson, A. W., 1996. Interdecadal variability over the north pacific in amulti-century climate simulation. Clim. Dyn. 12, 227–241, in the appendixthere is a very clear explanation of the MSSA and the Monte Carlo teststo separate signal from noise and identify modulated oscillations from amultivariate timeseries.R´o˙zy´nski, G., Larson, M., Pruszak, Z., 2001. Forced and self-organised shore-line response for a beach in the southern Baltic Sea determined throughSingular Spectrum Analysis. Coast. Eng. 43, 41–58.R´o˙zy´nski, G., 2003. Data-driven modeling of multiple longshore bars andtheir interactions. Coast. Eng. 48, 151–170.R´o˙zy´nski, G., 2005. Long-term shoreline response to a nontidal, barred coast.Coast. Eng. 52, 79–91.R´o˙zy´nski, G., 2010. Long-term evolution of Baltic Sea wave climate near acoastal segment in Poland; its drivers and impacts. Ocean Eng. 37, 186–199.Ruelle, D., 1990. Deterministic chaos: The science and the fiction. Proc.Royal Soc. London A 427, 241–248.Ruessink, B. G., Kroon, A., 1994. The behaviour of a multiple bar system inthe nearshore zone of Terschelling: 1965–1993. Mar. Geol. 121, 187–197.Ruelle, D., 1990. Deterministic chaos: The science and the fiction. Proc.Royal Soc. London A 427, 241–248.Ruessink, B. G., Wijnberg, K. M., Holman, R. A., Kuriyama, Y., Van Enck-evort, I. M. J., 2003. Intersite comparison of interannual nearshore barbehavior. J. Geophys. Res. 108 (C8), 3249.Shand, R. D., Bailey, D. G., 1999. A review of net offshore bar migrationwith photographic illustrations from Wanganui, New Zealand. J. CoastalRes. 15, 365–378. 31hand, R. D., Bailey, D. G., Shepherd, M. J., 1999. An inter-site comparisonof net offshore bar migration characteristics and environmental conditions.J. Coastal Res. 15, 750–765.Shand, R. D., 2007. Bar splitting: system attributes and sediment budgetimplications for a net offshore migrating bar system. J. Coastal Res. SI50 (Proceedings of the 9th International Coastal Symposium, Gold Coast,Australia), 721–730.Southgate, H. N., M¨oller, I., 2000. Fractal properties of coastal profile evolu-tion at Duck, North Carolina. J. Geophys. Res. 105 (C5), 11489–11507.Southgate, H. N., Wijnberg, K., M., L., Capobianco, M., Jansen, H., 2003.Analysis of field data of coastal morphological evolution over yearly anddecadal timescales. part 2: Non-linear techniques. Jour. Coast. Res. 19,776–789.Vautard, R., Ghil, M., 1989. Singular spectrum analysis in nonlinear dynam-ics, with applications to paleoclimatic time series. Physica D 35, 395–424.Welch, P. D., June 1967. The use of fast fourier transform for the estimationof power spectra: a method based on time averaging over short, modifiedperiodograms. IEEE Trans. Audio Electroacoustics 15 (2), 70–73.Wijnberg, K. M., Kroon, A., 2002. Barred beaches. Geomorphology 48, 103–120.Wijnberg, K. M., Terwindt, J. H. J., 1995. Extracting decadal morphologicalbehaviour from high-resolution, long-term bathymetric surveys along theHolland coast using eigenfunction analysis. Mar. Geol. 126, 301–330.Winant, C. D., Inman, D. L., Nordstrom, C. E., 1975. Description of sea-sonal beach changes using empirical eigenfunctions. Jour. Geophys. Res.80, 1979–1986.Wright, L. D., Short, A. D., 1984. Morphodynamic variability of surf zonesand beaches: a synthesis. Mar. Geol. 56 (C7), 11575–11590.Yiou, P., Sornette, D., Ghil, M., 2000. Data-adaptive wavelets and multi-scale singular-spectrum analysis. Physica D 142, 254–290.32 . Description of methods
A.1. SSA and SDA
The fundamental theory of SSA has been discussed by Broomhead andKing (1986), Vautard and Ghil (1989) and Ghil et al. (2002) in relation to dis-crete time series analysis. The following summarises the SSA methodology.Suppose the bed level variable z ( t ), sampled at times t = iτ s , i = 1 , , . . . , N ( N being the length of the time series and τ s the sampling interval), at a givencross-shore and long-shore location ( x, y ), characterises the seafloor dynam-ical system. First, the noise part of the bed level time series may be identi-fied by computing the statistical dimension, S , which is found by construct-ing a secondary, M − dimensional sequence, Z n = ( z n , z n +1 , . . . , z n + M − ), for n = 1 , , . . . , N − M +1, where M < N is the window length and z n = z ( nτ s ).Time scales of the dynamics that can be reconstructed from this time seriesare between τ s and τ w = M τ s , where τ w is the window time span. Then, theeigenvalue problem Ce ( k ) = λ k e ( k ) , (1)must be solved, where C is the covariance matrix of the sequence Z n and λ k are the eigenvalues associated with the eigenvector e ( k ) . The square rootsof λ k are called the singular values, with their set being called the singularspectrum . The largest singular values relate to the directions that resolvemost of the variance, with the rest being a representation of the noise com-ponent and appearing as an approximately flat floor in a plot of singularvalue vs. singular value rank (with the singular values ordered from largestto smallest). The dimension S , then, will be given by the number of singularvalues above this floor.Now, each singular value λ k has an associated principal component (PC),defined as A k ( iτ s ) = M (cid:88) j =1 z ( { i + j − } τ s ) e ( k ) j . (2)Once the PCs have been identified, the part of the original time series relatedto that PC, or to a combination of K PCs, may be obtained by computing RC ( iτ s ) = 1 M i (cid:88) k ∈ K U i (cid:88) j = L i A k ( { i − j + 1 } τ s ) e ( k ) j , (3)where RC, the reconstructed component , is the part of the original signalrelated to the PC combination. The values of the renormalisation factor,33 i , and the lower and upper bound, L i and U i , respectively, depend on thetime at which the reconstruction is being evaluated.The frequencies that can be resolved with SSA depend on the choice ofwindow length, M , which is a free parameter of the method. However, thepatterns extracted should be robust to small changes in window length. Herethe depth time series were divided into a trend and a detrended signal. Thetrend generally contains a part found by linear regression of the data, aswell as the RC corresponding to the first few eigenvectors of the SSA at agiven window length M . These first few eigenvectors generally contain alloscillations within the original time series that have periods above M , whichin this case was taken as four years. From the trend thus computed, long-period oscillations (LPOs) were extracted, while from the detrended signalshort-period oscillations (SPOs) may be characterised. Cycles of three tofive years may fall in either signal because of the window length selected forthe detrending analysis (this will be discussed in more detail in subsequentsections).In natural systems, there is the additional problem of red noise (whichis noise composed of many random frequencies but, in contrast to whitenoise where all frequencies have equal intensity, the intensity of the lowerfrequencies is stronger) that can be removed from the data by using eitherMonte Carlo or chi-squared tests. In this work a chi-squared test was appliedon 100 red-noise surrogates; these surrogates were projected onto the SSAeigenvector basis, and the distribution of projections was approximated aschi-squared with 3 M/N degrees of freedom. Then those projections that felloutside a percentile range of 2 .
5% to 97 .
5% were identified as being part ofthe signal, and could thus be dissociated from the red noise in the time series.Spectral density analysis (SDA) is a well established technique to anal-yse oscillatory patterns in time series. The spectral density was computedusing Welch’s (Welch, 1967) periodogram method and the Blackman-Tukeycorrelogram approach (Blackman and Tukey, 1958). The former consists ofdividing the input signal into short segments and computing a modified pe-riodogram for each segment, thus leading to a set of periodograms which arefinally averaged to obtain an estimate of the spectral density function. In thelatter a windowed Fast Fourier Transform of the autocorrelation function isequated to the spectral density.It is well known that the spectral density is a measure of the energycontained in the system. By analysing the spectral density it is possible tocharacterise how the total energy is distributed and, in particular, to find the34requencies that contribute the most to it.
A.2. MEOF/Correlation analysis
Multivariate empirical orthogonal functions can be used to study therelationship between different physical variables that are sampled at the sameset of times. It can be thought of as a prefiltering step for MSSA/MEEOFanalysis: instead of performing the MSSA directly on the data, one firstfinds the MEOFs and uses time-lagged copies of the corresponding PCs inthe MSSA, in order to reduce the dimension of the covariance matrix.We construct a matrix B whose elements B ij = x ( j ) i are the values of the j th variable sampled at time iτ s , where i = 1 , , ..., N and j = 1 , , ..., L , where each vector x ( j ) contains the full time series fora particular variable. The full set of variables might include samples ofthe same physical variable at different locations and/or different physicalvariables at the same location(s).We next subtract from each element of x ( j ) the mean value of all itselements, to form a matrix F consisting of columns y ( j ) with zero mean: F ij = y ( j ) i = x ( j ) i − N N (cid:88) i =1 x ( j ) i . Since the columns of F represent different physical variables, we renor-malise them so that the results of the analysis are not influenced by the unitsin which we choose to measure them (Mote et al., 2000). Our data consists ofmeasurements of the bathymetry at 1125 separate locations over Duck beach(covering the area between y = −
91 to 1097 m along the shore and between x = 80 to 520 m offshore), together with measurements of the monthly meanwater levels (MWL) and monthly wave heights (MWH) at a wave gauge lo-cated near Duck beach, and a measurement of the North Atlantic Oscillation(NAO) spatially averaged over the whole surveyed area, the NAO being thepressure difference between Iceland and the Azores. Since all of these forc-ings were simultaneously available for only the last 19 years of bathymetricsurveys, we reduced the analysis to this duration. We considered two possi-ble choices for the renormalisation factor: the range of each variable or itsstandard deviation. We decided to use the range because we were interested35n understanding the effect of each possible forcing independently, and usingthe range separated these out in the calculated MEOFs to a greater degreethan using the standard deviation. Finally, since we have 1125 pieces of datafor the bathymetry at each time iτ s and only one piece for each possibleforcing, we divide the bathymetry measurements by 1125 so that the effectof the potential forcings is comparable with that of the spatial variations inthe bathymetry. From now on we consider F to have been renormalised inthis manner.Next, a singular value decomposition (SVD) of F permits the identifi-cation of bathymetric locations where the bathymetry may be strongly cor-related with each of the potential forcings. Such decomposition is of theform F = U SV T , where the columns of V are the EOFs of the SVD, or the spatial eigenvectors,the columns of U are the PCs, or the temporal eigenvectors and S is the di-agonal matrix of eigenvalues; the ratio of each eigenvalue divided by the sumof all the eigenvalues gives the fraction of the total variance explained by eacheigenvector. The EOFs represent the stationary patterns in the bathymetryand potential forcings that explain most of the observed variability, with thegreatest variance being resolved by the EOF corresponding to the largesteigenvalue and so on. The PCs provide the variation over time of the signand amplitude of the associated EOF.Each spatial eigenvector giving a direction in which each potential forc-ing is dominant gives an excellent proxy of the relative correlation of thatphenomenon with that particular spatial mode of the bathymetry. In fact,computing the correlations directly from the data and comparing the corre-lation distributions to the spatial EOF contours would show that correlationcontours and EOF contours are equivalent. This is because of the dominanceof the phenomena along well defined EOF directions. Since linear correlationsare conceptually more straightforward than EOF distributions, we presentand discuss the correlation plots instead of the EOF distributions. However,we use MEOF as a filter for the MSSA. A.3. MEEOF/MSSA
Multivariate Extended Empirical Orthogonal Functions (MEEOF) is thename given to Multichannel Singular Spectrum Analysis (MSSA) when sev-eral different physical variables are included in the analysis. Such a technique36as been useful, in particular, in the analysis of the behaviour of climatic phe-nomena, with the aim of extracting modulated oscillations from the colourednoise commonly present in natural systems (Vautard and Ghil, 1989; Moteet al., 2000). Its use has focused on accurate identification and characteri-sation of the system’s temporal and spectral properties. MSSA is fully de-scribed in Allen and Robertson (1996) and references therein. Here a shortdescription is presented for the reader interested in the methodology and itscapabilities.Suppose we have the dataset, (cid:110) y ( j ) i ; i = 1 , . . . , N, j = 1 , . . . , L (cid:111) , discussedin the previous section. Several techniques to produce a matrix of time-laggedvectors may be proposed. Here the method used by Allen and Robertson(1996) is adopted. It consists of creating the matrices ˜ Y ( j ) , formed by set-ting the channel to j , 1 ≤ j ≤ L , sliding a vector of length M down from i = 1 to i = N (cid:48) ≡ N − M + 1 and ordering the resulting vectors from left toright as the columns of ˜ Y ( j ) . Thus each ˜ Y ( j ) has size M × N (cid:48) . This leads toa trajectory matrix , ˜ Y = ˜ Y (1) ...˜ Y ( L ) (4)of size M L × N (cid:48) .A singular value decomposition of ˜ Y may then be performed, such that˜ Y = P ˜ Y Λ / Y E T ˜ Y , (5)where P ˜ Y consists of the temporal principal components (T-PCs). As equa-tion 5 implies, these are the left-singular vectors of ˜ Y . The matrix E ˜ Y consists of the spatio-temporal empirical orthogonal functions (ST-EOFs)corresponding to the right-singular vectors of ˜ Y . Finally, Λ ˜ Y is the diago-nal matrix of variances associated with these orthogonal bases (Robertson,1996).The bases P ˜ Y and E ˜ Y are the eigenvectors of the covariance matrices C P ˜ Y = 1 M L ˜ Y ˜ Y T and C E ˜ Y = 1 N (cid:48) ˜ Y T ˜ Y , (6)respectively. However, it is only for the smallest of C P ˜ Y and C E ˜ Y that either P ˜ Y or E ˜ Y is full-rank (Allen and Robertson, 1996; Robertson, 1996). There-fore, the choice of M determines which of the eigenbases is to be analysed.37owever, the method of construction of ˜ Y implies that the column vectorsof ˜ Y T are time-lagged copies at a given channel, so it is convenient to choose N (cid:48) small enough to satisfy N (cid:48) < M L , so that the full rank matrix is C E ˜ Y .As mentioned in Sect. A.2, rather than performing the MSSA over thedata and having to resolve a very large covariance problem, it is better toperform a data reduction procedure by prefiltering the original data with theMEOF and apply the MSSA to the resulting EOFs. In this way, the mostimportant information is compressed into a small number of variables (Plautand Vautard, 1994) and the analysis of the coherent patterns is simplifiedsignificantly. In this case the first three eigenvalues of the MEOF resolvemost of the variance and also are the ones in the directions of the threepotential forcings considered, so the MSSA can be applied to these threeeigenfunctions rather than to the data.The first 14 MSSA components of the filtered system – that is, the systemconsisting of the first 3 MEOF PCs as input channels, were analysed. Eachof these components has its corresponding eigenvalue, spatiotemporal PC(ST-PC) and spatiotemporal EOF (ST-EOF). A window length of M = 120,equivalent to N (cid:48) = 9 years, was chosen. The criteria presented by Plaut andVautard (1994) were used to identify coherent patterns within the ST-PCs.The method consists of finding pairs of consecutive eigenvalues ( k, k + 1) forwhich: • the two eigenvalues are nearly equal, that is they are both within themargins of error identified by North’s rule of thumb (North et al., 1982), • the two corresponding time sequences described by the ST-PCs arenearly periodic, with the same period and in quadrature, • the associated ST-EOFs are in quadrature and • the lag correlation between ST-EOFs of order k and k + 1 itself showsoscillatory behaviour (Ghil and Mo, 1991).Pairs that pass all these tests are associated with quasi-oscillations withperiods equal to the period of the ST-PCs - which is also the period of theST-EOFs - in the embedding space, here spanning 121 months.The first test was to inspect pairs of consecutive eigenvalues λ ( k ) , λ ( k +1)to see if they satisfied North’s rule of thumb. It was found that eigenvalue λ ( k + 1) lies within the interval [ λ ( k )(1 − (cid:112) /N ); λ ( k )(1 + (cid:112) /N )], with38 = 228 months being the sampling duration, for all pairs considered, sothe 7 of them satisfy the first criterion. Then plots of the ST-PCs for the 14MSSA components at channel 1 showed that all pairs considered are quasi-periodic with same period and in phase quadrature so they also satisfy thesecond criterion. This led to a check of the quadrature criterion for the ST-EOFs associated with each ST-PC pair. Since the ST-EOFs are the samefor all the PCA channels in the prefiltering, only the results for channel 1 arenecessary. The ST-EOF plots are shown in Fig. 7 for the first 6 pairs, withthe solid line corresponding to the eigenvector with smaller index. All thepairs (including pair 13-14, not shown) pass the test. The final test consistsof analysing the lagged correlation plots for the pairs which have satisfied allprevious tests. The plots, shown in Fig. 8, should depict oscillatory functions.Again, all pairs pass this test. 39 igure Captions: Caption 1: Data contours for transects 58 (top, North of the pier) and 190(bottom, South of the pier) with: beach elevation (in metres) displayed inFigs. (a,i) and (b,i); trend signal elevation (in metres) displayed in Figs. (a,ii)and (b,ii) and; detrended signal elevation (in meters) displayed in Figs. (a,iii)and (b,iii). The upper panels in Figs. (a,ii), (a,iii), (b,ii) and (b,iii) show thepercentage of variance resolved by the corresponding signal.Caption 2: Quasi-oscillations identified within the trend (left) and the de-trended (right) signals together with the resolved variance, at window lengthsof 9 years, for transects 58 (top) and 190 (bottom).Caption 3: Spectral density plot for the reconstructed component of the SPOsat x = 160 m for Transect 58. The thin lines indicate the 95% confidenceinterval for the spectral estimates.Caption 4: SSA reconstruction of wave dynamics, with the first fundamentalpair shown at the top and the second fundamental pair shown at the bottom.Caption 5: Correlation Contours between the seabed and the NAO, theMWH and the MWL, respectively. Tr. 58 (top, Northern side) and Tr. 190(bottom, Southern side), are highlighted (solid lines). The pier is located at y = 513 m (dashed line).Caption 6: First six pairs of ST-EOFs of the MSSA components linked to po-tentially quasi-oscillatory pairs at all channels, with their associated resolvedvariance. The abscissa represents time (in months) and spans the windowlength M = 121 months.Caption 7: ST-PCs for consecutive MSSA components at channel 1; theST-PCs span the window length N (cid:48) = 108 months.Caption 8: Lagged correlations for all ST-PC pairs.40 igures: x (m) t i m e ( y ea r s ) (a,i)
200 40005101520 −2−1012 050100 V a r x (m) t i m e ( y ea r s ) (a,ii)
200 40005101520 −1−0.500.51 050100 V a r x (m) t i m e ( y ea r s ) (a,iii)
200 40005101520 −101x (m) t i m e ( y ea r s ) (b,i)
200 40005101520 −2−1012 050100 V a r x (m) t i m e ( y ea r s ) (b,ii)
200 40005101520 −1−0.500.51 050100 V a r x (m) t i m e ( y ea r s ) (b,iii)
200 40005101520 −101
Figure 1: Data contours for transects 58 (top, North of the pier) and 190 (bottom, Southof the pier) with: beach elevation (in metres) displayed in Figs. (a,i) and (b,i); trend signalelevation (in metres) displayed in Figs. (a,ii) and (b,ii) and; detrended signal elevation (inmeters) displayed in Figs. (a,iii) and (b,iii). The upper panels in Figs. (a,ii), (a,iii), (b,ii)and (b,iii) show the percentage of variance resolved by the corresponding signal. V a r x (m) x (m) t i m e ( y ea r s )
100 200 300 40005101520 −1−0.500.5 (a) LPOs, Tr. 58 V a r x (m) x (m) t i m e ( y ea r s )
100 200 300 40005101520 −1−0.500.5 (b) SPOs, Tr. 58 V a r x (m) x (m) t i m e ( y ea r s )
100 200 300 40005101520 −1−0.8−0.6−0.4−0.200.2 (c) LPOs, Tr. 190 V a r x (m) x (m) t i m e ( y ea r s )
100 200 300 40005101520 −1−0.500.51 (d) SPOs, Tr. 190Figure 2: Quasi-oscillations identified within the trend (left) and the detrended (right)signals together with the resolved variance, at window lengths of 9 years, for transects 58(top) and 190 (bottom). l og ( S D ) Figure 3: Spectral density plot for the reconstructed component of the SPOs at x = 160 mfor Transect 58. The thin lines indicate the 95% confidence interval for the spectralestimates. .2 8.3 12.5 16.7 −0.500.5 H e i gh t ( m ) −0.500.5 Time ( years ) H e i gh t ( m ) Figure 4: SSA reconstruction of wave dynamics, with the first fundamental pair shown atthe top and the second fundamental pair shown at the bottom.
100 200 300 400 50002004006008001000 x (m) y ( m ) s −0.2−0.15−0.1−0.0500.050.1 (a) NAO Correlations
100 200 300 400 50002004006008001000 x (m) y ( m ) s −0.15−0.1−0.0500.050.1 (b) MWH Correlations
100 200 300 400 50002004006008001000 x (m) y ( m ) s −0.3−0.2−0.100.10.20.3 (c) MWL CorrelationsFigure 5: Correlation Contours between the seabed and the NAO, the MWH and theMWL, respectively. Tr. 58 (top, Northern side) and Tr. 190 (bottom, Southern side), arehighlighted (solid lines). The pier is located at y = 513 m (dashed line). Time (years) S T − E O F Time (years) S T − E O F Time (years) S T − E O F (a) Pair 1-2 Time (years) S T − E O F Time (years) S T − E O F Time (years) S T − E O F (b) Pair 3-4 Time (years) S T − E O F Time (years) S T − E O F Time (years) S T − E O F (c) Pair 5-6 Time (years) S T − E O F Time (years) S T − E O F Time (years) S T − E O F (d) Pair 7-8 Time (years) S T − E O F Time (years) S T − E O F Time (years) S T − E O F (e) Pair 9-10 Time (years) S T − E O F Time (years) S T − E O F Time (years) S T − E O F (f) Pair 11-12Figure 6: First six pairs of ST-EOFs of the MSSA components linked to potentially quasi-oscillatory pairs at all channels, with their associated resolved variance. The abscissarepresents time (in months) and spans the window length M = 121 months. Time (years) S T − P C Time (years) S T − P C Time (years) S T − P C Time (years) S T − P C Time (years) S T − P C Time (years) S T − P C Time (years) S T − P C Figure 7: ST-PCs for consecutive MSSA components at channel 1; the ST-PCs span thewindow length N (cid:48) = 108 months. Lag Corr, Pair 1 and 2
Time lag (years) −2 −1 0 1 2−101
Lag Corr, Pair 3 and 4
Time lag (years) −2 −1 0 1 2−101
Lag Corr, Pair 5 and 6
Time lag (years) −2 −1 0 1 2−101
Lag Corr, Pair 7 and 8
Time lag (years) −2 −1 0 1 2−101
Lag Corr, Pair 9 and 10
Time lag (years) −2 −1 0 1 2−101
Lag Corr, Pair 11 and 12
Time lag (years) −2 −1 0 1 2−101
Lag Corr, Pair 13 and 14
Time lag (years)
Figure 8: Lagged correlations for all ST-PC pairs. ables: TN 190 188 62 58 y -91 0 1008 1097∆ y
32 17 98 91
Table 1: Correspondence between transect number and transect longshore position for thefour best surveyed transects at Duck, North Carolina, USA.
T N is the transect number, y the averaged longshore position ∆ y the absolute error around this average. ransect 58 Transect 190x (m) LPOs (yrs) SPOs (yrs) LPOs (yrs) SPOs (yrs)120 6 . .
7) 15 . .
5) 1(0 . . . . .
1) 8(1) 1(0 . . . . .
1) 1(0 .
2) 1 . . . .
5) 1(0 . . . . . . .
6) 1(0 .
1) 6 . .
5) 2 . . . . .
3) 1(0 .
1) 13(1), 1 . .
1) 2 . . . . . .
8) 1(0 .
1) 1 . .
3) 2 . . . . .
6) all below 1 5 . .
2) 2 . . . . . . .
4) 0 . . . .
1) 3 . .
4) 4 . . . . . . .
2) 2 . .
5) 10(1) 1 . . . . .
8) 1(0 . . .
2) 12(1) 2 . . . .
9) 1(0 . . .
1) 12(1) 2 . . . . . . .
7) 1 . . . .
6) 14(2) 2 . . . . . . . .
7) 1 . . . .
3) 5 . .
2) 1 . . . . . . . .
8) 1 . .
3) 6 . .
5) 3 . . . . . . . .
9) 6 . . .
4) 1(0 . . . . . . .
5) 1 . . . .
1) 6 . . . .
6) 1(0 . . . . .
2) 2 . .
6) 5 . .
4) 1(0 . . . . . . .
9) 1 . .
1) 4 . .
2) 1(0 . . . . . .
8) 2 . .
9) 3 . . . .
2) same as 410480 10 . .
2) 2 . . . .
8) 3 . .
2) same as 410490 7 . . .
5) 4 . .
2) 1(0 . . . . . Table 2: Periods of dominant oscillatory patterns at selected locations along transects 58and 190, with the uncertainty or confidence interval in brackets. OF values with all seabed surveysForcing: EOF1 EOF2 EOF3 EOF4NAO 0 .
998 -0 .
067 -0 .
012 -0 . .
067 0 .
998 0 . . .
012 0 .
001 -0 .
999 0 . .
998 -0 .
067 -0 .
013 -0 . .
067 0 .
998 0 . . .
012 0 .
001 -0 .