PUBLICATIONS Journal of Geophysical Research: Oceans RESEARCH ARTICLE Wave climatology in the Apostle Islands, Lake Superior 10.1002/2014JC010278 JoshuaD.Anderson1,ChinH.Wu1,andDavidJ.Schwab2 KeyPoints: (cid:2)WaveclimateoftheApostleIslands 1DepartmentofCivilandEnvironmentalEngineering,UniversityofWisconsin-Madison,Madison,Wisconsin,USA,2U-M inLakeSuperiorfor35yearwas WaterCenter,UniversityofMichigan,AnnArbor,Michigan,USA hindcast (cid:2)Statisticsofthewaveclimatereveal thespatialvariabilityofwave Abstract ThewaveclimateoftheApostleIslandsinLakeSuperiorfor35year(1979–2013)washindcast properties andexaminedusingathird-generationspectralwavemodel.WavemeasurementswithintheApostle (cid:2)AnincreasingtrendofSWHisfound duetoclimatechange IslandsandoffshoreNOAAbuoyswereusedtovalidatethemodel.Statisticsofthesignificantwaveheight, peakwaveperiod,andmeanwavedirectionwerecomputedtorevealthespatialvariabilityofwaveproper- Correspondenceto: tieswithinthearchipelagoforaverageandextremeevents.Extremevalueanalysiswasperformedtoesti- C.H.Wu, matethesignificantwaveheightatthe1,10,and100yearreturnperiods.Significantwaveheightsinthe [email protected] interiorareasoftheislandsvaryspatiallybutareapproximatelyhalfthoseimmediatelyoffshoreofthe islands.Duetoreducedwintericecoverandaclockwiseshiftinwinddirectionoverthehindcastperiod, Citation: long-termtrendanalysisindicatesanincreasingtrendofsignificantwaveheightsstatisticsbyasmuchas Anderson,J.D.,C.H.Wu,and D.J.Schwab(2015),Waveclimatology 2%peryear,whichisapproximatelyanorderofmagnitudegreaterthansimilaranalysisperformedinthe intheApostleIslands,LakeSuperior,J. globaloceanforareasunaffectedbyice.Twoscientificquestionsrelatedtowaveclimateareaddressed. Geophys.Res.Oceans,120,4869–4890, First,thewaveclimatechangeduetotherelativeroleofchangingwindfieldsoricecoversoverthepast35 doi:10.1002/2014JC010278. yearswasrevealed.Second,potentialblufferosionaffectedbythechangeofwaveclimateandthetrendof lowerwaterlevelsintheApostleIslands,LakeSuperiorwasexamined. Received1JUL2014 Accepted10JUN2015 Acceptedarticleonline15JUN2015 Publishedonline14JUL2015 1.Introduction Informationaboutthewaveclimateiscrucialtomanyenvironmentalandsocietalissues.Forexample,inthe ApostleIsland,areaofLakeSuperiorrecreationalboatinginvolvestensofthousandsofsailboats,motorboats, andkayaksperyear[Kraftetal.,2007].Extremewaveconditionscanposenavigationhazardsforrecreational watercraftandresultintragicdrowningincidents(DuluthNewsTribune,12September2010and9September 2011). Safe design of coastal and offshore structures also requires accurate knowledge of wave climate extremes[Panchangetal.,2013].Furthermore,waveclimatecanaffectmanygeomorphicandecologicproc- esses.Particularly,coastalwaveenergyhasbeendirectlyrelatedtoshorelinedamage,propertyloss,andbluff recession[Meadowsetal.,1997;Swensonetal.,2006;LinandWu,2014].Lakebedsedimentresuspensionhas beenshowntobedominatedbywave-generatedshearstressduringlargestormevents[Schwabetal.,2006]. Approximately3%ofthecoastlineintheApostleIslandsNationalLakeshoreofLakeSuperiorisdesignatedas wetland[Kraftetal.,2007],whichissignificantlyaffectedbywaveexposure[Thomasenetal.,2013].Lastbut notleast,laketroutareeconomicallyimportantspeciesthatpreferentiallyspawninathresholdofelevated wavedisturbanceorturbulence[FitzsimonsandMarsden,2014].Spatialinformationaboutthewaveclimate canprovideinsightintothehistoricalandpresentsuccessofembryonicsurvivaloflaketroutintheApostle Islands[CoberlyandHorrall,1980;Schrametal.,1995].Inviewoftheseimportantconsequences,itisimpera- tivetocharacterizewaveclimatologyaroundtheApostleIslandsinLakeSuperior. Determination of wave climate statistics often relies on analysis of historical wave records. Many of the aforementionedwave-dependentissuesareassessedusingthestatisticaldistributionofthreewavecharac- teristics:SignificantWaveHeight(SWH),MeanWaveDirection(MWD),andPeakWavePeriod(PWP),i.e.,the periodatthepeakoftheenergyspectrum.Toaccountforinterannualvariability,theWorldMeteorological Organization’s(WMO)standarddurationforanevaluationofclimateis30years[WorldMeteorologicalOrga- nization,2007].HistoricalmeasurementsbytheNationalDataBuoyCenterofNOAAarethebestcandidate to meet this criterion. In the Great Lakes, these buoys have been in operation since the early 1980s, i.e., approximately30years.However,thespatialcoverageofwavebuoysislimitedtoafewlocationsperlake. Thebuoylocationsaretypicallyoffshore,farawayfromcoastalareaswherewavecharacteristicsaresignifi- VC2015.AmericanGeophysicalUnion. AllRightsReserved. cantly altered by nearshore transformation processes. Consequently, wave hindcast through modeling is ANDERSONETAL. WAVECLIMATOLOGYINLAKESUPERIOR 4869 Journal of Geophysical Research: Oceans 10.1002/2014JC010278 usuallyemployedtofillinthepaucityofmeasuredwavedataandprovidethespatialvariabilityofwavecli- mate[CairesandSterl,2005;Jensenetal.,2012]. Inthepast,modelsforhindcastingwavedatasetshavebeensuccessfullyappliedtoobtainwaveclimate statistics such as annual means, seasonal variability, return level estimates, correlation with atmospheric cycles,andlong-termtrendsonbothglobal[SterlandCaires,2005;Chawlaetal.,2013;Stopaetal.,2013a] andregionalscales[Panchangetal.,2008;Watersetal.,2009;Dodetetal.,2010].IntheGreatLakes,hind- castingofthewaveclimatehasbeenledbytheU.S.ArmyCorpsofEngineersundertheWaveInformation Studies(WIS)project[ResioandVincent,1978;Hubertzetal.,1991;Jensenetal.,2012].Recently,theWISpro- jecthascompletedawavemodelhindcastforLakeSuperiorfortheperiodfrom1979to2012.Datafrom the WIS project include various products like historical time series of SWH, monthly means and maxima, extreme value analysis of SWH, and wave percent occurrence tables (http://wis.usace.army.mil/hindcasts. shtml).Nevertheless,dataarelimitedtotheselectedWISlocationsintheLakeSuperioranddonotinclude wavesinsidemanyshelteredareas,e.g.,theApostleIslands.Whileotheroperationalwavemodels,including the Donelan model developed by the Great Lakes Environmental Research Laboratory in the early 1980s [Liu et al., 1984] and Wave Watch III model developed by National Centers for Environmental Prediction withimprovedmodelphysics[Alvesetal.,2014],havebeenusedforLakeSuperior,theoverallgridresolu- tions range from 3 to 10 km. With the dimensions of the Apostle Islands as small as 250 m and adjacent islandsseparatedbyaslittleas1km,nowavehindcastdatainsidetheislandshavebeenreported. Giventherecentfindingofclimatechange,assessingthefutureofthewaveclimateisincreasinglyimportant. Intheocean,waveclimatetrendshavebeenquantifiedduetochangesinwindspeedanddirection[Caires andSterl,2005;Dodetetal.,2010].InLakeSuperior,increasingsummerwindspeedsoverLakeSuperiorhave been shown due to a decreased water-air temperature gradient [Desai et al., 2009], supporting the finding that water temperatures are rising faster than the air [Austin and Colman, 2007]. Liu and Ross [1980] showedthatunstableconditionscanenhancewavegrowthwhenwatertemperaturesarewarmerthanthe air.GiventheobservedriseinwatersurfacetemperaturesofLakeSuperiorandalengtheningofthestratified season[AustinandColman,2007],atmosphericunstabilitycanfurtherelevatethewaveclimate.Furthermore, trendsofthewaveclimateintheGreatLakescanbeaffectedbytheobservedlossoficecover.Schwabetal. [2006] showed that decreased ice cover in the southern basin of Lake Michigan allowed for greater occur- rence of large sediment resuspension events forced by waves. In the Apostle Islands, Howk [2009] showed thaticecoverhasbeendecreasingatahigherratethantheaveragerateofLakeSuperior,suggestingthat thewaveclimateintheApostlesisparticularlysusceptibletoclimatechange.Whilemuchevidenceexiststo indicatethewaveclimateintheApostleIslandsisincreasing,therehasbeennoquantificationofthesetrends intheApostlesoranywhereontheGreatLakes.Additionally,therehasbeenlittleexaminationofhowwave climatetrendsaffectedbyicelosscomparedtothosefoundinareasunaffectedbyice.Lastly,relativerolesof windsandicecoversinmodifyingthewaveclimateintheApostleIslandsareunclear. Bluffrecession,stronglycorrelatedwithwaveclimate[Swensonetal.,2006],isasignificantissueaffecting theeconomicdevelopmentandpropertysafelyinthecoastalcommunityintheGreatLakes[Angel,1995; Meadowsetal.,1997].ConservationoftheshorelineiscruciallyimportantintheApostleIslands,possessing thelargestcollectionoflighthousesintheNationalParkSystem[Kraftetal.,2007].Inthepast,characteriza- tion of coastal bluff stability and recession within the Apostles has had to be performed with the limited informationofthewaveclimate[Swensonetal.,2006;Pendletonetal.,2007].Asaresult,difficultiesinpro- jecting bluff recession in Lake Superior have been recognized and acknowledged. To remedy this issue, obtaining reliable wave climate are desperately needed since waves typically play a dominant role in the recessionratesofcoastalbluffs.Forexample,wavescanelevatebluffrecessionbyerodingnearbylakebed sediments,aprocessknownasdowncutting[LinandWu,2014].Additionally,wavescanerodesedimentsat thetoeof a bluff,resultingin asteeperslopethat becomesunstableandis morepronetofailure[Brown etal.,2005].Meanwhile,waterlevelscanalsoaffectblufferosionbyraisingorloweringtheimpactheightof wavesonablufftoe[Meadowsetal.,1997;Castedoetal.,2013].Brownetal.[2005]andSwensonetal.[2006] showed that bluff recession is correlated well with the magnitude-based wave impact height, defined as theelevationofwaverunupminustheelevationofablufftoe.WhilethemeanwaterlevelsonLakeSupe- rior have been steadily decreasing over the past 30 years [Assel et al., 2004; Gronewold et al., 2013], it is unclearifthelowerlakelevelsreducetheelevatedwaveclimateonoverallimpactofcoastalbluffsinApos- tleIslands. ANDERSONETAL. WAVECLIMATOLOGYINLAKESUPERIOR 4870 Journal of Geophysical Research: Oceans 10.1002/2014JC010278 The objective of this paper is to analyze the wave climate of the Apostle Islands in Lake Superior to the addressthetwosciencequestions.First,howwouldthewaveclimatechangeduetochangingwindfields or ice covers over the past 35 years? Second, how would the combined changing wave climate and the trend of lower water levels affect the coastal bluff erosion? To answer these questions, a wave model of LakeSuperiorwithjustifiablespatialresolutionintheApostleIslandsregionisdevelopedandappliedina 35yearhindcast(1979–2013).WavemeasurementscollectedwithintheApostleIslandsandfromoffshore NOAAbuoysareusedtovalidatethemodel.Statisticsofthewaveclimateincludingsignificantwaveheight, peakwaveperiod,andmeanwavedirectionareobtainedforaverageandextremeevents.Extremevalue analysisisperformedtoestimatethesignificantwaveheightforarangeofreturnperiods.Long-termtrends ofwaveheightarequantifiedtoevaluatetheeffectsofclimatechangeonwaveclimatestatistics.Finally, theeffectsofwaveclimateandlowwaterleveltrendonblufferosionareexamined. The paper is organized as follows. Description of the bathymetry and atmospheric climate in the Apostle Islandsisprovidedinsection2.Section3describesthedatasourcesformodelinputsofwindandice,wave modelperformance,thestatisticalmethodsappliedintheanalysisofthehindcastwaves,andthemethod forbluffrecessionanalysis.Section4presentstheresultsofthehindcastincluding:statisticsofthreewave climate parameters (SWH, MWD, and PWP), return value estimates of SWH, long-term trends in SWH, and trendsinwaveimpactsoncoastalbluffs.Section5discussestheperformanceandlimitationofwavemod- eling, effects of ice-time on wave statistics, the role of ice and wind on long-term wave climate statistics, and the role of water level and wave climate on bluff recession. Finally, a summary of the major findings fromthisstudyisgiveninsection6. 2.StudySite TheApostleIslandsarelocatedofftheBayfieldPeninsula(BP)inthewesternarmofLakeSuperior.Including all21islandsandtheadjacentmainland,thesitecontainsover300kmofGreatLakescoastline[Kraftetal., 2007]. Water depths within the archipelago are typically less than 60 m except for a deep channel in the east,whichrangesfrom100to140m(seeFigure1).Alongthenorthwesternboundary,watersareshallow withdepthsrangingfrom40mtolessthan20m.Multipleshoalsandheadlandsexistpredominantlyatthe outeredgesoftheislandnetworkwithnortheastalignment,mostnotablytheGullShoal(GS)attheeastern exterior.ANOAAmeterologicalstation(DISW3),locatedonthenorthernboundaryofthearchipelago,has collectedhourlywindmeasurementssince1983.ThehistoricalrecordofDISW3indicatesthattheaverage windspeedis5.03m/swithameandirectiondirectlyfromthewest.Thestrongeststormwindsblowfrom thenortheastbutlargewesterlywindeventsareequallycommon.Icecoverhasahighinterannualvariabili- tyandmaybepresentfromDecembertoMay[Assel,2003;Kraftetal.,2007].IcecoverattheCityofBay- field,locatedontheeasterncoastoftheBP,hasbeendecreasingatarateofapproximately14.7daysper decadesince1975[Howk,2009]. 3.Methods 3.1.Data 3.1.1.WindField In this study, gridded 10 m winds (U ) from the Climate Forecast System Reanalysis (CFSR) were usedto 10 drivethewavemodel.TheCFSRproductsareproducedbytheNationalCenterforEnvironmentalPrediction (NCEP)usingacoupledglobalatmospheric-oceanic-ice-landmodelwithadvanceddataassimilationtechni- quesandanextensivedatabaseofmeteorologicalobservations[Sahaetal.,2010].TheoriginalCFSRdata setspansfrom1979to2010andcontinuestooperateasthesecondversionoftheClimateForecastSystem (CFSv2)withmultipleimprovementsoverCFSR,includingahigherspatialresolution[Sahaetal.,2014].Spa- tialresolutionoftheCFSRisapproximately38km((cid:3)20kmforCFSv2)andtemporalresolutionis1h,asig- nificant upgrade from previous global reanalysis products. The spatial resolution rivals focused regional products like the North American Regional Reanalysis. Data can be accessed through the National Center forAtmosphericResearch(http://rda.ucar.edu/pub/cfsr.html). ComparisonofCFSRwindswithmeasurementsfrommeteorologicalstationsatthestudysitesandoffshore buoyswasmadehere.Atmosphericstability,whichcanaffectwavegrowth,wasaccountedforbyadjusting CFSR gridded winds to an equivalent neutrally stable wind speed [Liu and Schwab, 1987]. Specifically, all ANDERSONETAL. WAVECLIMATOLOGYINLAKESUPERIOR 4871 Journal of Geophysical Research: Oceans 10.1002/2014JC010278 Figure1.ApostleIslandsbathymetrymapandlocationonLakeSuperior.Sitesidentifiedwithacircle(•)aremeteorologicalorwaveobser- vationstations.Sitesidentifiedwithatriangle(~)areUSACEWISstations.Sitesidentifiedwithasquare((cid:2))arelocationsanalyzedfor waveimpactsonbluffs.TheunstructuredSWANmeshappliedinthehindcastisoverlaidonthemap.Thenonhydrostaticphase-resolving modeldomainisidentifiedasaboxfilledwithdiagonallinesatapproximately478Nand290.98E. measured data wereconvertedto 10 m height using the Monin-Obukhov similarity theory.Figure 2 shows statisticalcomparisonsofCFSRandmeasuredwindsattwocoastalstations(DISW3andSXHW3)andtwooff- shorebuoylocations(St.45006andSt.45001).Ingeneral,CFSRwindsmatchwellwiththeobservations,espe- ciallyatlarge windspeeds,but CFSRwinds slightly overpredictmeasuredvalues(positivebias) with higher biasesattheoffshorebuoys.Ofthefourstationsshown,onlySXHW3waslocatedinalandcellintheCFSR model.Consequently,CFSRunderpredictsobservedwindspeedspresumablyduetoahigherdragcoefficient implementedforland,whichisconsistentwithfindingsoftheCFSRwindsneartheland-seainterfaceinthe ocean[Chawlaetal.,2013].Toaccountforthisunderprediction,weextrapolatedCFSRwindmagnitudesfrom sea/watergridstolandgridsandretainedtheoriginaldirectionofthelandgrid.Abetteragreementbetween theextrapolatedandmeasuredwindsatthehigherquantilesisshowninFigure2c(plusmarkers). Water temperatures were obtained from the ‘‘Great Lakes Surface Environmental Analysis’’ product of the NOAACoastWatchnodeatGLERL.Dailyaveragedlakewidesurfacetemperaturesfrom1995to2013were availableandthelong-termaveragedtemperatures(1992–2013)wereusedforyearspriorto1995.Gridded air temperatures were from the CFSR data archive. An extrapolation technique similar to the technique usedforgriddedwinddatawasemployedforairtemperaturesattheland-seainterface. 3.1.2.IceField Icethatcanimpedewaveformationandtransmissionwasconsideredinthewavemodelhindcast.Theinter- actionbetweenwavesandiceisacomplexprocessthatshouldbeaddressedinwavehindcast.Inthisstudy, weemployeda commontechniquebymaskingmodelgridswithlandwheniceconcentrationexceededa particularthreshold.Iceconcentrationistheareaoficecoverdividedbythetotalareaofwater.Previousstud- ieshaveimplementedarangeoficeconcentrationthresholdsrangingfrom30to50%[Hubertzetal.,1991; Bennington et al., 2010; Tuomi et al., 2011]. Furthermore we implemented a land mask by lowering water depthstozerowheniceconcentrationsexceededathresholdof30%undertheassumptionthatwavescan be consideredunaffected by icefor ice concentrations 30% andbelow[Tuomi etal.,2011]. Great Lakes ice data prior to 2003 have been digitized to gridded fields, which are available at http://www.glerl.noaa.gov/ ANDERSONETAL. WAVECLIMATOLOGYINLAKESUPERIOR 4872 Journal of Geophysical Research: Oceans 10.1002/2014JC010278 Figure2.ScatterdensityplotsofCFSRandobservedU10(m/s)atselectlocations.Colorscalingislogarithmicandrepresentsthepercentageofhourlyrecordsoccurringwithina0.53 0.5m/ssquare.AlsoshownareQ-QplotsforrawCFSRandobservedU10(solidcircle).Atstation(c)SXHW3,aQ-QplotfortheextrapolatedoverseaCFSRwindsandobservedwindsis shown(plus).Theverticalandhorizontaldashedlinesrepresent99.9%quantileforeachsource.Thethickgreylineisthelineofequivalence.Collectionperiodsfortheavailableobserva- tionsareshownaboveeachplot. data/ice/atlas[Assel,2003].After2003,griddedfieldsoficeconcentrationhavebeenproducedbytheNOAA/ National Ice Center and can be obtained online (http://www.natice.noaa.gov/products/great_lakes.html). Griddediceconcentrationsrangedfrom0to99%andthespatialresolutionrangedfrom1.25to2.5km.Ice concentrationdatawerelinearlyinterpolatedintimeto adailyresolution [Assel,2005].Spatialinterpolation fromgriddedicedatatomodelgridswasaccomplishedbythenearestneighbortechnique. 3.1.3.WaveMeasurements Offshorewavedataatthethreebuoys, Stations(St.) 45001, 45004, and45006 (see Figure1), weredown- loadedfromtheNOAANationalDataBuoyCenter.Datesofthemeasurementsvariedslightlyforeachbuoy but nearly covered the full duration of the model hindcast (Table1). Each buoy measured SWH and PWP duringtheircollectionperiods,butonlySt.45001beganmeasuringMWDin2004.Furtherwavemeasure- mentswithintheApostleIslandswerecollectedatthreesites(Figure1)atvarioustimesbetween2011and 2013(Table1).AttheGSandOakShoal(OS)sites,a1MHzNortekAcousticWaveandCurrent(AWAC),man- ufacturedbyNortekAS,wasdeployedatapproximately12mwaterdepth.TheAWACusesacombination ofacousticsurfacetrack,pressuresensor,andacousticvelocimetrytomeasureSWH,PWP,andMWD[Nortek AS, 2005]. At the sea cave (SC) site, wave measurements were obtained using a single Acculevel pressure sensor, manufactured by Keller America, at 1.8 m above the bed in a water depth of 4 m. Pressure data wereconvertedtobulkwavestatisticswithapplicationofthelineartransferfunctionandspectralanalysis ANDERSONETAL. WAVECLIMATOLOGYINLAKESUPERIOR 4873 Journal of Geophysical Research: Oceans 10.1002/2014JC010278 Table1.SWANModelHindcastValidationStatisticsatAvailableLocations SWH PWP Direction Site RecordDuration Bias(m) RMSE(m) SI R Bias(s) RMSE(s) SI R Bias(8) RMSE(8) SC 5Jun2011to11Oct2013 20.02 0.12 0.47 0.85 20.27 1.20 0.35 0.63 OS 2Jul2013to1Oct2013 20.03 0.10 0.44 0.77 20.12 0.90 0.34 0.60 9.6 41.2 GS 11Oct2011to14Jan2012 0.06 0.18 0.26 0.92 0.04 1.00 0.25 0.76 3.2 34.8 St.45006 22Jun1981to30Nov2013 20.02 0.27 0.42 0.88 20.25 0.99 0.25 0.76 St.45001 3May1979to31Dec2013 20.02 0.30 0.37 0.90 20.19 0.80 0.18 0.82 7.7a 32.0a St.45004 26Apr1980to23Apr2013 0.04 0.29 0.37 0.91 20.08 0.79 0.18 0.81 ST95231 1Jan1979to31Dec2012 0.04 0.15 0.38 0.94 0.22 0.83 0.25 0.85 2.5 35.7 ST95216 1Jan1979to31Dec2012 0.02 0.15 0.29 0.96 0.20 0.75 0.21 0.89 21.5 28.7 ST95205 1Jan1979to31Dec2012 0.01 0.14 0.53 0.88 0.11 1.19 0.38 0.77 0.1 43.3 aDirectionaldatalimitedto21April2004to31December2013. [Sorensen, 2006; Jones and Monismith, 2007], which has been shown to produce wave height estimates accuratetowithin5%[BishopandDonelan,1987].MWDestimatescouldnotbedeterminedfromthepres- suresensor,butMWDismostlikelyperpendiculartotheshoreduetorefractiveeffects.Timeseriesplotsof themeasured SWHs and PWPsduring the most energetic conditions foreach site will be shownand dis- cussedinsection3.2.2. 3.2.Modeling 3.2.1.WaveModel The third-generation spectral wave model, Simulating WAves Nearshore (SWAN), was used to estimate waveparametersinthisstudy.AfullyimplicitnumericalschemewasimplementedtotheSWANmodelto remove the Courant stability criterion in explicit time stepping schemes [Booij et al., 1999]. To further improvetheefficiency,recentmodificationsincludethedevelopmentofaparallelizedcodestructureand unstructuredgrids[Zijlema,2010].SWANhasbeensuccessfullyappliedtomanycoastallocationstoexam- inewaveclimate[Gormanetal.,2003;Panchangetal.,2008;Watersetal.,2009;Stopaetal.,2013b]. Inthisstudy,weconstructedanunstructuredmeshforthewholeofLakeSuperiorwithhigherresolutionin theApostleIslands.Themeshconsistedofapproximately27,000triangularelementsrangingfromacharac- teristiclengthof5kmthroughoutthemajorityofthelakeareato100–500mintheApostleIslandsregion (Figure1).Discretizationinspectralspacewassetupwith50logarithmicallyspacedfrequencybandsrang- ing from 0.05 to 5 Hz and 36 evenly spaced directional bands of 108. The model time step was set to 10 min.Allavailablesourcetermformulationsandtunableparametersweresettotheirdefaultvalues[SWAN Team, 2013] except for wind growth and whitecap dissipation, which were set to the formulations com- monly referred to as WAM cycle 4 [Booij et al., 1999]. Currents and fluctuating water levels were omitted duetotheuncertaintyinthewindforcing,whichismostcommonlythedominantsourceoferrorinwave modeling [Jensen et al., 2012]. To implement a land mask for ice cover, water surface elevations were set belowbedelevationswhereicewaspresent.Modelgridswithzerowaterdepthwereeffectivelyremoved fromthecomputationaldomainofthewavemodelandbehavedsimilarlytolandboundaries. 3.2.2.CalibrationandValidation TheSWANmodelwascalibratedbyvaryingtwoparametersinthewhitecapdissipationformulation[SWAN Team,2013].Theoverallrateofwhitecappingdissipation,C ,wasvariedwithincrementsof0.5from1.0to ds 6.0,andtheparameterwhichdeterminesthedependencyofwhitecappingbasedonwavenumber,d,was varied by increments of 0.2 from 0.0 to 1.0. Evaluation of the results was based on minimizing the root- mean-squareerror(RMSE)betweenthemodelandmeasuredSWHattheSCandGSsitesfortheperiod10– qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 31October2011.TheRMSEisdefinedbyRMSE5 1 PN ðM2OÞ2,whereNisthenumberofsamples,M N i51 i i is model data, and O is observation data. Numerical tests revealed that both sites exhibited a similar response in RMSE to variations in the calibration parameters and that RMSE quantities were minimized whenC 53.5andd51.0.DefaultvaluesfortheSWANmodelareC 54.5,d50.5;however,itisrecom- ds ds mendedtoused51.0forbetteragreementofPWPsduringlowenergyeventswhileretuningC [SWAN ds Team,2013].AnapplicationofSWANintheGulfofMexicofoundC 52.0andd50.7providedthebestfit ds [Siadatmousavietal.,2011]. ANDERSONETAL. WAVECLIMATOLOGYINLAKESUPERIOR 4874 Journal of Geophysical Research: Oceans 10.1002/2014JC010278 Figure3.ScatterdensityplotsoftheSWANandWISmodeledSWHs(m)atthreelocations(seeFigure1).Colorscalingislogarithmicandrepresentsthepercentageofhourlyrecords occurringwithina0.130.1m2.AlsoshownareQ-QplotsforSWANandWISSWHs(solidcircle).Theverticalandhorizontaldashedlinesrepresentthe99.9%quantileforeachsource. Thethickgreylineisthelineofequivalence Validation of model performance was based on statistical quantities commonly applied in wave modeling studies[Gormanetal.,2003;Stopaetal.,2013b].Specifically,thebias51 PN ðM2OÞ,RMSE,scatterindex N i51 i i (SI)5 RMSE ,andcorrelationcoefficients(R)forSWH,PWP,andMWDwereevaluated.ForMWD,alimitof N1 PNi51Oi 1808wasappliedtoerrorestimates(M-O)duetothecircularnatureoftheproperty.Validationstatisticswere computedwithtwosourcesofdata.First,wecomparedtheSWANmodelresultstotheUSACEWIShindcast results, which were computed with the WAM Cycle 4.5 wave model [Jensen et al., 2012]. In this case, WIS modelresultsreplacedobservedvalues(O)inthecomputationofthestatistics.Table1showsthevalidation statisticsatthreeWISstations(ST95231,ST95216,andST95205),locatedaroundtheperimeteroftheApostle Islands(seeFigure1).Inaddition,weexaminedthescatterdensityandQ-QplotsofSWHtocomparemodel resultsforextremeevents(Figure3).Whiledifferencesexistbetweenthetwohindcastresults,thetwodata setscorrelatewell.Second,wecomparedtheSWANmodelhindcastresultstoavailableinsitumeasurements. Table1summarizestheresultsforthethreesitesintheApostleIslandregionandthethreeNOAAbuoys.In general,thewavemodelperformedreasonablywellinpredictingSWHwithlowbiasesandRMSEs.SIvalues for the SWHs are slightly higher than those reported in other studies [Stopa et al., 2013; Rusu and Guedes Soares, 2012],but this maybe explained bythecomparatively smallmean observed SWHsfrom this study, whichisknowntoinflateSIvalues[Risetal.,1999].ValidationstatisticsforSWHswerecomparablebetween measurementsitesintheApostleIslandsandatopenlakebuoyswiththeexceptionofslightlybettercorrela- tioncoefficientsatthebuoys.Duetoambiguityinthespectralpeakandwavedirectionduringlowenergy conditions,estimatesforPWPandMWDwerelimitedtotimeswhenSWH>0.15m.ModeledPWPswerebet- teratoffshorelocationswithlowerSIvaluesbyasmuchasafactorof2.Validationstatisticsfortheoffshore buoyswerealsocomputedoverthesametimeperiodasthenearshoreApostleIslanddataandasimilartrend wasfound.ThenegativebiasinPWPisduetoawell-knownlimitationofSWANduringlowenergyconditions [SWANTeam,2013].Anindicationofbettermodelperformanceduringmoreenergeticconditionsisseenfor theGSsite,whichwasdeployedduringthefallandrarelyexperiencedcalmconditions.Astatisticaltrendsim- ilartoPWPisseenforMWD.Overall,themodelperformanceisconsistentwiththeothermodelingstudiesin thatthevalidationstatisticsofwavemodelingexhibitslowerqualityinshelteredandnearshorelocationsdue to the more complex wave processes [Rusu et al., 2008]. In addition, the accuracy at offshore locations increasessincewavecharacteristicslikePWPandMWDarebetterdefinedatlargervalues.Forinstance,when SWH>0.5m,wavedirectionRMSEvaluesdecreaseto16.88and23.28fortheOSandGSsites,respectively. Figure 4 shows time series plots for SWH at the three measurement sites inside theApostle Islands. Each timeseriesspans60daysofmeasurementsandcontainsstormswithdominantwindsfromtheNEandNW (windsnotshownhere).Modelwaveheightsmatchwellwithobservationsforboththepeakandduration ofstormevents.Figure5showstimeseriesplotsofmodelandmeasuredPWPforthesamedatesshownin ANDERSONETAL. WAVECLIMATOLOGYINLAKESUPERIOR 4875 Journal of Geophysical Research: Oceans 10.1002/2014JC010278 Figure4.SWHcomparisonsofmodel(solidline)andmeasurements(solidcircle)collectedwithintheApostleIslands.Duetothedeploy- mentschedule,collectionyearsvarybetweensites:2011forSC,2013forOS,and2011–2012forGS. Figure 4. ObservedPWP measurements were excluded if SWH<0.15 m. Overall, the SWAN model does a reasonablejobcapturingPWPmagnitudes,especiallyduringenergeticconditions.However,afewinstan- cesofrelativelylargeerrorsdoexist(Figure5aon5October,Figure5con9November),whichmaypartially explainthelargererrorstatisticsatsiteswithintheislands(Table1).Shortinstancesofswellarewellpre- dictedatthetailendofstormevents,whichcanbeseenmostvisiblyinFigure5aon1October.Themodel performsconsistentlyinpredictingstormeventswithvariouswinddirections,demonstratingthereliability andaccuracyofthemodelinsimulatingwindwavesthroughouttheApostleIslands. 3.3.Analysis 3.3.1.WaveStatistics Statistics for the wave climate were based on the hourly model output over the duration of the hindcast withtheomissionofperiodsofice.Effectsofomittingicecoveredperiodsinthecalculationofwaveclimate statistics[Tuomietal.,2011]willbeaddressedinthediscussionsection.Inthisstudy,wavestatisticsinclud- ingmean,0.90and0.99quantiles,andmaximumvaluesofSWH,PWP,andMWDwerecalculated.Thedeci- mal of the quantile specifies the fraction of the data that do not exceed the quantile value. Unlike calculatingstatisticsofSWHfromtheentirehindcastrecord,statisticsofPWPandMWDwerebasedoncon- ditionspresentduringequivalentstatisticsofSWH.ThisdefinitionwaschosenbecausethemaximumPWPs inthehindcastoccurredasintermittentswellduringsmallSWHs,e.g.,seeFigures4and5fortheGSsiteon 31 December. So, the maximum PWP and MWD correspond to the values present during the maximum ANDERSONETAL. WAVECLIMATOLOGYINLAKESUPERIOR 4876 Journal of Geophysical Research: Oceans 10.1002/2014JC010278 Figure5.PWPcomparisonsofmodel(solidline)andmeasurements(solidcircle)collectedwithintheApostleIslands.EstimatesofPWP duringlowwaves(H <0.15m)wereexcludedfromanalysis.YearsoftheplotsvarybetweensitesandarethesameasFigure4. S SWH.SincearangeofPWPsandMWDsoccurredforagivenquantileofSWH,arepresentativevaluewas calculated by averaging PWPs or MWDs when the SWH was within 61 cm of the quantile value. This methodeffectivelydeterminesthePWPsandMWDsthataremostprobableforagivenquantileofSWH. 3.3.2.ExtremeValue EstimatesofSWHatvariousreturnperiods(i.e.,1,10,and100years)werecalculatedateachelement(loca- tion)inthemodel.Thepeaks-over-threshold(POT)methodwasappliedtoselectindependenteventsfrom thecontinuoustimeseriesofSWHs.Consecutiveexceedancesofthethreshold(i.e.,clusters)wereconsid- eredasasingleindependenteventfromwhichonlythemaximumvaluewasretained.Duetothevariability ofwindsinasinglestorm,werequiredthetemporalboundariesofconsecutiveclusterstobeseparatedby at least 48 h [Caires and Sterl, 2005; Aarnes et al., 2012] to further ensure only independent events were selected.DatasetsselectedbythePOTmethodwereassumedtobeindependentlydistributedwithagen- eralizedParetodistribution(GPD)withacumulativedistributionfunctiongivenby 8 >>><12(cid:3)11nrx(cid:4)21=n if n6¼0 FðxÞ5 (1) >>>:12exp(cid:5)2x(cid:6) if n50; r wherexvaluesarethePOTselecteddatasetminusthethresholdandtherangeofxisð0;1Þifn(cid:4)0and ð0;r=nÞifn>0[Coles,2001].Theparametersnandrarereferredtoastheshapeandscale,respectively. ANDERSONETAL. WAVECLIMATOLOGYINLAKESUPERIOR 4877 Journal of Geophysical Research: Oceans 10.1002/2014JC010278 Shape and scale parameter estimates were fitted using the maximum likelihood approach. Evaluation of the fits was accomplished based upon the Anderson-Darling goodness-of-fit-test at the 5% significance level[ChoulakianandStephens,2001].ThetestallowedustodeterminethevalidityofthefittedGPDand appropriatethresholdvalues. Inthisstudy,weusedthe95%quantileofthecompleteSWHtimeseriesforthethreshold,whichyieldeda valid fit of the GPD at 85% of model grids. For locations with invalid fits using the 95% quantile as the threshold,thequantilepercentagewasincrementallyincreasedupto99.5%untilavalidfitwasachieved. Thequantilepercentagesappliedinsettingthethresholdsherearewithintherangeusedinotherstudies: 93–97%[Caires andSterl,2005] and99.7%[Aarneset al.,2012]. Afterdetermining thethreshold value, we estimatedthereturnvalueas 8 rh i <l1 ðkNÞn21 if n6¼0 x 5 n (2) N :l1rlogðkNÞ if n50; wherelisthethreshold,kistheaveragenumberofexceedanceeventsperyear,andNisthereturnperiod inyears. 3.3.3.Long-TermTrends Tocalculatelong-termtrendsoverthe35yearhindcast,theleastsquareslinearregressionfitwasapplied toannualstatistics.Duetothewintericeseason,weappliedthecommonpracticetocalculateannualsta- tisticsbasedupontheclimateyear,definedasbeginningandendingonJunefirst[Schwabetal.,2006].The Mann-Kendalltest[McLeodetal.,1990;WangandSwail,2002]wasusedtoassessthesignificanceofatrend andtrendswereconsideredstatisticallysignificantwithapvaluelessthan0.05,i.e.,95%confidence.Long- termtrendsofmultiplevariableswereexaminedinthisstudy:meanSWH,0.90and0.99quantilesofSWH, icecoverduration,meanwindspeed,andmeanwinddirection. 3.3.4.NormalizedCumulativeWaveImpactHeight Coastalbluffrecessionoccurswhenwaveactionerodestheblufftoe.Acommonindexforbluffrecessionin the Great Lakes is thewave impact height (WIH), definedas the difference between wave runup and the blufftoeelevation: WIH5MWL1R1S2TOE (3) whereMWListhemeanwaterlevelelevation,Risthewaverunupheight,Sisthewindsetup,andTOEis theelevationoftheblufftoe[Brownetal.,2005].Swensonetal.[2006]proposedcumulativewaveimpact height(CWIH),whichisthepositiveWIHsintegratedoveratimeperiod,toaccountforthefrequency,dura- tion,andmagnitudeofallwaveeventsthatacttoerodeablufftoe. Inthisstudy,weexaminethetrendsinCWIHatsixsites,asseenassolidsquaresinFigure1.Thethreesites onthemainland(B1–B3)arecolocatedwithonesanalyzedinSwensonetal.[2006],andthethreesitesare locatedontheislands(B4–B6)identifiedashavinghighcoastalchangepotentialbyPendletonetal.[2007]. SitesB5andB6havehistoricallighthousesatopthebluff.ToestimatetheWIHateachlocation,equation(3) was applied. Mean water levels of Lake Superior were obtained at a monthly interval from the USACE- Detroit District over the hindcast period. Wave runup was calculated using the modified Mase method [Melby,2012]withdeepwaterwaveheightsandpeakperiodstakenfrommodelgrids(cid:3)1kmoffshoreofa site.BeachslopesandtoeelevationsusedinrunupcalculationswereacquiredfromSwensonetal.[2006] forsitesB1–B3andestimatedfromphotographsoftheshoresduringvisitstositesB4–B6.Windsetupand setdownwereestimatedwithalinearforcebalanceapproach[Sorensen,2006].TheWIHwasestimatedat each1houtputfromthehindcast.TheCWIHwastherebycalculatedbyintegratingovereachclimateyear andnormalizingbythenumberofdaysinthatyear.Asimilarlong-termtrendanalysiswasperformedon theCWIHsoverthehindcastperiodasdescribedinsection3.3.3. 4.Results Inthissection,theresultsfromthe35yearhindcastarepresented.Contourplotsweregeneratedtoillus- trate the spatial variability of the results within the Apostle Islands. To improve readability, contour plots weresmoothedusingamedianfilterovera1km2area. ANDERSONETAL. WAVECLIMATOLOGYINLAKESUPERIOR 4878