Mapping recreational visits and values of European National Parks by combining statistical modelling and unit value transfer

(1)

ContentslistsavailableatScienceDirect

Journal

for

Nature

Conservation

jou rn a l h o m e p a g e :w w w . e l s e v i e r . d e / j n c

Mapping

recreational

visits

and

values

of

European

National

Parks

by

combining

statistical

modelling

and

unit

value

transfer

Jan

Philipp

Schägner

a,∗

_,

_Luke

_Brander

b

_,

_Joachim

_Maes

a

_,

_Maria

_Luisa

_Paracchini

a

_,

Volkmar

Hartje

c

a_Joint_Research_Centre,_European_Commission,_Italy

b_Institute_for_{Environmental}_Studies,_VU_University_Amsterdam,_The_Netherlands c_Chair_of_{Environmental}_and_Land_Economics,_Technical_University_Berlin,_Germany

a

r

t

i

c

l

e

i

n

f

o

Articlehistory:

Received2September2015

Receivedinrevisedform1March2016 Accepted2March2016

Keywords:

Ecosystemservicemodelling Recreationaldemandmodelling Ecosystemservicemapping Ecosystemservicevaluation Recreationalvisitornumbers Protectedarea

a

b

s

t

r

a

c

t

Recreationisamajorecosystemserviceandanimportantco-beneﬁtofnatureconservation.The recre-ationalvalueofNationalParks(NPs)canbeastrongargumentinfavourofallocatingresourcesfor preservingandcreatingNPsworldwide.ManagingNPstooptimizerecreationalservicescantherefore indirectlycontributetonatureconservationandbiodiversityprotection.Understandingthedriversof recreationaluseofnationalparksiscrucial.

Inthisstudyweuseacombinationofprimarydataonannualvisitorcountsfor205EuropeanNPs,GIS andstatisticalregressiontechniquestoanalysehowcharacteristicsofNPsandtheirsurroundings inﬂu-encetotalannualrecreationalvisitornumbers.Thestatisticalmodelcanbeusedforland-useplanning byassessingtheimpactofalternativeconservationscenariosonrecreationaluseinNPs.Therecreational useofnewNPscanbeestimatedex-ante,therebyaidingtheoptimisationoftheirlocationanddesign.

Weapplythemodelto:(1)maprecreationalvisitstopotentialnewNPsacrossEuropeinorderto identifybestNPlocation;(2)maprecreationalvisitstoaproposednewNPinthewestofGermanyin orderestimatemonetaryvaluesandtoshowhowvisitsaredistributedacrossthesite;and(3)predict annualvisitstoallNPsof26Europeancountries.Totalannualvisitsamounttomorethan2billion annually.Assumingameanvaluepervisitderivedfrom244primaryvalueestimatesindicatesthatthese visitsresultinaconsumersurplusofapproximatelyD 14.5billionannually.

1. Introduction

NationalParks(NPs)areprotectedareasfortheconservationof extraordinarylandscapeandwildlifeforposterityandasasymbol ofnationalpride.NPscontributetostoppingthelossofbiodiversity, maintainingthenaturalnessandbeautyofourlandscapeandthe supplyofecosystemservices.Thereby,NPscontributeto achiev-ingthetargetsdeﬁnedinEUbiodiversitystrategy2020,suchas “haltingthelossofbiodiversityandthedegradationofecosystem services”(EC,2011), and theAichitargets,suchas“toimprove thestatusofbiodiversitybysafeguardingecosystems,speciesand geneticdiversity”(CBD,2013).

However,ﬁnancialresourcesandpolitical supportfornature conservation are limited and halting ecosystem degradation remainsagreatchallenge.Inthepast,majorpolicygoalson

biodi-∗ Correspondingauthor.

E-mailaddress:philipp.schaegner@gmx.net(J.P.Schägner).

versityprotectionhavetypicallynotbeenmet,suchasthosesetby theConventiononBiologicalDiversity,ratifiedaftertheglobal sum-mitinRiodeJaneiro(1992)(Barbault,2011;Leadleyetal.,2010). Andstill,thefutureoutlookrevealsthatbiodiversityremainsunder threatandsubstantialactionneedstobeundertaken(SCBD,2014). Onemajorco-benefitofnatureconservationisthesupplyof recreationalopportunities.NPsprovideopportunitiesforvisiting, experiencing, enjoyingandlearning aboutnatureand biodiver-sity,andthuscontributetohumanwell-beingandenvironmental awareness. Nature recreationand tourismpresent a great eco-nomicvalueandanopportunityforruraleconomicdevelopment bygeneratingincomeandemploymentthroughvisitors’ expendi-tures.Thevalueofnaturerecreationanditseconomicopportunities canbeusedasastrongargumentinfavourofallocatingfinancial resourcestowardsnatureconservationatdifferentspatialscales (Balmfordetal.,2015).

Natureconservationshouldnotonlyfocusonbiodiversityand habitatprotection,butshouldalsotakerecreationalco-beneﬁts into account. Efﬁcient land-use planning needs to consider all http://dx.doi.org/10.1016/j.jnc.2016.03.001

(2)

ecosystemservicessupplied.For allocatingresourcesfornature conservation, it can be important to know how recreational co-beneﬁtsof nature conservationcan beoptimized. The most importantindicatorofthecontributionofrecreationtothelocal economyisthenumberofvisitors(Jones,Bateman,&Wright,2003; Bateman,Day,Georgiou,&Lake,2006).Therefore,understanding thedriversthatdeterminethenumberofvisitorstoprotectedareas iscrucialforprotectedareamanagementandforprotectedarea designation.

TheaimofthisstudyistoanalysetheeffectsofNPcharacteristics andtheirspatialcontextontotalannualvisitsthatareconsidered themaindeterminantofrecreationaleconomic value(Bateman etal.,2006).Tothisend,we developregressionmodelsof visi-tornumbersusingprimarydataforEuropeanNPscombinedwith additionalspatialvariablesderivedfromGISdata.Theestimated modelsgiveinsightsintothedrivers ofrecreationalusewithin EuropeanNPsandthusallowthepredictionofvisitornumbersfor designatednewNPsandalternativemanagementscenarios. Simi-lartothestudyofBalmfordetal.(2015),wecombineourpredicted visitornumberswithameanvalueestimateperrecreationalvisit, butderivedfromamuchlargersetofprimaryvaluationstudies. Thereby,therelativeimportanceofrecreationalservicesis high-lightedascomparedtootherecosystemservicesandman-made goods.

Severalstudies have modelled visitor numbersof protected areasornatureareasbasedonspatialvariables.Onewidelyapplied approachistousechoicemodelstopredictrecreationalbehaviour at the individual level. Typically, such studies use survey data containinginformationontheoriginanddestinationofan individ-ualrecreationaltrip.However,suchdatasetsaretime-consuming to develop and are usually only available for relatively small areas(Pouta&Ovaskainen,2006;Batemanetal.,2011;Hausman, Leonard,&McFadden,1995;Jones,Wright,Bateman,&Schaafsma, 2010;Loomis,1995;Featheretal.,1995;Parsons&Hauber,1998; Senetal., 2013; Shaw &Ozog,1999;Termansen,Zandersen,& McClean,2008).Thepurposeofthepresentstudyistoinvestigate thedeterminantsofrecreationaluseofNPsataEuropeanscaleand thereforeweusedatafromvisitormonitoringstudiesforNPsacross Europe.Someexisting studies have usedsimilar approaches in ordertoinvestigatedriversofrecreationalparkvisits.Forexample, Neuvonen,Pouta,Puustinen,andSievänen(2010)analyseeffects ofparkcharacteristicsonvisitationratesfor35FinnishNPs.Mills andWestover(1987)modelthevisitationratesfor121Californian StateParksusingfourpredictorsrepresentingparkcharacteristics andthedistancetothenearestpopulationagglomeration.Hanink andWhite(1999)modelrecreationaldemandfor36USNational Parksusingageandsizeasvariablesfordescribingthepark,its dis-tanceandthepopulationoftheclosestmetropolitanarea,aswellas substituteavailabilityascontextcharacteristics.HaninkandStutts (2002)model thedemandfor19recreationalbattlefieldsinthe US.Theyuseasubstituteavailabilityindicatorweightedby indi-vidualsubstitute’scharacteristics.Loomis,Bonetti,andEchohawk (1999)findasignificanteffectofGDPpercapitaandofavailability ofwildernessonthenumber ofrecreationaltripstowilderness areasper capitaintheUS. Ejstrud(2006) useanumber of GIS indicatorsformodellingvisitorfrequencyto10Danishopen-air museumsusingsixpredictorvariables,butdonotreportwhether theyshowsignificanteffects.Theonlystudyusinginternational visitordataisfromBalmfordetal.(2015),whichusesvisitordata ofprotectedareasworldwide.Theirstudyusesonlyalimited num-berofrelativelysimplepredictorvariablesandfindsfewsignificant effects.Theirmodelmaybeappropriatetoassessoveralltrends inprotectedareavisitationrates,butmayhavefewsitespecific implications.Loomis(2004)usesregressiontechniquestoestimate theeffectofelkandbisonpopulationsonvisitationratesinGrand TetonNationalPark,US,usingexplanatoryvariablesonhowthe

parkchangesovertime,butdoesnotcompareeffectsofalternative sites’characteristics.

Allexceptoneoftheabovementionedstudiesusenationaldata onlyfortheirstatisticalanalysis.Thereby,thenumberofprimary observationsisingeneralrelativelylow.Thepurposeofthepresent studyistoinvestigatedriversofrecreationaluseforNPs Europe-wide and therefore, use visitordata from NPs in 21 European countriescomprising205casestudyareasintotal.Consequently, wecan includemore predictors inourinitialmodel and tryto estimateamorerobustmodel.Forexample,nationalstudyareas arerelativelysmallandthereforeclimaticconditionsareoftentoo similartobeconsideredasapredictorinarecreationaldemand model.Furthermore,weusemorerefinedsiteandcontext char-acteristicsaspredictors inourmodel,whichare computedand extractedfromEurope-wideGISdatalayers.Asallourpredictors arederivedfromlargescaleGISdatalayers,thefinalmodelcan easilybeusedtomakepredictionsofvisitors’frequencyforany potentialNPinEurope.Thus,recreationalusecanbemappedfor anylocationinEuropewithouttheneedforanadditionalcollection ofinformationonthepredictorvalues.Ourspatialassessmentcan therebybeusedforecosystemservicemappingasrequiredbythe EUBiodiversityStrategy2020,improvingresourceallocationand calculatingagreenGDP(UN,2014;Maesetal.,2012).Finally,we useanumberofdifferentstatisticalregressiontechniquestodeal withspatialautocorrelationforamorein-depthidentificationof thespatialdimensionofrecreationaluse.

Thispaperisorganizedasfollows:insectiontwowedescribe thedataweuse,ﬁrsttheprimarydataofvisitormonitoringstudies andsecondthepredictorsusedinourmodels.Insectionthreewe explainthestatisticalregressiontechniquesappliedandpresent theestimatedvisitormodels.Theresultsarepresentedand dis-cussedinsectionfourandﬁve,withconclusionsprovidedinsection six.

2. Data

2.1. Primarydata

Ourprimarydataare205totalannualvisitorestimatesto Euro-peanNPsand245estimatesofmonetaryvaluesperrecreational visit for147separatenature areasin Europe.We collectedthe datathroughinternetsearches,reviewofrelevantliteratureand bycontactingresearchersinvolvedinthisﬁeld,NPadministrations andrelevantgovernmentalbodiesinallEUcountries.Thedatais describedmoreindetailin(Schägneretal.,submitted).

Forthevisitordatatobeincluded,werequiredasaminimum qualitycriteriathatthetotalannualvisitorestimatesarebasedon someformofon-sitevisitormonitoring,whichisthenscaledup totheentireareaandtheentireyear.Inordertocheckwhether thequalitycriteriaismet,weanalysedtherelevantpublications onthevisitormonitoringprograms.Incasesinwhichthe informa-tionwasnotavailableornotaccessibleduetolanguagebarriers,we contactedtheauthorsandrelevantinstitutions.Intotalwecould obtainannualvisitorobservationsfor205separatecasestudyareas withinEurope,whichareeitheranentireNPorasubsectionofa NP(seeFig.1).Allcollecteddatawereattachedasattributestoa spatiallayerinvectorformat,containingtheboundariesofNPsor oftheirsurveyedpart.WeobtainedNPpolygonsfrom(WorldData BaseofProtectedAreas)andtheCDDA(CommonDatabaseon Des-ignatedAreas)(IUCN&UNEP,2015;EEA,2013)andfromnational agencies.Ifcasestudyareasdifferedfromtheavailablepolygons, wetriedtoobtainpolygonsfromtheauthorsofthestudies,the parkmanagementorotherstakeholders.Insomecaseswe manu-allydrawpolygonswithArcGIS,basedoninformationavailablein thecasestudypublicationsorsuppliedbytheauthors.Ifmultiple

(3)

Fig.1.LocationofvisitorcountsacrossEurope.

observationsofvisitornumbersareavailableforthesamestudy area,weusedtheaverage.

NPandcasestudyareacharacteristicsdifferwidelyinterms ofsize,location,visitationrateandecosystemcharacteristics.The smallestcasestudyareaisaninehectarebeachwithintheWadden SeaNPinGermany,whereasthelargestcasestudyareaisthe Cairn-gormsNPinScotlandcomprising3816km2_._Most_of_the_case_study

areasinourdatabasearelocatedinNorthernEurope.For South-ernEuropewecouldobtainvisitornumbersforallSpanish,most ItalianandFrenchNPs.Forourstatisticalanalysiswedividedthe totalannualvisitornumbersbythetotalterrestrialareaofthe sin-glestudyareas1_and_thereby_obtained_total_annual_visitor_densities perhaasourdependentvariableinourmodels.Visitornumbers rangefrom0.03visitors/ha/yearinthelargeSarekNPinnorthern Swedenupto56,680visitors/ha/yearonasmallbeachwithinthe WaddenSeaNP.Thetotalmedianandmeanis13and368with standarddeviationof3962visitors/ha/year,indicatinga skewed distributionwithatailofveryhighvisitationrates.Themean rela-tivedeviationisabout167%.Formoreinformationontheprimary data,itcanbeaccessedviatheESPVisualisationTool(Drakouetal., 2015).

Forourstatisticalanalysiswedividedthetotalannualvisitor numbersbythetotalterrestrialhasizeofthesinglestudyareas andtherebyobtainedtotalannualvisitordensitiesperhaasour dependentvariableinourmodels,whichiscommonwithinspecies distributionmodelling.

ThevaluationstudiesuseeitherTravelCostMethod(TCM)(57%) orContingentValuation Method(CVM)(43%).For thevaluation studies,wetransferallvalueestimatestoEuro2013pricelevel usingpurchasingpowerparityandcountryspeciﬁcinﬂationdata. Weexcludeoneoutlierwithanextremedeviationof60timesthe meanvalue.TheremainingvalueestimatesrangefromD 0.16to 64.7pervisitwithameanofD7.17,amedianofD2.8,astandard deviationof11andameanrelativedeviationof95%.Moststudy sitesarelocatedinWesternEurope(51%).TheUKhasthehighest

1_We_used_the_terrestrial_area_not_including_area_covered_with_water_because_some

NP–inparticularmarineNP–comprisemainlyofwater.Includingtheareaofwater wouldbiasouranalysissincethisareaishardlyvisited.

numberofobservations(81),followedbyItaly(32),Ireland(28), Finland(27)andGermany(22).

2.2. Explanatoryvariables

Explanatory variables used tomodel visitationrates can be divided into three categories: (1) site characteristics, which describetheNPitself;(2)contextcharacteristics,whichdescribe thespatialcontextoftheNP;and,(3)studycharacteristics,which describethemethodologyofprimarydatacollection.Theselection ofvariableswasbasedonareviewoftheliteratureonrecreational demandmodellingandenvironmentalrecreationalvaluetransfer studies.However,limitationsintheavailabilityofcomprehensive andconsistentEurope-widedatasetsandintheinformation pro-videdinvisitormonitoringpublications restrictedourchoiceof predictors.A completelistofallpredictors usedinouranalysis ispresented inTable1.Detaileddescriptionispresentedin the followingsections.Eachvariableisavailableingeospatialraster format,thereforesiteandcontextcharacteristicsforeachsitecould beeasilycalculatedinaGISenvironment.Weextractedmean val-uesof allpredictorvariables for each case studyareausingan automatedmodelbuiltinArcGIS,includingtheuseofthezonal statisticstool(ArcGIS10.1).Therasterlayersofthepredictorswere eithertakenfromavailableGISdatasetsorwecomputedthemby reprocessingorcombiningexistingdatasetsusingArcGIS(ArcGIS 10.1).Thenweconductedanexplorationofourdatafollowingthe recommendationsof(Zuur,Ieno,&Elphick,2010)inordertogain initialinsightsintodistributionsanddependencies.Forsome pre-dictorsweusedlogarithmicorsquareroottransformationseither becausetheyshowedarelativelyskeweddistributionorbecause wewantedtoapproximatelylinearizeanexpectednon-linear rela-tionship. We tested all ourpredictors for multicollinearity, but couldnotidentifyanythingofconcern.

2.2.1. Sitecharacteristics

The following site characteristics are used to model visita-tionrates:(1)Shareoflandcover/use:WeusedtheCORINEland cover/usedataset(EEA,2006)todeterminethesharesof differ-entlandcover/useclassesandaggregatesofsinglelandcover/use classesforeach NP.Inparticularwefocusedonnatural vegeta-tioncover.We do not,however,havestrongpriorexpectations regarding the signs of these land cover predictors. In general, onemayassumethat naturalvegetationsupportsnature recre-ation.However,NPstypicallyofferplentifulnaturalvegetationand thereforeadditionalnaturalvegetationofanykindmaynot neces-sarilyattractadditionalvisitors.Ouranalysisofthedifferentland covershasanexploratorycharacteranddoesnotaimtotest spe-ciﬁc hypotheses. Theseparateclasses and aggregated areasare presentedinTable1.(2)Waterbodies:Wecomputeda300m res-olutiongridoftheshareofsurfaceareacoveredwithrivers,lakes oroceanusingtheEuroRegionalMapasinputdataset(EG,2010). Thenweapplieda kerneldensityfunctiontool (ArcGIS10.1)to computetheamountofsurfacecoveredwithwaterwithina3km radiusofeachpixel.Thedensityfunctionallowswaterareathatis furtherawaytobeweightedlessthanwaternearbyandthereby incorporatesadistancedecayeffect.Thepresenceofwaterbodies inaNPareexpectedtohaveapositiveimpactonrecreationaluse (Termansenetal.,2008).

Weexpectthatmorediverselandscapesareperceivedasmore beautiful(Dramstad&Tveit,2006)andtherebyattractmore visi-tors.Basedonthebasiceconomicprincipleofdecreasingmarginal utilityandratesofsubstitution,diversitytendstoberatedhigher than uniformity (Mankiw, 2001).In order toaccount for land-scapediversitywecomputedthreedifferentindicators.(3)Three dimensionality: We computedthe areavisible from each pixel within a 30km radius using the view shed tool (ArcGIS 10.1)

(4)

Table1

Listofpredictorsusedinthemodels.

Type Variables Explanationa _{Mean/Standard}_Deviation

SiteCharacteristics: Sqrt(grassland) Shareofgrasslandscoverofthestudyarea(100mresolution raster)

0.2/0.24 Sqrt(wetland) Shareofwetlandscoverofthestudyarea(100mresolutionraster) 0.14/0.23 Sqrt(water) Shareofwaterbodiesofthestudyarea(300mresolutionraster) 0.23/0.26 Log(broadleaf) Shareofbroadleafforestofthestudyarea(100mresolutionraster) 0.73/0.86 Conifer Shareofconiferforestofthestudyarea(100mresolutionraster) 4.44/4.63 Log(forestedge) Transitionareabetweenforestandotherlanduse/cover(25m

resolutionraster)

0.83/0.4 Sqrt(landcoverdiversity) SimpsonDiversityIndexofCorinelanduse/coverwithina3km

radius(100mresolutionraster)

1.61/0.22 Log(viewshed) Areavisiblefromeachlocationwithinina30kmradius(1km

resolutionraster)

5.43/0.69 Log(redlistspecies) Totalnumberofredlistspeciesfoundinstudyarea 2.65/0.84 Temperature Totalnumberofdayswithmaximumtemperatureabove5◦

Celsius(10kmresolutionraster)

256/57.5 NPage YearssinceNPfoundationuntil2015 40.6/26.94 Log(trails) Traildensityusingdensityfunctioninordertoaccountfor

distancedecayeffect

5.69/1.87 Log(roads) Densityofminorroadsusingdensityfunctioninordertoaccount

fordistancedecayeffect(100mresolutionraster)

0.9/0.83 Studyareakm2 _Size_of_the_study_area_in_km2 _352/621

ContextCharacteristics: Log(NPsubstitutes) AreaofNPwithin130kmradiusofthestudyareausingaGaussian weightfunctioninordertoaccountfordistancedecay(1km resolutionraster)

11.35/1.5

Log(Population50km2₎ _Population_living_within₅₀_km_radius_of_the_study_area_using_a

Gaussianweightfunctioninordertoaccountfordistancedecay (100mresolutionraster)

12.88/1.75

GDP/capita GDP/capitaintheNUTS2or3regioninwhichthestudyareais located

21,856/7713

StudyCharacteristics: Surveyyear Yearofvisitormonitoringsurvey 2005.6/4.16 Surveyquality Qualityofthevisitormonitoringsurveymethodologyandstudy

areadeﬁnition

7.17/1.53

a_For_all_predictors_mean_values_per_study_area_were_computed.

anda1000mresolutiondigitalelevationmapfromtheEuropean EnvironmentalAgency(EEA,2015a).Webelievethatvisitors pre-ferthree-dimensional landscapesofferinggreat views. (4)Land use/coverdiversity:BasedontheCORINElanduse/coverdataset wecomputedtheSimpsonDiversityIndex(Magurran,1988)ofland use/coverwithina3kmradiusforeachpixeloftheCORINEmap.In theirstudy(Neuvonenetal.,2010)usethenumberofbiotopesasa diversityindicatorandfindasignificantpositiveeffectonvisitation frequencyinFinnishNPs.However,thenumberofbiotopesmaybe positivelycorrelatedwiththestudyareasize.Therefore,this pre-dictormaypickuppartofthesizeeffect.Furthermore,largerNPs mayhavemorebiotopeseveniftheirlandscapeisnotmorediverse. (5)Forestedges:UsingtheJointResearchCentreforestcovermap (EC,2006),wecomputedthenumberofforestpixels(25m res-olution)thatarenotclassifiedasforestcore.We considerthese forestpixelsasthetransitionareabetweenforestandotherland use/coverandtherefore,asamajorvisiblechangeinthe ecosys-temtype(EC,2006).(6)Temperature:Weappliedadatasetfrom (Biavetti,Karetsos,Ceglar,Toreti,&Panagos,2014)indicatingthe numberofdayswithmaximumtemperatureabovefivedegrees Celsius.Duetothepredominanceofsouthboundtourismfluxesin Europe,weexpecttemperaturetohaveapositiveeffecton visi-tationrates.(7)Regions:Siteswerefurtherclassifiedaccordingto theirmembershipofbio-geographicalandgeographical regions. We do not have expectations regarding thesigns of these fac-torvariables,but mightdiscoversomeculturaleffects.(8) Trail density:Weusedtraildensityasproxyforoverallrecreational facil-ities,whichmayattractvisitors.FromtheOSM(OpenStreetMap) dataset(OSM,2012),weextractedallvectorelementsthatcanbe classifiedasnon-motorizedtrafficinfrastructure.WeusedfiveOSM classes:trails;footpaths;bikepaths;bridlepathsand,steps.On a100mresolutionweappliedthelinedensitytool(ArcGIS10.1)

tocomputeanindicatorfortrailavailability.Again,trailsthatare furtherawayfromapixelwereweightedlessthantrailscloseby. Otherstudiesfoundsignificantpositiveimpactsoftrails(Neuvonen etal.,2010)orrecreationalfacilitiesingeneral(Mills&Westover, 1987),buttheyusedindividualparkdataandnocomprehensive largescaleGISdatasets.(9)Streetdensity:Similartotraildensity wecomputedanindicatorforstreetavailabilityforallminorroads (TeleRoadAtlasroadclasses4–6)basedontheTeleRoadAtlas dataset(TS,2006).Roadsareanimportantinfrastructurefor access-ingremotelocationsandtherebyareexpectedtoincreasevisitor numbers.However,ifroadsaretooabundant,theymaynegatively affectthequalityperceptionofnaturerecreationinaNPandthus, detervisitors.(10)Studyareasize:Weexpectthatareasizehas anegativeimpactonthemeannumberofvisitorsperhabecause oftworeasons:First,largerstudyareasactasasubstituteinitself, becausevisitorscanbedistributedacrossalargerarea.Second, visi-torcountingtendstoresultinlowermeanvisitornumbersforlarger areas.Ifavisitorhikesthroughalargestudyarea,heiscounted once.Ifthesamestudyareaissplitintoseparatestudyareas,the samevisitormayeventuallybecountedseveraltimes.Most exist-ingstudiesofNPvisitsusetotalvisitornumbersasthedependent variableandthereforefindapositiveinfluenceofstudyareasize onvisitornumbers(Hanink&Stutts,2002;Hanink&White,1999; Mills&Westover1987).However,byworkingwithlinearmodels theypotentiallymissoutthatvisitornumbersdonotincreasein directunitaryproportiontothesizeofthestudyarea.(11)Ageof NP:Finally,wecharacterizedeachNPbyitsage(numberofyears sincefoundationuntil2015).Existingstudieshavefounda posi-tivecorrelationbetweenparkageandvisitornumbers(Neuvonen etal.,2010;Mills&Westover,1987;Hanink&Stutts,2002;Hanink &White,1999).Thismaybecausedbythegeneraltendencythatthe mostattractivelocationsweredesignatedasprotectedareasearlier

(5)

orthatolderNPshavehadmoretimetoestablishrecreational facil-ities.ThedesignationofaNPmaycreateanadvertisementeffect andestablisha goodreputationincreasingtheparkspopularity overtime.(12)Biodiversity:Inthiscaseweusedthetotal num-berofredlistspeciesencounteredinastudyareaasanindicator forbiodiversity(IUCN,2013).

2.2.2. Contextcharacteristics

Ascontextcharacteristicsweusedthefollowingvariables:(1) Accessibility:Weexpectthatthenumberofpeoplethatcanaccess acertainlocationwithinacertaintimeislikelytohaveapositive effectonthevisitationrate.Wedefinethis variableasthetotal populationlivingwithina50kmradiusaroundthesite,using pop-ulationdatafrom(BatistaeSilva,Gallego,&Lavalle,2013).Inorder toaccountfordistancedecay,weappliedaGaussianweight func-tion,whichcausesthepopulationthatisfurtherawayfromtheNP tobeweightedlessthanthepopulationnearby.Theweight func-tionwascalculatedsothat95%ofitsintegralwaslocatedwithin the50kmradius.Otherstudiesfindsignificantpositiveeffectsof accessibilityonvisitornumbers.Theyuseforexampledistanceto nearesttowns(Mills&Westover,1987)orconsiderthe popula-tionofmetropolitanareas(Hanink&Stutts,2002;Hanink&White, 1999)anddonotincludedistancedecayeffects(Neuvonenetal., 2010).(2) NPsubstitutes:We computedaraster inwhich each pixelis thesumofareasclassifiedasNP within130km radius. TheEurope-wideNPdatasetwasacombinationofsitesfromthe WDPAandCDDAdatabases.Inordertoaccountfordistancedecay, weusedthesamemethodologyasforpopulation.Asaresult,large NPsandNPswithsmalldistancefromeachotherhavearelatively highavailabilityofsubstitutes.Otherstudieshavefoundnegative influencesofsubstituteavailabilityonvisitornumbers.Theyuse forexampledistancestocompetingrecreationalsites(Hanink& White,1999;Hanink&Stutts,2002)orthenumberofparkswithin acertaindistance(Neuvonenetal.,2010).Theydonot,however, accountforthesizeofsubstituteareas.(3)Finally,weintroduce GDPpercapitaasaproxyofvisitorincome,whichweextracted fromtheEurostatdatabase(EC,2013).Wetookthemeanvalues ofthelasttenyears(asfarasavailable)andthehighestdata reso-lutionavailable,whichiseitherNUTS2orNUTS3level.Weexpect thatvisitationratesarelikelytobehigherinlocationswithhigher percapitaGDP.Existingstudieshaveobservedthatpeople engag-inginnaturerecreationhaveaboveaverageincomes(Loomisetal., 1999).

2.2.3. Studycharacteristics

Initially,weconsideredcollectingdetailedinformationonstudy characteristicsdescribingthemethodologyofthevisitor monitor-ingprocedureforeachcasestudyarea.Inthatway,wehopedto identifytheinfluenceofdifferentvisitormonitoringtechniqueson thefinaltotalannualvisitorestimate.Similarattemptshavebeen successfullyimplementedinmeta-analysisstudiesof environmen-taleconomicvaluationstudies(Zandersen&Tol,2009;Brouwer, Langford,Bateman,&Turner,1999).However,weencountered dif-ficultiesincodingsuchmethodologicalstudycharacteristicsdue tothelanguageandincompletereportingintheunderlyingcase studypublications.Therefore,weonlyintroducetwostudy char-acteristicsaspredictorsinouranalysis:(1)theyearofthevisitor monitoringsurveyforwhichweusedthemeanvaluesoftheyears inwhichvisitormonitoringtookplace.(2)Furthermore,we clas-sifiedallvisitormonitoringstudiesaccordingtodifferentlevelsof primarydatacollectionqualityfromoneforthelowestandtenfor thehighestquality.Thequalityjudgmentrepresentsacomposite indicatorofdifferentqualitydimensions:thetypeofpublication (scientificvs.grey literature);thevisitormonitoringstudy pur-pose(scientificvs.political);theinstitutionconductingthestudy (academic, NP management, others); themethodological

docu-Fig.2. Bubbleplotofthespatialdistributionofthefullmodel’sresidualwithout spatialcorrelationstructure.

mentationofstudy(full,incomplete,none).Ifthedocumentation forthestudywasavailable,weassessedthequalityof method-ologiesbasedondetailssuchasthetemporalandspatialcounting resolution,manualorelectroniccountingdevicesandthe tempo-ralandspatialup-scalingmethodology.Finally,averyimportant aspectforthevisitormonitoringstudiesqualityisthedescription ofthestudyarea.Somepublicationsdonotsupplymapsandonly roughdescriptionsofthestudyarea.Iftheareaofthestudyareais uncertain,thenthenumberofvisitorsperhectareisuncertainas well.

3. Methodology

Weappliedanumberorregressiontechniquesinordertomodel thetotal annualvisits perhatoEuropeanNPsusingtheabove describedpredictors. Allmodelswereestimatedusingtheopen sourcestatisticalsoftwareR.Westartedouranalysiswitha sim-plelinearregression,butitshowedastrongspreadoftheresiduals forlargerﬁttedvaluesandthereforeaviolationofthe homogene-ityassumption.We tried tocontrolthis effectbyintroducinga numberofdifferentvariancestructures,butwerenotsuccessful ineliminatingtheheterogeneitytoanacceptabledegree.

Asourdependentvariableisacount,wecontinuedour analy-siswithgeneralizedlinearmodelsusingaPoissonandanegative binomialdistribution(usingR-packageglmmADMB,MASS,lme4, nlme andgamlss(Bolkeretal., 2012;Ripleyet al.,2015; Bates etal.,2015;Pinheiroetal.,2015;Stasinopoulos,Rigby,Voudouris, Akantziliotou, &Enea, 2015), which are typicaldistributions of countdata(O’Hara&Kotze,2010).However,modelresultsshow spatialresidualpatternssimilartotheonedisplayedinFig.2.The negative(greybubbles)andpositiveresiduals(blackbubbles)are clustered,whichisaviolationoftheindependenceassumptionof generallinearregressionanalysis.Inordertoovercomethis prob-lemweaddedaspatialresidualstructure,eitherbyaspatialrandom effectoraspatialautocorrelation,butweranintonumerical con-versionproblemsoftheoptimizationalgorithmtryingtosolvethe complexstatisticalmodel.Wethereforeabandonedthisapproach anddonotpresenttheinterimresultsoftheseattempts.

Because ourcountdatashows relatively largevalues (mean value367),logtransformationisanalternativeapproach,which shouldhave a negligible effectontheparameterestimates but

(6)

decreasesthemodelprocessingcomplexitysubstantially(O’Hara &Kotze,2010).Wethereforecontinuedouranalysiswithlinearlog transformedmodelofthefollowingform:

log (Vi)=˛+ˇ∗Xi+i where i∼N(0,2)

Vstandsforthedependentvariable(inourcasethetotalannual visitsperha),␣isaconstant,␤representsavectorofparameters, Xisavectorofexplanatoryvariablesand␮istheresidual,whichis normallydistributedwithmeanofzeroandvariance␴.Again,we hadtodealwithspatialresidualpatterns,whichwetriedtocontrol forusingaspatialrandomeffectinamixedmodel2_and_by_a resid-ualspatialautocorrelationstructure.Wetriedanumberofdifferent randominterceptsandrandomslopesinthemixedmodelandalso anumberofspatialautocorrelationstructures.3 _We_investigated allestimatedmodelsonhowsuccessfultheyareincontrollingfor thespatialresidualpatternsandontheirAICandBICscores(as criteriaformodelselection).Thebestmodelcontainedaspatial sphericalcorrelationstructure,whichmodelstheresiduals’ cor-relationacrossspaceasphericalfunctionofdistance.Themodel formularemainsthesameasbefore,butthistimeweassumethat theresiduals␮iofdifferentlocationsarecorrelatedbasedonthe

functionfandtheirdistance. cor (a,b)=

1 ifa=b

f (a,b,) else

Weusedthismodelasastartingpointandconductedstepwise modelselectionbydropping theleastsignificantpredictoruntil everypredictorwassignificant.Wedeterminedstartingvaluesfor therange(maximumdistanceofspatialcorrelation)andthenugget (oneminusthecorrelationoftwoarbitrarilycloseobservations)of spatialcorrelationstructurebasedoninterpretationofvariogram andspatialresidualplotsinordertoimproveconsistencyacross thedifferentmodels.Inthefollowingsectiononresults,wepresent detailedresultsonourinitiallogtransformedmodel,thestarting modelincludingthespatialsphericalcorrelationstructureandon thefinalmodelafterstepwisemodelselection.Wevalidatedour finalmodelagainsttheassumptionsoflinearregressionanalysis. Therefore,weplottedourresidualagainstfittedvaluesandagainst eachpredictor.Wecouldnotidentifyanylinearornon-linear pat-ternsofconcern.Topresentacomparablemeasureofthegoodness offitofallmodelswecomputetherootmeansquaredeviation (RMSE)andthecoefficientofvariationoftheRMSD(CVRMSE).

Weuseourﬁnal model (1)tomake predictionsofthetotal annualvisitstoallEuropeanNPswithinthecountriescoveredby ourexplanatoryvariablelayers,(2)tomapthetotalannualvisitsto aﬁctivenew80km2_NP,_located_anywhere_in_the_European

coun-triescoveredbyourexplanatoryvariablelayersand(3)tomapthe distributionofthepredictedtotalannualvisitstoaproposednew NP(TeutoburgerforestandSenneheathland)inthewesternpart ofGermany.

InordertopredictthenumberofvisitstoNPsofmostEuropean countries,weextractedallshapeﬁlesfromtheWDPAandtheCDDA (IUCN&UNEP,2015;EEA,2013),whichfallintotheIUCNcategory II(NationalPark).Furthermore,weaccessednationaldatabasesto obtainshapesofNPs,whichweremissinginthosetwodatabases. Intotalweincluded449separateNPsareas.Itistobenotedthatnot allofthesesitesfallintoIUCNcategoryII.Nouniformdeﬁnitionof thetermNPsexistsanditwasusedlongbeforetheIUCNcategories systemwascreated.ManyexistingNPsallovertheworldare

dif-2 _In_other_disciplines,_mixed_modelling_is_also_referred_as_to_multilevel_analysis,

nesteddatamodels,hierarchicallinearmodels,andrepeatedmeasurements.

3 _For_an_introduction_into_mixed_modelling_we_would_like_refer_the_reader_to

(Zuur,Ieno,Walker,Saveliev,&Smith,2009)andforinintroductionintospatial autocorrelationto(Bivand,Pebesma,&Gómez-Rubio,2013).

ferentlymanagedthandemandedbytherequirementsofcategory II((InternationalUnionforConservationofNature)IUCN,2008), butarestillcalledNPbasedonthedecisionofgovernmentsand otherlocalstakeholders.WeusedvectorlayerofallNPs bound-ariesandzonalstatistics (ArcGIS10.2)todrivemeanvaluesfor theexplanatoryvariables.Predictionsweremadeusingtherms R-package(Harrell,2015).Inordertoimproveourpredictionsand accountforunobservedeffectsonvisitationrates,wekrigedthe residualsofourmodelacrosstheentirestudyareausingthegstat, GeoRandrasterR-packages(Pebesma&Graeler,2015;Hijmans etal.,2015;Diggle,Ribeiro,&Peter,2015).Wethenaddedtheresult tothepredictionofeachNP.

ForpredictingthenumberofvisitsofamarginalincreaseofNP area,weassumeaﬁctivelycreatedmediumsizeNPof80km2_._We

thencreatedexplanatoryvariablerasterlayersaccountingforthe averagesubstituteeffectofthenewNPandthesizeofthenewNP. Thequalityofthevisitormonitoringmethodology,whichisone explanatoryvariableinourmodel,wassettothehighestquality availableinourprimarydatabase(9.5).TheNPagewassettozero. Wethenusedthemodeltomaptheannualnumberofvisitsfor each1km2_resolution_grid_cell_across_Europe,_as_though_it_is_part_of

thenewlycreatedNP.Themappingwasconductedusingtheraster, gstatandgeoRR-packages.Again,weaddedthekrigedresidualsto ourpredictions.

Inordertotestourvisitormappingprocedureinarealisticpolicy setting,weappliedittoaproposednewNPinthewesternpart ofGermany(TeutoburgerforestandSenneheathland).Thearea ofthisproposedNPisapproximately20,000haandcomprisesa forestedmountainrangeandaheathland,whichhadbeenusedas anarmybaseinthepast.Itisalreadylargelyprotectedandhasbeen proposedforNPdesignation(NABU,2015).Wemadepredictions on1haresolutioninordertoestimatetotalvisitstotheareaand showhowvisitorsdistributeacrossthearea.

Finally,wecombinethepredictednumberofvisitswitha mon-etaryvalueestimate,derivedbytakingtheoverallmeanvalueper visit(7.17D)fromthe244valueestimatesdescribedabove,which isalmostthesamevalueestimateappliedina similarstudyby Balmfordetal.(2015)(7$),butbasedonmuchlargerprimary val-uationdatabase.This approach,so-called unit valuetransferor averagevaluetransferandisacommonapproachusedforvalue transferandecosystemservicevaluemapping(Schägner,Brander, Maes,Hartje,2013;Rosenberger&Loomis,2001;Balmfordetal., 2015)andamethodconsideredforaggregatingecosystemservice valuestodevelopaSystemofEnvironmental-Economic Account-ing(SEEA)(UN,2014).Itassumesaconstantvalueperrecreational visitacrossspace,which isindeedasimpliﬁcation.However,as thevalueperrecreationalvisitvariesbyfarlessacrossspacethan thenumberofrecreationalvisits(Batemanetal.,2006;Jonesetal., 2003),itseffectontheoverallrecreationalvalueofanareais rela-tivelysmall.GiventhefactthatwefocusonlyonNPsandonanarea ofrelativelysimilarsocioeconomicandculturalcharacteristics,we considerunitvaluetransferasagoodapproximationforthecase studypresented(seediscussionforfurtherdetails).

4. Results

The results of the statistical NP visitor model using a log-transformeddependentvariablearepresentedinTable2.14ofthe 19predictorsshowstaticallysigniﬁcantcoefﬁcientsandthe mul-tipleR2_of_0.68_indicates_relatively_high_explained_variance._Most

coefﬁcientshavetheexpectedsign.However,theresidualplotsof themodelshowsomespatialpatterns,whicharetobecontrolled for.TheresidualbubbleplotinFig.2showsthespatialdistribution ofthefullmodel’sresidualwithoutspatialcorrelationstructure showsclusteringofpositiveandnegativeresidualsacrossEurope.

(7)

Table2

Nationalparkvisitormodel.Dependentvariableisthelogofannualnumberofvisitorsperhectare.Spatialpatternsintheresidualsarenotcontrolledfor.

Variable Coefﬁcient p-value

(Intercept) 15.64 0.79

Sqrt(grassland) −0.75 0.19

Sqrt(wetland) −1.05 3.49E-02 *

Sqrt(water) 1.32 4.50E-03 **

Log(broadleaf) −0.51 3.70E-03 **

Conifer −0.04 0.18

Log(forestedge) −0.48 0.13

Sqrt(landcoverdiversity) 1.47 3.60E-03 **

Log(viewshed) 0.34 3.28E-02 *

Log(redlistspecies) −0.39 3.71E-02 *

Days>5◦ _6.70E-03 _9.40E-03 _**

NPage 8.08E-03 3.44E-02 *

Log(trails) 0.47 0.00E+00 ***

Log(roads) 0.38 5.00E-03 **

Studyareakm2 _−4.91E-04 _7.00E-03 _**

Log(NPsubstitutes) −0.25 1.13E-02 *

Log(Population50km) 0.48 0.00E+00 ***

GDP/capita −3.50E-05 3.02E-02 *

Surveyyear −1.12E-02 0.70

Surveyquality −2.93E-02 0.68

MultipleR2_:_0.68 _RMSE:_1.21 _AIC:_796.9

AdjustedR2_:_0.65 _CV(RMSD):_0.45 _BIC:_864.5

Signiﬁcantcodes:“***”≤0.001,“**”≤0.01,“*”≤0.05,“.”≤0.1.

Weappliedanumberofdifferenttechniquestocontrolforthese

patterns.

Firstweaddeddifferentregionalfactorvariablestothemodel,

inordertoexplainthespatialpatterns.Wetriedbio-geographical

regions,geographical regionsandcountries4 _as_factor _variables.

However,adding oneof thesevariablesreduced thedegrees of freedomand increasedthecomplexityofthemodel tosuchan extentthatweendedupwithmodelshavingalotofnon-signiﬁcant variables.Alsomostofthedifferentlevelsoftheregionalfactor vari-ablesdidnotshowanysigniﬁcanteffect.Inaddition,AICandBIC valuesdidnotshowanyfavourablescoresforthemodels.

Then, we tried to implement a mixed model byadding the regionalvariablesasarandompartinordertocontrolforthespatial patternsintheresiduals.Wetriedvariouscombinationsofrandom interceptandrandomslopemodels,whichsignificantlyimproved themodelintermsofAICandBICvalues,butaconsiderablespatial residualpatternstillremained.Finally,wetrieddifferentspatial autocorrelationstructures,whichimprovedthemodel’sAICand BICvaluessubstantially,beyondallthemodelswetriedbefore.The bestmodelintermsofAICandBICvaluesaswellasincontrolling forthespatialresidualpatternsappliedisasphericalspatial corre-lationstructure.Theresultofthefullmodelincludingthespatial autocorrelationstructureisshownTable3.Intotal,13predictors ofthefullmodelshowasignificantcorrelationwithtotalannual numbersofvisitsperha.Afterstepwiseeliminationoftheleast significantvariableuntilonlysignificantpredictorsremained(at leastatthe0.1level),weendedupwiththesame13significant predictorsasbeforeandsubstantiallylowAICandBICvalues(see Table4).

Ourﬁnalmodelsshowaspatialautocorrelationbetween sin-gleobservationsuptoarangeof530kmforthefullmodelandup toarangeof580kmfortheﬁnalmodel.Thenuggetrefersto dif-ferencesbetweenobservations,whichcanneitherbeexplainedby themodelnorbythespatialautocorrelationduetomeasurement errorsormicrovariability.

4_For_the_country_variable_we_combined_some_countries_to_one_region_in_order_to

reducethelevelsofthefactorvariable,suchasBeneluxcountries,Alpinecountries andBalticcountries.

Astrongpositiveandhighlysignificantinfluenceisshownfor the presence of water bodies, both in thefull and in the final model. The betacoefficients indicate that it is thefourth most importantpredictorforexplainingrecreationaluseinourmodels. Interestingly,eventhoughwedidnothavestrongprior expecta-tionsregardingthesignsofpredictorsrepresentingthetypenatural vegetation,allofthem–broadleafand coniferousforest, grass-landandwetlands–shownegativesignsinthefullmodel.Only broadleafforestandwetlandsshowasignificanteffectinthefull modelaswellasinthefinalmodel.Alsothevariableforestedges, contrarytoourexpectations,showsanegativeandsignificantsign. However,forestedgesarestronglycorrelatedwithtotalforest(the sumofbroadleafandconiferousforest).Therefore,forestedgesmay pickupsomeofthenegativeimpactsofforestcoveronrecreational useinourmodel.Both,broadleafandconiferousforestshave neg-ativesigns,evenifonlybroadleafforestshowsasignificanteffect. Weinitiallythoughtthatwecouldseparatetheeffectofforestson thenumbersofvisitsfromtheeffectofforestedgesbyincluding singlepredictorsforbroadleafandconiferforest.Oneexplanation ofthenegativesignsofthevegetationcoverpredictorscouldbethat NPsdohavenaturalvegetationtosuchanextent,thatitbecomes abundantandthereby,moreofitdetersvisitors.The transforma-tionofthepredictorsindicatesthattheirnegativeeffectonthe numberofvisitsdecreaseswiththeirincreasingshareoflandcover. Nevertheless,thebetacoefficientsofthesinglevegetation-cover predictorsindicatethattheyonlyhavearelativesmalleffecton thetotalvisitationrate.Alsothepredictormeasuringlandcover diversityshowsasignificantpositiveeffect.Onthecontrary,the predictorviewshedandredlistspeciesabundancedonotproveto haveasignificanteffect.Redlistspeciesabundancehasanegative sign,whichiscontrarytoourexpectations.Nevertheless,both vari-ablesdropoutofthemodelduringthevariableselectionprocedure. Wealsofindapositiveeffectofthenumbersofdayswitha max-imumtemperatureabovefivedegrees.Anotherpredictor,which showsasignificantpositivebutrelativelysmalleffectonthe num-berofvisitsistheageofthenationalpark.Themostimportantand highlysignificantpredictoristheavailabilityoftrails.Inthefinal model,itexplainsalmost17%ofthenumberofvisits.However, thequestionofcorrelationandcausalityisinparticularrelevant forthispredictor.Towhatextenttrailsattractvisitorsandtowhat extenttrailsareputinplaceduetohighvisitornumberscannot

(8)

Table3

Fullmodelincludingsphericalspatialcorrelationstructure.

Variable Coefﬁcient p-value Betacoefﬁcient

(Intercept) −9.30 0.88 2.02%

Sqrt(water) 1.61 3.00E-04 *** 7.26%

Sqrt(grassland) −0.61 0.27 2.53%

Sqrt(wetland) −0.98 3.77E-02 * 3.98%

Log(broadleaf) −0.39 2.72E-02 * 5.74%

Conifer −0.03 0.37 2.31%

Log(forestedge) −0.55 6.91E-02 . 3.84%

Sqrt(landcoverdiversity) 1.34 5.50E-03 ** 5.02%

Log(viewshed) 0.13 0.40 1.54%

Log(redlistspecies) −0.24 0.28 3.46%

Days>5◦ _7.43E-03 _3.59E-02 _* _7.40%

NPage 1.07E-02 5.10E-03 ** 4.98%

Studyareakm2 _−5.69E-04 _3.40E-03 _** _6.12%

Log(trails) 0.44 0.00E+00 *** 14.24%

Log(roads) 0.50 1.70E-03 ** 7.16%

Log(NPsubstitutes) −0.30 1.57E-02 * 7.82%

Log(populationwith50km) 0.37 6.00E-04 *** 11.19%

GDP/capita −1.00E-06 0.96 0.13%

Surveyyear 2.40E-03 0.94 0.17%

Surveyquality −0.12 0.10 . 3.11%

Sphericalspatialcorrelationstructure RMSE:1.26 AIC:768.7

Range:530km,nugget:0.40 CV(RMSD):0.48 BIC:842.7

Table4

Finalmodelafterstepwisemodelselectionincludingsphericalspatialcorrelationstructure.

Variable Coefﬁcient p-value Betacoefﬁcient

(Intercept) −3.35 0.11 2.26%

Sqrt(water) 1.8 0.00E+00 *** 9.29%

Sqrt(wetland) −0.83 4.81E-02 * 3.84%

Log(broadleaf) −0.31 3.41E-02 * 5.18%

Log(forestedge) −0.57 3.32E-02 * 4.53%

Sqrt(landcoverdiversity) 1.32 4.70E-03 ** 5.65%

Days>5◦ _6.89E-03 _3.72E-02 _* _7.83%

NPage 1.07E-02 4.30E-03 ** 5.72%

Studyareakm2 _−5.14E-04 _5.20E-03 _** _6.31%

Log(trails) 0.46 0.00E+00 *** 16.95%

Log(roads) 0.44 2.90E-03 ** 7.26%

Log(npsubstitutes) −0.36 2.50E-03 ** 10.81%

Log(Populationwith50km) 0.32 1.20E-03 ** 10.98%

Surveyquality −0.11 0.1 . 3.38%

Sphericalspatialcorrelationstructure RMSE:1.29 AIC:727.5

Range:580km,nugget:0.38 CV(RMSD):0.48 BIC:782.8

beansweredbythisanalysis.Thesamemayapplytothe

avail-abilityofminorroads,whichalsoshowasigniﬁcantpositiveeffect

butbeinglessimportantforexplainingtheobservedvisitor

num-bers.Asigniﬁcantnegativeimpactcanbefoundforthesizeofthe

studyareaofthevisitormonitoringstudy,butalowbetacoefﬁcient

indicatesarelativelylowimportance.Astrongerandsigniﬁcant,

butnegativeimpactshowstheavailabilityofothernationalpark

areaswithintheregion.Itisthethirdmostimportantvariablein

ourmodels.Thesecondmostimportantvariableinexplainingthe

observednumberofvisitsisthepopulationlivingintheregionof

thestudyarea,whichshowsasigniﬁcantpositiveeffect.Aminor

negativebutnotsigniﬁcanteffectisfoundfortheGDP/capitaand

theyearofthevisitormonitoringsurvey.Thisiscontrarytoour

initialexpectations.Itcouldbethatotherculturalaspectsinterfere

withthiseffect.ItmayalsobethatSouthernEuropeancountries

withlowerGDP/capita(e.g.ItalyandSpain)receivemorevisitors

inNPsbecauseofhightouristvisits,whereasricherNorthern

Euro-peancountries(e.g.Scandinaviancountries)receivefewervisitors

becauseoflowertouristnumbers.Attheedgeofthe0.1signiﬁcance

level,thepredictormeasuringthequalityofthevisitor

monitor-ingstudyshows arelatively smalland negativeeffect.Initially,

thisvariablewasconsideredforexplainingresidualpatterns.We

expectedthatvisitormonitoringstudieswithalowerquality

judg-mentwouldresultinlessprecisevisitorestimatesandthereforein

higherresiduals.However,inourpre-analysiswecouldnotﬁnda

signiﬁcanteffectofthevisitormonitoringqualityontheresiduals.

Moreover,weﬁndthatvisitormonitoringstudiesoflowerquality

tendtooverestimatevisitornumbers.Thiscouldbecausedbythe

incentiveofNPmanagerstohighlighttheimportanceoftheirNP

andtherebyuseassumptionsmadewithinthevisitormonitoring

studyinfavourofhighervisitornumbers.Visitormonitoring

stud-iesofhigherqualitymayallowforlessoftheseassumptionstobe

made(bymorecompletecountingandlessup-scaling).

Further-more,completereportingoftheassumptionsmademaystimulate

morerealisticjudgments.

WeusedourﬁnalmodeltomakepredictionsforallNPssitesin

ourprimaryvisitordatabaseandalsoforallNPsinmostoftheEU5

aswellasinNorwayandSwitzerland.Comparingourpredictions

5_We_could_not_make_predictions_for_some_EU_countries_for_which_we_are_missing

rasterlayersoftheexplanatoryvariablesinourmodel.ThesecountriesareBulgaria, Croatia,Cypress,IslandandMalta.

(9)

withourprimarydata,weestimateanaveragerelativeprediction errorofabout185%(thefullmodel174%),whichseemsreasonably good.Interestingly,thefourobservationscontributingmosttoour relativepredictionerrorarealllocatedinItaly.

Usingourmodeltopredictthenumberofvisitstoall449NPs acrossourstudyarea,weestimateatotalannualnumberofvisits ofmorethan2billion(2,016,028,000;lowerandupper95% con-ﬁdenceinterval:1,217,818,000;3,404,254,000).6_Combining_this estimatewiththeaveragemonetaryvaluepervisit(7.17_D,prices 2013),whichweextractedfromameta-analysisofrecreation valu-ationstudies,thetotalrecreationalvalueofthe449NPsamountsto D14.5billionannually.Theresultcompareswelltotheestimates ofBalmfordetal.(2015),whoestimate3.8billionvisitsannually andavalueof$US26.9billionforallprotectedareaswithinEurope, notonlyNPs.Ouraggregatedestimatespercountryareshownin Table5.

Most visits are received by British NPs, which resultsfrom thelargetotal areaof NPs,highpopulationdensity and inten-siverecreationalfacilities interms of traildensities.Alsoother denselypopulatedcountriessuchasDenmark,Belgiumandthe Netherlandsshowrelativelyhighvisitornumbers.Onthecontrary, countriessuchasSweden, FinlandandNorway showrelatively lowvisitornumbersfortheirlargeandmainlyforestedNPsinthe lowpopulatednorth.Germanyshowsexceptionallyhighvisitor numbersconsideringtherelativelysmallNParea.However,these numbersaredominatedbyonelargeNP,forwhichourmodelmay overpredictthetotalnumberofvisits.TheWaddenSeaNP–an UNESCOnaturalheritage–isbyfarthelargestNPofGermanyand stretchesalmostallalongtheNorthSeashoreofGermany.Thearea liesinthecatchmentoflargecitiessuchasHamburgandBremen. Itisatouristichotspotreceivingbyfarthehighestnumberofday andovernightvisitsofGermanNPs(Job,Woltering,Harrer,2010). Allvariablesusedinourmodel,exceptsize,showvaluesinfavour ofhighvisitornumbersfortheWaddenSeaNP.Thiscombination ofsuchvariablevaluesisexceptionalinourdataandmaycausean unreasonableoverprediction.

OurpredictionsofvisitsperhaforamarginalincreaseofNP supplyinEuropeareshowninFig.3.Weassumeahypothetical newlycreatedNPof about80km anywherethroughout Europe anestimatethenumberofvisitsitwouldreceive.Allurbanareas areexcludedfromthisprediction(EEA,2015b),asitseems unre-alisticthatsuchareaswouldbeconvertedintoaNPandbecause urbanareasaretypicallycharacterizedbyexplanatoryvariable val-uesthatliebeyondtherangeoftheexplanatoryvariablevaluesof ourprimarydata.Themapshowsvaluesfromalmostzeroupto themaximumofabout147,000annualvisitsperha.Lownumbers ofvisitsarepredictedforremoteareas,whicharecharacterized bylowpopulationandlittleaccessinfrastructure.Themaximum predictedvisitsof147,000perhaseemshigh,but34visitorsfor anaveragedaylighthourmaynotbeunreasonableforapopular visitorhotspotinaNP.However,itshouldbeconsideredthatthe predictedvisitornumbersarestronglyskewedwithameanand medianvaluesofabout87and4.8.Morethan90%ofthepixels receivevisitorsoflessthan100visitorsayearandanythingabove 2000istobeexpressedinpermile.Amappresentingthespatially expliciteconomicvaluescanbefoundinAppendixof Supplemen-tarymaterial(Fig.S1andS2).

Toexemplifyourmodelforarealisticsettingwechosethe Teu-toburgerforestandtheSenneheathlandinthewestofGermany, whichisproposedforNPdesignation.Fig.4showshowthe pre-dictedannualvisitsperhadistributedacrossthearea.Onaverage, weexpectabout283annualvisitsperhaforthearea.Thehighest

6_We_used_the_rms_R-package_for_estimating_conﬁdence_intervals.

Fig.3.PredictedvisitsperhaandyearforapotentialnewNationalParkofabout 80km2_.

Notethatthepredictedtotalvisitornumberoftheentireareaislessthanthesum ofthepredictedvisitorsforeachhabecauseoftworeasons:visitorsmaycrossmore thanonehaduringavisitanditisnotpossibletotakethelinearmeanofamodel containingnon-linearvariables.

Fig.4.PredictedvisitsperhaandyearforapotentialNationalParkinthe Teuto-burgerforestandtheSenneheathland(westofGermany).

visitationrateispredictedintheperipheralareas,closetothe pop-ulationcentresofcitiesofDetmoldandPaderbornreceivingupto 24,000visitsperhaandyear.Incontrast,thecenteroftheproposed NP,whichishardlyaccessible,ispredictedtoreceivelessthanone visit/ha/year.In totalwe predictabout5.8millionannualvisits

(10)

Table5

EstimatesoftotalannualvisitstoNationalParksinEuropeancountriesandtheirestimatedmonetaryvalue.

Country Km2_of_NP _Predicted_Visits _95%_Conﬁdence_Interval_{(lower/upper)} _Monetary_Value

Austria 3098 24,098,000 14,001,000/41,660,000 172,684,000D Belgium 3200 63,569,000 32,294,000/125,388,000 455,527,000D Switzerland 170 135,000 72,000/256,000 969,000D CzechRepublic 3543 32,835,000 17,148,000/63,127,000 235,290,000D Germany 2363 534,188,000 309,773,000/921,987,000 3,827,911,000D Denmark 846 77,623,000 55,797,000/108,203,000 556,236,000D Spain 10,450 121,666,000 89,810,000/170,467,000 871,840,000D Estonia 1618 2,182,000 1,561,000/3,078,000 15,635,000D Finland 8196 6,427,000 4,564,000/9,456,000 46,054,000D France 13,565 71,408,000 36,506,000/140,680,000 511,700,000D UnitedKingdom 21,754 700,862,000 429,126,000/1,162,686,000 5,022,270,000D Greece 4677 14,713,000 10,287,000/21,934,000 105,432,000D Hungary 6234 18,543,000 11,457,000/30,336,000 132,878,000D Ireland 2221 3,510,000 2,447,000/5,070,000 25,152,000D Italy 17,419 145,719,000 93,198,000/231,777,000 1,044,203,000D Lithuania 1345 2,398,000 1,482,000/3,909,000 17,186,000D Luxembourg 465 2,912,000 1,560,000/5,441,000 20,866,000D Latvia 3201 3,711,000 2,508,000/5,538,000 26,592,000D Netherlands 1889 93,133,000 48,749,000/182,005,000 667,375,000D Norway 30,696 2,150,000 1,821,000/2,602,000 15,404,000D Poland 10,168 46,227,000 25,125,000/85,506,000 331,254,000D Portugal 930 15,245,000 10,006,000/23,227,000 109,244,000D Romania 5670 2,662,000 1,565,000/4,546,000 19,077,000D Slovakia 7679 18,218,000 9,079,000/37,180,000 130,544,000D Slovenia 1157 4,121,000 2,425,000/7,004,000 29,531,000D Sweden 8370 7,773,000 5,457,000/11,191,000 55,700,000D

fortheentirearea(95%conﬁdenceintervallowerbound3.38and

upper9.91million),7_which_accounts_for_an_annual_monetary_value

ofapproximatelyD 41.5million.Amappresentingthespatially expliciteconomicvaluescanbefoundintheAppendixof Supple-mentarymaterial(Fig.S1andS2).8_Negative_impacts_on_the_number ofvisitsincludetherelativelylowpresenceofwaterbodies,high forestcover,lowtrailavailabilityandthelowageofthepotential newNP.PositiveimpactsincludethesmallsizeoftheNP,thehigh populationpressure,lowsubstituteavailabilityandthehighland coverdiversity.Thenumberofvisitsisexpectedtoincreasewith theageoftheNPandifrecreationalfacilitiesareestablished. 5. Discussion

5.1. Spatialeffectsandmodelling

Ourestimatedmodelﬁtsthedatareasonablywellandtherefore offersvaluableinformationonthemaindriversofrecreationaluse withinEuropeanNPs.Allpredictorswithstatisticallysigniﬁcant effectsonthenumberofrecreationalvisitshavesignsthatarein linewithourinterpretationsandtheoreticalexpectations.

Nevertheless,therearesomeuncertaintiesin themodel and predictionaccuracywhichmaybeimprovedbyfurtherresearch. Thequestionremains,whatmaybethesourceofthespatial auto-correlation.Inanoptimalstatisticaltextbookworld,introducing spatialautocorrelationinamodelwouldnotinﬂuenceparameter estimates,butonlyreducethedegreesoffreedomofthemodel. However,lookingat realworldspatial data,this is hardlyever thecase.Ifparameterestimatesareaffectedasinourcase, this mayindicatesomecommonspatialeconometricproblems,suchas missingpredictors,whicharepickedupbythespatialerrorterm,a spatialweightmatrixoranon-linearrelationship(Diggle,Morris,&

7 _Note_that_the_map_displaying_the_recreational_ecosystem_service_values_should

beinterpretedwithcaution,becausewedonotaccountforspatialvariationswithin thevalueperrecreationalvisit,whichmayalterthetotalvalueestimateforcertain locationsconsiderably.

8 _Note_that_for_illustrative_purpose_the_color_scheme_is_set_to_display_the_same

amountofpixelspercolorshade.

Wakefield,2000;Smith&Lee,2011;Fingleton&LeGallo,2010).A likelyexplanationcouldbethatunobserveddeterminantsof recre-ationalvisitsexist,whicharespatiallyrelated.Suchdeterminants couldbemanifoldandincludeeverythingfromsite,contextand methodologicalstudycharacteristicsaswellastheirinteractions. Oneimportantaspectcouldberelatedtothesocial-culturalcontext andpathdependencies,whichmayresultinspecificrecreational patternsincertaincountriesandregions.Alsodifferingproperty rightscouldplay animportantrole. Investigatinghuman recre-ationalbehaviouracrossastudyareaasbigasEuropeissucha complexissuethatalloftheseeconometricproblemsmayarise. Theremayhardlybeanymodelthatcanincorporateallrelevant driversofrecreationaluse,theirinteractionsandnon-lineareffects. Encounteringsuchproblemsiscommonformodellingspatial dataandtherefore,wehavetobecautiousininterpretingp-values andparameterestimates.Anoption togainfurtherinsightsand confidenceinmodelresultinterpretationsistotrydifferent spa-tialmodellingapproachesandcomparetheirresults.Inparticular, comparetheconfidenceintervalsoftheparameterestimates.There isnumberofmodelsetups,whichwouldqualifyforevaluatingsuch spatialdatasets.Sinceinthisstudyweareanalysingcountdata, oneoptionwouldbetouseanegativebinomialoraquasi-Poisson distribution,eventhoughitshouldnotchangethemodelresultstoo much(O’HaraandKotze,2010).However,thereareonlyavery lim-itednumberofstatisticalR-packages,whichallowforcombining thesedistributionswithspatialautocorrelationsandaswestated abovewehadproblemsinsolvingthemaximizationalgorithms forthesemodels.Anoption wouldbetotryalternativemodels incorporatingthespatialautocorrelationeitherwithinthefixed ortherandompartofthemodel,suchasaspatiallagmodel,a Durbinmodel,spatialautoregressivemodelswithautoregressive disturbances,geographicallyweightedregressionsorevenbyusing Bayesianapproaches. However,thereis noconsensusonwhich modeltousebestforthisspecificpurpose.Fittingalloratleast someofthesemodelsandcomparingtheirresultsmaybesubject tofurtherresearch(Bivand,2011;Elhorst,2010;Gerkman,2011; Brunsdon,StewartFotheringham,&Charlton,1996).

Nevertheless, consideringthe complexity ofthe spatial pro-cessesdrivinghumanrecreationalbehaviour,wecanconﬁdently

(11)

saythatwemodelrecreationalusereasonablywell.Noneofthe predictors’signsdifferacrossthedifferentestimatedmodels, nei-therformodelswithoutautocorrelationnorforthemixedmodels, whichindicatestherobustnessofouranalysis.Anyhow,other pub-licationsconductingspatialmodellingofrecreationalusedonotat allengagetosuchadepthinthespatialdimensionsnordotheytake intoaccountsuchconsiderationsonuncertainties,potential alter-nativeregressiontechniquesandmodelsetups(Neuvonenetal., 2010;Mills&Westover,1987;Hanink&Stutts2002;Hanink & White1999).

Futureresearchonthisissuemaybenefitfromgreaterandmore reliableprimarydata availability.Errorsin primarydata collec-tionimposehugedifficultiesfor identifyingrelevant predictors. In recent years visitor monitoringstudies encountered a huge dynamicintermsofinterestandtechnicaladvancement.Recent remotecontrolledelectronicvisitorcountersallowfarmore accu-ratevisitorestimatesatlowercostsascomparedtoconventional personalcounting.MorerefinedGISdatasetsmayallowformore accurate, detailed and comprehensive predictors for modelling recreationaldemand.

5.2. Valuationofrecreationalservices

Anotheraspectofimprovementmaybetoaccountforspatial variationsinthevalueperrecreationalvisitbyapplyingavalue functiontransfer(suchasmeta-analyticvaluetransfer).Using a unitvaluetransferformappingecosystemservicevalues across alargerareaisassociatedwithtransfererrors,inparticularwith so-calledgeneralisationerrors.Nevertheless,thevalueofa recre-ationalvisit variesacrossspaceduetodifferencesinecosystem characteristicsandthelocalpopulation’spreferences,differences thatarenotaccountedforinaunitvaluetransfer(Rosenberger &Loomis,2001).Valuefunctiontransferallowsadjusting trans-ferredvaluestositespeciﬁccircumstancesandmaythereforebe moreaccurateforecosystemservicevaluemapping.However,even though,value function transfer is consideredtoproduce lower transfer errors, there is no consensus onwhich value transfer methodis bestfor speciﬁc circumstances.Evidence ontransfer errorsshowmixedresultsandunitvaluetransfermaybesuperior toothervaluetransfertechniquesforsomeapplications(Navrud& Ready,2007;Rosenberger&Phipps,2007;Johnston&Rosenberger, 2010;Brouwer2000; Rosenberger&Stanley,2006;Lindhjem & Navrud,2008).Inecosystemservicevaluemapping,theunitvalue transfermethodismostcommon(Schägneretal.,2013).Itisalso proposedfortheaggregationofvaluestosetup,forexample,a Sys-temofEnvironmental-EconomicAccounting(SEEA),eventhough aggregationoveralargeareaiscontroversialandshouldbe inter-pretedwithcaution(UN,2014;Costanzaetal.,1998).

Inthecaseofrecreationalservices,meta-analysisofrecreational valuationstudiesshowthatmostof thevariationsinthevalue pervisit result fromdifferentvaluation methodologiesand not from site specific circumstances, indicating large measurement errors.Moreover,it remainsdifficulttoidentifyrobust relation-shipsbetweenspatial explanatoryvariablesand thefinal value estimate.Meta-analysisonrecreationalvaluationstudiesidentify onlyfewsignificantandtypicallyweakeffectsofbiophysical, socio-economic and regional or national dummyvariables (Shrestha, Rosenberger, & Loomis,2007; Zandersen &Tol, 2009; Brander, Eppink,Schägner,vanBeukering,&Wagtendonk,2015;Senetal., 2013; Sen et al., 2011; Rosenberger & Loomis, 2001; Londo ño &Johnston,2012).By usingthemean valueof a largenumber ofprimaryvaluationstudies,weaimataveragingout measure-menterrors withinourvaluetransfer (Johnston,Elena,Ranson, 2006),whichmayresultinlowertransfererrorsascomparedto theusageofsinglestudiesorregionalsubsets,eventhough cul-turaldifferencesacrosscountriesmayaffectvalueperrecreational

visit(Ready&Navrud,2006;Shrestha&Loomis,2001;Lindhjem &Navrud,2008;Kaul,Boyle,Kuminoff,Parmeter,&Pope,2013; Hynes,Norton,&Hanley,2012).

Finally,theoverallrecreationalvalueofasiteispredominantly determinedbyspatialvariationsinthenumberofrecreational vis-its.Spatial variationsin valueper recreationalvisit playonly a minorrole(Batemanetal.,2006;Jonesetal.,2003).Thisinsight isalsosupportedbymeanrelativedeviationsofourprimarydata, which is considerably higher for the visitor numbers as com-paredthevaluepervisitestimates.Inconsequence,accuratevisitor estimatesarebyfarmoreimportantfordeﬁningtheoverall recre-ationalvalueofacertainlocationthanaccurateestimatesofthe recreationalvaluepervisit.Ascomparedtometa-analysisof recre-ational valuationstudies, the explanatorypower of ourspatial variablesexplainingvisitornumbersishigh.Wearetherefore con-ﬁdentthatwecapturethemainspatialvariationsintheoverall recreationalvalueNPrecreationandthatthevalueestimatesgive agoodindicationoftherelativeimportanceofEuropeanNP recre-ation as compared to otherecosystem services and man-made goods.

5.3. Policyimplications

Themodelcanbeusedforanumberofpolicyapplications:(1) ThemodelmaycontributetothefulfilmentsoftheEUbiodiversity strategy2020,whichrequireofEUmembersstatesto“mapand assessthestateofecosystemsandtheirservicesintheirnational terri-toryby2014,assesstheeconomicvalueofsuchservices,andpromote theintegrationofthesevaluesintoaccountingandreportingsystems atEUandnationallevelby2020”(EC,2011)andtheachievement oftheAichiTargets,whichaimat“reflectingthevaluesof biodiver-sityinspatialplanningandresourcemanagementexercisesincluding throughthemappingofbiodiversityandrelatedecosystemservices” (CBD,2013).(2)Themappedrecreationalvisitornumbersandthe relatedeconomicvalueofrecreationalESScanactasaspatialvalue databasethatcanbeusedforvaluetransfers.Policymakerscan quicklyderiveavalueestimateoftherecreationalservicesofany NPacrossEuropebyconsultingthemap.(3)Themapsmay con-tributetoanefficientresourceallocationbyallowingpolicymakers toprioritizeareasforconservationduetotheirhighrecreational value.Inaddition,recreationalinfrastructuremaybedesignedto matchtheneedsoftheexpectedvisitornumberswithinagivenNP. Furthermore,itmaybevaluabletocomparethemodel’s predic-tionswithrealworldobservationsonrecreationaluseandvalues (ifavailable)and,forexample,investigatewhysomeNPsmight remainbelowtheirrecreationalpotentialandhowtherecreational useanditsvaluecouldbeincreased.However,itshouldbenoted thatthemodelallowsonlyforassessmentsofNP.Evenif predic-tionscanbemadeforanewhypotheticalNP,noconclusioncanbe madeonwhetherNPdesignationresultsinanincreaseordecrease ofrecreationaluseanditsvalues.(4)Themodelallowstoevaluate theeffectoflandusepolicieswithinEuropeanNPonrecreational servicesandvalues.(5)Finally,theestimatedrecreationalservice valuesmaycontributetothesettingupofagreenGDPora Sys-temofEnvironmental-EconomicAccounting(SEEA)asproposed bytheUN(2014),whichmayactasacounterparttotraditional GDPaccountsandrepresentanadditionalmeasurefortheimpacts ofhumanactiononhumanwell-being.

6. Conclusion

We model recreational use of European NPs using a large numberofspatiallyvariablepredictors.Ourmodelﬁtsthedata rea-sonablywellandweidentifythemaindeterminantsofvariationin recreationaluseinEuropeanNPs.Amonganalysedvariablestrails

(12)

density,populationdensity,presenceofsubstitutes,presenceof waterbodiesandnumberofdayswithtemperatureabove5◦are thosethatshowahigherexplanatoryvalue.Themodelallowsthe estimationandvaluationoftotalrecreationaluseofexistingand plannedNPs.ForourstudyareacoveringmostofEuropeandin total449NPs,weestimatemorethan2billionrecreationalvisits ayear,withaneconomicvalueofapproximatelyD 14.5billion. Thelatterinformationisparticularlyrelevanttosupportthetask thatEUcountriesshouldfulﬁlby2020,accordingtoEC(2011)of assessingtheeconomicvalueofecosystemservicesandintegrate suchvaluesintoaccountingandreportingsystemsby2020.

Since allourpredictors areobtained fromGIS raster layers, whichcoverentireEurope,themodelcanbeappliedforex-ante evaluationofalternativepolicyscenariosofchangeforexistingNPs andonthecreationofnewNPsataEuropeanscale.This informa-tionmaybeusefulinplanningthesupplyofrecreationalfacilities suchasparkingandaccommodation.Furthermore,NPlocations anddesignfeaturesoptimizingrecreationalusecanbeidentified. Thereby,themodel hasimplicationsfor NP policyofEuropean countries.Basedonourfindings,wecanconcludethattoensure highnumbersofrecreationalvisits,potentialnewNPsshouldbe locatedincloseproximitytopopulatedareasbutnotclosetoother NPs.Thetotalconservationareashouldbeusedforalargernumber ofsmallparksratherthanforasmallernumberoflargeones.The availabilityofwaterbodiesandthediversityofthelandcover con-tributetohighervisitationrates,whereasextensiveforestcover tendstodetervisitors.However,itshouldbekeptinmindthatthe mainpurposesofNPsarenottosupplyrecreationalservicesbut preserveabeautifulandnaturallandscapeaswellasbiodiversity forposterity.Recreationalopportunitiesareaco-benefitofNPs, whichcanbeusedasanargumentforallocatingresourcestowards NPcreationandconservation.

Acknowledgements

We would like to thank everybody who supplied primary dataon visitorcounts of EuropeanNP. In particular we would liketothankIgnacioPalomofromtheAutonomousUniversityof MadridSpain,LaurenceChabanisfromParcsNationauxdeFrance, HubertJobandManuelWolteringfromtheUniversityofWürzburg Germany,BertDeSomvielefromOrganisationforForestin Flan-dersandJeroenGilissenfromNPHogeKempenBelgium,Martin GoossenfromAlterraNetherlands,Arne Arnbergerand Thomas SchauppenlehnerfromBOKUAustria,DavidBaumgartnerfromNP HoheTauernAustria,MetteRohdefromVisitdjurslandDenmark, BoBredal Immersen from NP Thy Denmark, Camilla Nasstrom fromNaturvardsverket Sweden, KajalaLiisa fromMetsähallitus and Marjo Neuvonenfrom Natural Resources InstituteFinland, KrystynaSkarbekfromtheMinistryoftheEnvironmentPoland, AnnamáriaKopekfromDirektiondesNPBalaton,IrenaMuskare fromNatureConservationAgencyLatviaandeverybodyweforgot tomentionhere.TheresearchwasfundedbytheJointResearch Centre,EuropeanCommission.

AppendixA. Supplementarydata

Supplementarydataassociatedwiththisarticlecanbefound,in theonlineversion,athttp://dx.doi.org/10.1016/j.jnc.2016.03.001. References

Balmford,A.,Green,J.M.H.,Anderson,M.,Beresford,J.,Huang,C.,Naidoo,R.,etal. (2015).Walkonthewildside:estimatingtheglobalmagnitudeofvisitsto protectedareas.PLoSBiology,13(2),e1002074.http://dx.doi.org/10.1371/ journal.pbio.1002074

Barbault,R.(2011).2010:anewbeginningforbiodiversity?ComptesRendus Biologies,BiodiversityinFaceofHumanActivities/LaBiodiversiteFaceAux

ActivitesHumaines,334(5–6),483–488.http://dx.doi.org/10.1016/j.crvi.2011. 02.002

Bateman,I.J.,Abson,D.,Beaumont,N.,Darnell,A.,Fezzi,C.,Hanley,N.,etal.(2011).

Chapter22:economicvaluesfromecosystems.InUKnationalecosystem assessment:understandingnature’svaluetosociety,technicalreport,1466. Cambridge:UNEP-WCMC.

Bateman,I.J.,Day,B.H.,Georgiou,S.,&Lake,I.(2006).Theaggregationof environmentalbenefitvalues:welfaremeasures,distancedecayandtotal WTP.EcologicalEconomics,EnvironmentalBenefitsTransfer:Methods, ApplicationsandNewDirectionsBenefitsTransferS.I.,60(2),450–460.http://dx. doi.org/10.1016/j.ecolecon.2006.04.003

BatesD.,MaechlerM.,BolkerB.,WalkerS.,BojesenChristensenR.H.,SingmannH., DaiB.,GrothendieckG.(2015).lme4:LinearMixed-EffectsModelsUsing “Eigen”andS4(version1.1-8),https://cran.r-project.org/web/packages/lme4/ index.html.

BatistaeSilva,F.,Gallego,J.,&Lavalle,C.(2013).Ahigh-Resolutionpopulationgrid mapforeurope.JournalofMaps,9(1),16–28.

Biavetti,I.,Karetsos,S.,Ceglar,A.,Toreti,A.,&Panagos,P.(2014).European meteorologicaldata:contributiontoresearch,developmentandpolicy support.ProceedingsofSpie:TheInternationalSocietyforOpticalEngineering, 9229,922907.http://dx.doi.org/10.1117/12.2066286

Bivand,R.(2011).After‘Raisingthebar’:appliedmaximumlikelihoodestimationof familiesofmodelsinspatialeconometrics.SSRNscholarlypaperID1972278. SocialScienceResearchNetwork:Rochester,NY.http://papers.ssrn.com/ abstract=1972278

Bivand,R.S.,Pebesma,E.,&Gómez-Rubio,V.(2013).Appliedspatialdataanalysis withR(2nded.).NewYork:Springer.

BolkerB.,SkaugH.,MagnussonA.,PetersenA.H.(2012).GettingStartedwiththe glmmADMBPackage.

Brander,L.M.,Eppink,F.V.,Schägner,J.P.,vanBeukering,P.J.H.,&Wagtendonk, A.(2015).GIS-Basedmappingofecosystemservices:thecaseofcoralreefs.In beneﬁttransferofenvironmentalandresourcevalues.InR.J.Johnston,J.Rolfe, R.S.Rosenberger,&R.Brouwer(Eds.),Theeconomicsofnon-marketgoodsand resources14(pp.465–485).Netherlands:Springer.http://link.springer.com/ chapter/10.1007/978-94-017-9930-020

Brouwer,R.(2000).Environmentalvaluetransferstateoftheartandfuture prospects.EcologicalEconomics,32,137–152.

Brouwer,R.,Langford,I.H.,Bateman,I.J.,&Turner,R.K.(1999).Ameta-analysisof wetlandcontingentvaluationstudies.RegionalEnvironmentalChange,1(1), 47–57.

Brunsdon,C.,StewartFotheringham,A.,&Charlton,M.E.(1996).Geographically weightedregression:amethodforexploringspatialnonstationarity. GeographicalAnalysis,28(4),281–298.http://dx.doi.org/10.1111/j.1538-4632. 1996.tb00936.x

CBD(ConventiononBiologicalDiversity).(2013).Aichibiodiversitytargets.aichi biodiversitytargets.https://www.cbd.int/sp/targets/

Costanza,R.,d’Arge,R.,DeGroot,R.,Farber,S.,Grasso,M.,Hannon,B.,etal.(1998).

Thevalueofecosystemservices:puttingtheissuesinperspective.Ecological Economics,25(1),67–72.

Diggle,P.,Ribeiro,J.,&Peter,J.(2015).geoR:analysisofgeostatisticaldata(version1. 7-5.1).https://cran.r-project.org/web/packages/geoR/index.html

Diggle,P.J.,Morris,S.E.,&Wakeﬁeld,J.C.(2000).Point-sourcemodellingusing matchedcase-controldata.Biostatistics,1(1),89–105.http://dx.doi.org/10. 1093/biostatistics/1.1.89

Drakou,E.G.,Crossman,N.D.,Willemen,L.,Burkhard,B.,Palomo,I.,Maes,J.,etal. (2015).Avisualizationanddata-sharingtoolforecosystemservicemaps: lessonslearnt,challengesandthewayforward.EcosystemServices,http://dx. doi.org/10.1016/j.ecoser.2014.12.002

Dramstad,W.E.,&Tveit,M.S.(2006).Relationshipsbetweenvisuallandscape preferencesandmap-Basedindicatorsoflandscapestructure.Landscapeand UrbanPlanning,78(4),465–474.http://dx.doi.org/10.1016/j.landurbplan.2005. 12.006

EC(EuropeanCommission)(2006).ForestMapping.JRC(JointResearchCenter) ForestCoverMapshttp://forest.jrc.ec.europa.eu/download/data/. EC(EuropeanCommission)(2011).TheEuBiodiversityStrategyto2020,

Luxembourg.

EC(EuropeanCommission)eurostat.(2013).Eurostat:yourkeytoeuropean statistics.http://ec.europa.eu/eurostat/home

EEA(EuropeanEnvironmentAgency).(2006).CORINElandcover.http://www.eea. europa.eu/publications/COR0-landcover

CDDA(CommonDatabaseonDesignatedAreas).(2013).CDDA(Commondatabase ondesignatedareas).http://www.eea.europa.eu/data-and-maps/data/ nationally-designated-areas-national-cdda-4

DigitalElevationModeloverEurope(EU-DEM).(2015).Digitalelevationmodelover europe(EU-DEM).http://www.eea.europa.eu/data-and-maps/data/

eu-dem#tab-european-data

EEA(EuropeanEnvironmentAgency)(2015).Urbanmorphologicalzones2000. Data.http://www.eea.europa.eu/data-and-maps/data/

urban-morphological-zones-2000-2.

EG(eurogeographics)(2010).EuroRegionalMap.http://www.eurogeographics.org/ products-and-services/euroregionalmap.

Ejstrud,B.(2006).Visitornumbersandfeasibilitystudies.predictingvisitor numberstodanishopen-airmuseumsusingGISandmultivariatestatistics. ScandinavianJournalofHospitalityandTourism,6(4),327–335.http://dx.doi. org/10.1080/15022250600929270