Unravelling the relationship between vegetation change and geodiversity in the Dutch dunes: a case study

(1)

Unravelling the relationship between

vegetation change and geodiversity in the

Dutch dunes: a case study

Figure 1. Kennemerland Zuid dunes (van den Berg, 2021).

Thijmen van Hessen 12372560 03–06-2021 Amsterdam

Harry Seijmonsbergen & Jim Groot BSc Future Planet Studies

(2)

Abstract

Climate change and human expansion are changing coastal dune environments. To manage and preserve these vulnerable regions conservationists increasingly rely on remote sensing techniques for their work. The link between biotic and abiotic features is hereby an important relationship to take into account. This study aims to assert whether there is a relationship between geodiversity and vegetation cover in the Dutch dune area Kennemerland Zuid. To do this the geodiversity index (GDI) is compared to the change of the Normalized Difference Vegetation Index (NDVI) between 2010 and 2018. The workflow of this study is divided into three parts and all automated in ArcGIS Pro. The results of this study showed no significant relationship between the GDI and change in NDVI. This unexpected result and the discussion still gave information on the spatial temporal complexity of the dunes. Future research can use this paper for further linkage between conservation and geodiversity, and expand the method for index based geodiversity assessments.

(3)

1. Introduction:

Many natural environments are likely to be affected by climate change effects and sea levels will continue to rise in the coming decades due to the melting of polar ice caps and thermal expansion of ocean water (Maul, 1993; Nicholls & Cazenave, 2010). Of all natural environments in the world coastal dune systems are among the primary regions that are and will be threatened (Martinez et al., 2006). To protect the vulnerable coastal regions the Netherlands has one of the most extensively managed dune systems in the world after in 1990 the Dutch government adapted a “Dynamic Preservation” strategy that aims to maintain the position of the shoreline through tactics such as artificial sand nourishment (Arens et al., 2013).

Coastal dune environments are unique with a high diversity in terms of morphology, vegetation and dynamics (Anfuso & Del Pozo, 2009). The elevation and relief of dunes are extremely variable compared to other natural landforms (in the Netherlands) due to the many natural processes such as aeolian sand deposition and coastal erosion (Łabuz, 2016) which are reinforced by the relative lack of vegetation cover (Delgado-Fernandez et al., 2019). Dune dynamics are influenced by vegetation cover and vegetation cover influences dune dynamics (Stallins & Parker, 2003). This relationship is very complex and a new perspective on this linkage is provided by this study by looking at the relationship between ‘geodiversity’ and vegetation cover.

Conservationists have increasingly recognized a new way of preserving ecosystems and landscapes namely through ‘geodiversity’ (Erikstad 2008; Koster, 2009; Zarnetske et al., 2019). Geodiversity defined by Gray (2004) is “the natural range (diversity) of geological (rocks, minerals, fossils), geomorphological (land form, processes) and soil features”. Through the use of the

geodiversity, conservationists can better preserve and manage the ecosystems since these are highly connected and influenced by each other (Brilha et al., 2018; Chakraborty & Gray, 2020).

1.1 Goal and research questions

For the sustainable management and planning of coastal dune environments it is of increased importance that the link between geodiversity and vegetation cover is understood. In order to better understand this relationship this study aims to map the geodiversity and vegetation cover on a local scale in the designated study area in the Natura 2000-area Kennemerland-Zuid, and analyzing the relationship between geodiversity and vegetation cover. This resulted in the main research question.

“

How is the geodiversity in the Natura 2000-area Kennemerland-Zuid related to the

changes in vegetation cover between 2010 and 2018?”

To answer this question two sub-questions were formulated. 1. What is the geodiversity of the study area?

2. Where do the changes in vegetation cover occur within the study area?

To answer these questions, geodiversity features and changes in vegetation cover will be calculated and mapped out in a rasterized grid to analyze. This research paper is structured in the following way; first geodiversity will be defined, secondly the research method will be described using a workflow after which the results will be presented and thirdly these results will be discussed in relation to relevant literature to answer the research questions. This paper will then end with a

(6)

1.2 The concept of geodiversity

Even though conservationists, ecologists and policymakers have long recognized geodiversity as a key part in the development of nature-based solutions to environmental challenges, natural

resources demand and ecosystem preservation (Lawler et al., 2015), in practice the use of geodiversity for these solutions leaves much to be desired and is mostly applicated indirectly (Comer et al., 2015). As geodiversity is a relatively new concept the different approaches and ways to quantify geodiversity make room for broad purposes and applications but also creates discrepancies as it is still going through a process of conceptual self-affirmation (Araujo & Pereira, 2017). There are some

geodiversity case studies which have successfully been incorporated into conservation plans presented in Anderson et al. (2015) and Hjort et al. (2015) and all of which found some sort of correlation with the distribution of vegetation. These studies all used similar primary elements to construct geodiversity (e.g. topography and soils) which are also used in Seijmonsbergen et al. (2018) that forms the basis of this research.

(7)

2. Methods

2.1 Research area

The chosen study area is the Natura 2000-area Kennemerland-Zuid in the provinces of Noord Holland and Zuid Holland in the Netherlands shown in Figure 2. This research area was chosen because it is a dune ecosystem that has had some big changes in management and is one of the biggest consecutive dune areas in the country. Since 2009 conservationists have set up several projects to improve dune mobility such as making notches in the foredune, removing organic top soil layers in dune slacks and turfing. The hope is that succession of vegetation is slowed down which would result in an increase of white dunes, wet dune valleys and dune migration (Provincie Noord-Holland, 2018).

Figure 2. Map of the Natura 2000-area Kennemerland-Zuid located in the Netherlands.

2.2 Workflow

The methodology used in this study is divided into three parts and a workflow was designed to streamline process of this research. The workflow of the method is presented in Figure 3 and consists of three main parts; 1) pre-processing, 2) analysis and 3) deliverables and results, all done in ArcGIS Pro (Version 2.4; Esri Inc, 2020). In the Digital Appendix (Appendix D) the ArcGIS Pro Project can be found with the Geodatabase and the automated models used for the workflow (also found in Appendix C).

The pre-processing phase consists of data acquisition, transforming the data to the same coordinate system, creating a grid and making sure the layers cover the whole study area and rasterizing the maps. The analysis phase consists of four steps; 1) calculating the geodiversity index (GDI) using three sub-indices, 2) calculating the vegetation cover in 2010 and 2018 and the change in vegetation cover using the Normalized Difference Vegetation Index (NDVI), 3) reclassifying the GDI

(8)

and the NDVI maps 4) correlating the GDI to the changes in NDVI. In the third phase the deliverables will be presented which consist of the geodiversity map, the vegetation change map and three

correlation matrices.

Figure 3. The workflow for assessing geodiversity and vegetational changes in Kennemerland-Zuid.

2.3 Defining geodiversity

Geodiversity is a new, hard to grasp concept and there is no one-size-fits-all method agreed upon by a majority of researchers (Gray, 2008; dos Santos et al., 2019). There have been multiple methodologies proposed for assessing geodiversity using a quantitative approach focusing on the spatial diversity (Gonçalves et al., 2020). However, researchers have found it difficult to find a replicable methodology for different areas as the scale of the analysis, study area dimensions and availability of cartographic data are problems that arise (Gonçalves et al., 2020; Zwoliński et al., 2018).

This study has chosen for a quantitative index-based approach to map the spatial variability of geodiversity elements through sub-indices. A geodiversity index map will be constructed similar to the method used in Seijmonsbergen et al. (2018). The sub-indices that were included in the construction of the GDI equation were chosen based on Seijmonsbergen et al. (2018) and availability of the data that could be found online. Consequently the GDI consists of the topographical diversity index (Tdi), the geomorphological diversity index (Gmdi) and the pedological diversity index (Pdi). As the study area is too small for multiple geological features, geological diversity was not taken into account for the

(9)

GDI. The hydrological diversity was also omitted because the research focuses on the relationship between geodiversity and vegetation cover (which hydrological features do not have).

The resulting GDI equation is as followed;

𝐺𝐷𝐼 = 𝑇𝑑𝑖 + 𝐺𝑚𝑑𝑖 + 𝑃𝑑𝑖 (1)

2.4 Pre-processing

2.4.1 Collection of data

The first step in the workflow was finding and selecting the necessary data from freely available sources. The data that was used is shown in Table 1 and includes Landsat imagery, the Algemeen Hoogtebestand Nederland 3 (AHN3), a geomorphology map, a pedology map, a hydrology map and a map of the Natura 2000-area boundaries. These maps can also be found in Appendix A. Table 1. Metadata of the datasets

Layer Data description Data type Coordi nate system Scale/ cell size Publication date Source

Landsat images Bands 3 & 4 from the Landsat TM 5 Bands 4 & 5 from Landsat OLI 8 Raster EPSG 32631

30m 2010 & 2018 U.S. Geological Survey & EarthExplorer (2010) U.S. Geological Survey & EarthExplorer (2018) Topography (AHN3) Surface elevation (0.5m DTM LiDAR) Raster EPSG 7415 0.5m 2018 Stuurgroep AHN & Rijkswaterstaat Geomorphology Geomorphology map Polygon EPSG 28992 1:50.000 2019 Wageningen Environmental Research GeoDesk Pedology Soil type map Polygon EPSG

28992

1:50.000 2014 - 2018 Wageningen Environmental Research GeoDesk Hydrology Waterbody map Polygon EPSG

28992

1:10.000 2015 Planbureau voor de

Leefomgeving Area boundaries Natura

2000-areas boundaries Polygon EPSG 28992 1:10.000 2019 Waterschap Hollandse Delta & EEA

(10)

2.4.2 Coordinate system & grid creation

The second step in the workflow consists of three parts. First, all collected data was uploaded into ArcGIS Pro after which the Landsat and AHN3 data was converted to the same projection RD New (EPSG 28992) using the Project tool. This coordinate system was chosen because it is the commonly used projection for the Dutch government and most Dutch research institutions.

Secondly, from the Natura 2000-areas dataset (Waterschap Hollandse Delta & EEA, 2019) the Kennemerland Zuid polygon was selected and using the Erase tool the waterbody polygons

(Planbureau voor de Leefomgeving, 2015) were removed from the study area extent to prevent surface water being part of vegetation cover calculation.

The third part was the creation of a grid to calculate the GDI. An optimal grid size had to be chosen and according to Hengl (2006) this should be a compromise between the coarsest and finest resolution, taking into account the original scale and properties of the datasets. This method leads to the following equations:

Finest legible resolution:

𝑃 ≥ 𝑆𝑁 ∗ 0.0001 (2)

Coarsest legible resolution:

𝑃 ≤ 𝑆𝑁 ∗ 0.0025 (3)

Recommended compromise:

𝑃 = 𝑆𝑁 ∗ 0.0005 (4)

Where:

p = grid resolution and SN = scale number

When taking the largest scale number of the different datasets (1:50.000), SN would be 50,000. This would mean the finest grid size would be 5m, the coarsest grid size 125m and the recommended compromise grid size 25m. Since the Landsat images have a resolution of 30m the working grid size was decided to be 90m to fit exactly nine pixels inside one grid cell. This grid size lies in between the recommended and coarsest grid size as recommended by Hengl (2006) and additionally maintains the variety within the sub-indices. Further testing with different sizes showed diversity data was lost or inaccurate.

A fishnet grid was created with a cell size of 90m x 90m using the Create Fishnet tool and clipped to the extent of the study area using the Clip tool. The cells intersecting the borders of the study area were removed from the grid to prevent reshaping of these cells using the Select Layer By Location tool.

2.4.3 Clipping, merging and rasterization of layers

The third step in the workflow consists of several steps to clip and rasterize the data for the sub-indices and to merge and clip the Landsat imagery for the vegetation cover analysis. The geomorphology map (Wageningen Environmental Research GeoDesk, 2018) and pedology map (Wageningen Environmental Research GeoDesk, 2019) were clipped by means of the Clip tool and

(11)

using the extent of the study area. The ten AHN3 raster layers (Stuurgroep AHN & Rijkswaterstaat, 2018) were merged together to cover the study area using the Mosaic to New Raster tool, which was also used for the two Landsat 5 raster bands (U.S. Geological Survey & EarthExplorer 2010)

as these both covered only a part of the study area. These merged rasters and the two Landsat 8 raster bands (U.S. Geological Survey & EarthExplorer 2018) were subsequently clipped by means of the Clip Raster tool and using the extent of the study area. The geomorphological map did not cover the whole study area which meant two small polygons were manually added to fill up these gaps. After comparing the satellite imagery inside of the two polygons and comparing this to the surrounding imagery the polygons were classified accordingly.

As the geomorphology and pedology map have been downloaded as shapefiles the data needed to be converted to raster data for the grid analysis. The tool Feature to Raster was used to transform the feature layer shapefiles (soil types and landforms) to raster data layers with cells that are 30m by 30m. This resolution was chosen after tests with different scales because it is in between the coarsest and recommended legible resolution as proposed by Hengl (2006), ensures the data is accurate and is the same size as the Landsat raster bands.

2.5 Analysis

2.5.1 Calculation of sub-indices & NDVI

A. Sub-indices

The fourth step of the workflow was to calculate the diversity sub-indices the geomorphology and pedology rasters were used as input value raster for the Zonal Statistics tool where the fishnet grid was used as input raster. This tool calculates the number of unique values (variety) per grid cell (90m x 90m) to create diversity rasters.

Topographic complexity often is an indicator for topographic diversity and complexity can be measured through topographic variables such as slope and elevation. The authors of Seijmonsbergen et al. (2018) found that slope and elevation are highly correlated and therefore this study will just use the standard deviation of the slope as a measure for topographical diversity. This was done by using the Slope tool and subsequently the Resample tool to obtain a slope raster with the same spatial resolution as the geomorphological diversity and pedological diversity rasters. The Zonal Statistics tool was then also used for the slope data to calculate the standard deviation per grid cell.

The sub-indices were then each classified into three classes ranging from ‘Low’ to ‘High’ and named the geomorphological diversity index (Gmdi) and pedological diversity index (Pdi). Three classes were chosen because the variance per grid cell for both geomorphology and pedology did not exceed three features. The topographical diversity was however reclassified to the topographical diversity index (Tdi), as the standard deviation of the slope per cell ranged from 0m to above 18m (diversity maps are located in Appendix B). The reclassification was done through the Reclassify tool and the Natural Breaks method (Jenks, 1967). This method ensures that the variance within the classes is minimalized while the variance between the classes is maximized (Jenks, 1967). This method showed more variety than the equal interval method. The details on how the sub-indices were constructed and the values per grid cell are shown in Table 2.

(12)

Table 2. Values per grid cell per sub-index

Sub-index Derived from Classification Value per grid cell Tdi AHN3 DTM (Stuurgroep AHN &

Rijkswaterstaat); Standard deviation

Low (1) Medium (2) High (3) ≤ 4,48m ≤ 7,58m ≤ 18,18m Gmdi Landform type (Wageningen

Environmental Research GeoDesk, 2018); Variety Low (1) Medium (2) High (3) 1 feature type 2 feature types 3 feature types Pdi Soil types (Wageningen

Environmental Research GeoDesk, 2019); Variety Low (1) Medium (2) High (3) 1 feature type 2 feature types 3 feature types

B. NDVI

To measure and visualize the change in vegetation cover between 2010 and 2018 the NDVI for each year has to be calculated. To do this the Raster Calculator tool was used where the input is the NDVI equation;

NDVI = 𝑁𝐼𝑅 − 𝑅𝑒𝑑

𝑁𝐼𝑅 + 𝑅𝑒𝑑 (5)

Where:

NDVI = Normalized Difference Vegetation Index NIR = Near Infrared band

Red = Visible Red band

To calculate the NDVI in 2010 the Band 4 raster layer is the Near-Infrared band and the Band 3 raster layer is the Visible Red band. To calculate the NDVI in 2018 the Band 5 raster layer is the Near-Infrared band and the Band 4 raster layer is the Visible Red band. All steps taken to calculate the sub-indices and the NDVI rasters are located in Appendix C.

2.5.2 Calculating the GDI and NDVI change

A. GDI calculation & classification

To calculate the GDI, the sum of the values of the sub-indices according to Eq. (1) are calculated per grid cell using the Raster Calculator tool. The values for the GDI range between three and seven as no grid cells have a minimum value of three and no cells have a value of eight or nine. These five values are reclassified using the Reclassify tool and manually named ‘Very low’, ‘Low’, ‘Medium’, ‘High’ and ‘Very high’.

B. NDVI change calculation & classification

To calculate the temporal change of the NDVI in the study area the 2010 NDVI raster was subtracted from the NDVI 2018 raster using the Raster Calculator tool. This NDVI change raster was then transformed to the same extent as the GDI raster through the use of the Zonal Statistics tool and the fishnet. This raster shows the mean change in NDVI value per grid cell. This grid was then reclassified using the Natural Breaks method (Jenks, 1967) and slightly adjusting the values so that

(13)

there is no class with both negative and positive values for change. Table 3 shows the NDVI change classes and their upper values.

Table 3. NDVI change classification data

Map Derived from Class (new value) Upper value per

grid cell NDVI

Change

Much vegetation loss (1) Vegetation loss (2)

Little vegetation loss to no change (3) Vegetation gain (4)

Much vegetation gain (5)

≤ -0,114 ≤ -0,063 ≤ 0,00001

≤ 0,024 ≤ 0,484

2.6 Deliverables & results

2.6.1 Geodiversity map

The geodiversity map is based on the GDI and is shown in Figure 4. Constructed through the summation of three sub-indices, the map gives a rasterized overview of geodiversity ranges between ‘Very high’ and ‘Very low’. The map and the areal percentages of the five geodiversity classes that are listed in Table 4, indicate that the study area predominantly consists of ‘Very low’ to ‘Low’

geodiversity area (67,84%). However, the northwestern region, the edges near the beach and a small area next to Overveen display a bit higher geodiversity than the rest while the southeastern area is much less geodiverse.

(14)

Figure 4. The geodiversity map of the Kennemerland Zuid study area, classified into five categories. Table 4. Geodiversity classes with their value, number of grid cells, area and percentage

Geodiversity class Value Grid cells Area (ha) Percentage area (%)

Very low 1 2132 1726,92 27,65 Low 2 3099 2510,19 40,19 Medium 3 1961 1588,41 25,43 High 4 470 380,7 6,10 Very high 5 48 38,88 0,62 Total - 7710 6245,1 100

The correlation between the GDI and the sub-indices are given in Table 5 and calculated through the Band Collection Statistics tool. This tool gives a correlation matrix that shows the values of the correlation coefficients that depict the relationship between cell values from one layer to cell values from another layer. From this matrix it can be concluded that all three sub-indices show positive correlation with the GDI. The Tdi clearly has the strongest correlation with the GDI (0.836)

(15)

with the Pdi being a bit less correlated (0.295) and the Gmdi showing the least correlation. That the Tdi influences the GDI the most can also be seen when looking at the Tdi map (Appendix B1) where there is a high topographical diversity in the northwest and low diversity in the southeast. The correlation between the individual sub-indices is very small and thus insignificant.

Table 5. The geodiversity correlation matrix with sub-indices

GDI Topographical Index (Tdi) Geomorphological Index (Gmdi) Pedological Index (Pdi) GDI 1.000 Topographical Index (Tdi) 0.836 1.000 Geomorphological Index (Gmdi) 0.295 -0.016 1.000 Pedological Index (Pdi) 0.517 0.0391 0.107 1.000

2.6.2 NDVI change map

The NDVI change map is shown in Figure 5 and depicts the temporal changes of vegetation cover between 2010 and 2018. The rasterized map shows the changes in NDVI value per grid cell that ranges between ‘much vegetation loss’ and ‘much vegetation gain’. The map and the data depicted in Table 6 show that most cells show vegetation losses to no change (78,89%). The most vegetation loss occurred in the southeastern and middle regions while most vegetation growth is seen in the northwest, southern tip and at the edges of the study area.

(16)

Figure 5. The NDVI change map of the Kennemerland Zuid study area, classified into five categories. Table 6. Change of the NDVI depicted in classes with their value, number of grid cells, area and percentage

NDVI change class Value Grid cells Area (ha) Percentage area (%)

Much vegetation loss 1 2132 1726,92 4,79

Vegetation loss 2 3099 2510,19 23,81

Little vegetation loss to no change 3 1961 1588,41 55,08

Vegetation gain 4 470 380,7 11,69

Much vegetation gain 5 48 38,88 4,63

Total - 7710 6245,1 100

(17)

2.6.3 GDI & NDVI change correlation

Analysing the relationship between the GDI and the NDVI change is the goal of this study. To achieve this the Band Collection Statistics tool is used to calculate the relationship between the GDI raster layer and the NDVI change layer. The results of this tool are displayed in Table 7 and shows that these layers have a correlation coefficient of 0.063 which signifies a very small correlation. Table 7. The correlation matrix between GDI and NDVI change

GDI NDVI Change

GDI 1.000

(18)

3. Discussion

3.1 Discussion of the results

The main research question of this study is: how is the geodiversity in the Natura 2000-area Kennemerland-Zuid related to the changes in vegetation cover between 2010 and 2018?. To answer this question two sub-questions were derived, from which the first question is: what is the geodiversity of the study area?.

This question is answered in the form of the GDI map (Figure 4) and diversity sub-index maps (Appendix B). The geodiversity of the study area can be explained through the correlation matrix (Table 5). First of all, it is notable that the Tdi correlates most with the GDI. This is caused by the fact that in general there is much more topographical diversity than there is geomorphological and

pedological diversity (Appendix B). As there are only a limited amount of landforms (8) and soil types (11) (Appendix A) in the study area, most grid cells cover only one feature and thus have a low diversity value. For example, the majority of the geomorphological map is made up of ‘coastal dunes with associated plains and depressions’. This also explains why the majority of the area consists of low GDI values. At the same time there is a lot of slope variation (Appendix A1) which means there are much more cells with a ‘Medium’ to ‘High’ Tdi value. This means that the GDI is relatively more influenced by the Tdi values than by the lower Gmdi and Pdi values.

The lower GDI in the south and north-east are the result of the fact that this region is older and mainly consists of older more gentle dunes and planar dune valleys (Provincie Noord-Holland, 2018) which means less slope variety. The dunes are a gradient-ecosystem where the horizontal gradient, depending on the distance to the shore, is directly related to the soil salinity, calcium concentration and age of the terrain and other environmental factors (Figure 6). The GDI is higher near the beach

because embryotic dunes and foredunes are located here, these dunes are steeper and higher than older inland dunes (Hesp, 2002; Muñoz-Vallés & Cambrollé, 2016). The high GDI hotspot next to

Overveen can be explained by the fact that there are three different soil types, two landforms and steep slopes.

Figure 6. Main environmental gradients on coastal dunes (Muñoz-Vallés & Cambrollé, 2016). That the GDI in the study area is relatively low is not something unexpected. This study is done on a local scale with a relatively small grid resolution which generally leads to lower GDI values (Zarnetske et al., 2019). The impact of grid resolution on the outcome of the GDI is large and a change in resolution would lead to significantly different results (Erikstad, 2013; Hengl, 2006).

(19)

The second sub-question reads: where do the changes in vegetation cover occur within the study area?. This question is answered through visualization of the change in NDVI values between 2010 and 2018 in a rasterized map (Figure 5). The majority of the NDVI change values are between -0,114 and zero (Tables 3 & 6) which indicates there is mostly no change in vegetation cover. That most vegetation loss occurred in the southeast, middle and some spots in the north (Figure 5) can be explained by several factors. Most hotspots where there has been ‘much vegetation loss’ are due to conservational interventions like opening the foredunes and removing shrub vegetation for increased aeolian sand transportation as part of the ‘third generation dune dynamics experiments’ (Kooijman et al., 2012; Kuipers, et al., 2016; van der Valk & Arens, 2013). Other conservational practices have also been focused on removing vegetation by means of grazers and nature bridges. The vegetation loss east of Zandvoort is for example a direct result of these bridges where animals cross and remove vegetation (Figure 7).

Figure 7. A comparison of two satellite images (2010 & 2018) of the area around the Zandvoort nature bridge.

Only 16,32% of the study area experiences some type of vegetation growth (Table 6) and the spots where this happens are mostly located on the edges. A likely possible reason for this increase in vegetation is the increased nitrogen deposition coming from highways, industry, agriculture, ships on the Noordzeekanaal (Provincie Noord-Holland, 2018).

To answer the main research question the correlation between the GDI and the NDVI change was measured as 0.063 (Table 7). This signifies a very weak positive relationship and is negligible. The result that there is no relationship between the GDI and NDVI change was not expected but can be explained by looking at the sub-questions and the method. First of all, many changes in vegetation cover can be linked back to human interventions on locations that were chosen and not naturally selected, which makes a relationship between natural abiotic features and vegetation cover harder to find.

Additionally, as the GDI is largely dependent on the Tdi, it can be assumed that the NDVI changes are also not correlated to the Tdi. A possible explanation for this is that in rough dune areas with a steeper slopes, the vegetation mainly consists of small grasses and succession/change of this

vegetation cover is slow and minimal (Arens, 2018; Provincie Noord-Holland, 2018). Earlier research has also shown that vegetational patterns in dunes are largely influenced by environmental gradients (Stallins, & Parker, 2003) and it can be hypothesized that these gradients are of more influence than geodiversity. However, future studies and field surveys should be implemented to control and study if this effect is true.

Nature bridge, Zandvoort

(20)

3.2 Discussion of the method

Geodiversity is a relatively new field of research and further development of methodologies to quantify and study geodiversity are much needed (Alahuhta et al., 2018; Erikstad, 2008). This study uses the method used in Seijmonsbergen et al. (2018) to quantify geodiversity.

However, as mentioned before in section 4.1 and mentioned by previous authors (Erikstad 2013; Serrano & Ruiz-Flaño, 2007),the scale of the study matters a lot for the results. The small study area and thus relatively small grid size has consequences for the results of this study. No geological diversity was measured because of the scale choice. Fortunately this study was done through use of a model and the size of the grid is easily adjustable. At the same time hydrological diversity was omitted due to the vegetation cover analysis. This has probably influenced the results because hydrological features have previously been found to have significant influence on the geodiversity (Hjort & Luoto, 2012; Tukiainen et al., 2017). Additionally, the removal of boundary cells from the fishnet grid has caused an incomplete analysis of the study area.

Another limitation of this study was the accuracy of the data. The geomorphological and hydrological features did not always align with the landforms and waterbodies that were seen on satellite imagery during the research. Additionally, the input data maps were published in between 2010 and 2018 which means the temporal scale of the GDI differs from that of the NDVI maps. However, geodiversity is relatively stable on short timescales and the topographical diversity was found to be a reliable predictor of the GDI, as seen in previous studies (Hjort & Luoto, 2012; Seijmonsbergen et al., 2018). It is also true that when comparing vegetation cover maps produced at different times to detect change, the errors from the classification process (that can occur due to shadows/surface water) will show up as change (Zhu & Woodcock, 2014). It is recommended that future studies use data that is as accurate as possible.

Changes in vegetation were measured through the NDVI because it is a easy-to-use and common tool to do such an analysis. This also comes with limitations however; as coastal dune vegetation is short and sparse and generally less noticeable through remote sensing (Yousefi Lalimi et al., 2017). Additionally, the two NDVI maps that were compared were both taken only a few days apart but the summer of 2018 was the hottest Dutch summer in three centuries (KNMI & Huiskamp, 2018) with low precipitation all spring and summer. This has probably influenced the development of vegetation and thus the NDVI change map. Future studies can take this into account by comparing more NDVI rasters. Another shortcoming of using the NDVI is that it just measures vegetation cover and does not take into account the differences between separate plants. When assessing the relationship between geodiversity and vegetation, researchers might want to look at the influence geodiversity has on individual plant species through means of field study.

(21)

4. Conclusion

Geodiversity is a relatively new concept which can be used for conservational practices all over the world. The link between geodiversity and biotic elements is especially relevant for this purpose. This research has aimed to find such a link by calculating the correlation between a rasterized geodiversity map and rasterized Normalized Difference Vegetation Index change map within the Natura 2000-area Kennemerland Zuid.

The method used in this paper is build up of three parts; pre-processing, analysis and results. The method is easily adjusted and can be used for other areas because of the model workflow used in ArcGIS Pro. The geodiversity was measured using a 90m by 90m grid and summing the sub-indices of the geodiversity elements within each grid cell. These geodiversity elements were topography, pedology and geomorphology. Five geodiversity classes were made for the geodiversity index (GDI) ranging from ‘Very low’ to ‘Very high’. Additionally the vegetation cover was calculated by subtracting the 2010 NDVI from the 2018 NDVI in the same grid size as the GDI.

It can be concluded from the correlation matrix that no significant relationship between the GDI and NDVI change was found. This is probably related to 1) the changes in vegetation in the study area are mainly caused by anthropogenic factors such as nitrogen deposition and removal of vegetation for conservation, and 2) vegetation cover is mainly dependent on the horizontal environmental gradients and less on geodiversity elements such as topographical diversity. Additionally, this method has some shortcomings of which one is the local scale and thus low relative geodiversity. Although this research has its limits and did not find any correlation between the GDI and the NDVI changes, the workflow used in this paper provides further elaboration on quantitative remotely sensed geodiversity and the relationship with biotic features, as this study is only a small step in the further development and use of GDI assessments. .

(22)

References:

Alahuhta, J., Ala-Hulkko, T., Tukiainen, H., Purola, L., Akujärvi, A., Lampinen, R., & Hjort, J. (2018). The role of geodiversity in providing ecosystem services at broad scales. Ecol. Indic., 91, 47– 56. doi: 10.1016/j.ecolind.2018.03.068

Anfuso, G., & Del Pozo, J. Á. M. (2009). Assessment of coastal vulnerability through the use of GIS tools in South Sicily (Italy). Environmental Management, 43(3), 533-545.

Araujo, A. M., & Pereira, D. Í. (2018). A new methodological contribution for the geodiversity assessment: applicability to Ceará State (Brazil). Geoheritage, 10(4), 591-605.

Arens, S. M., Mulder, J. P., Slings, Q. L., Geelen, L. H., & Damsma, P. (2013). Dynamic dune management, integrating objectives of nature development and coastal safety: examples from the Netherlands. Geomorphology, 199, 205-213.

Arens, B. (2018) Monitoring Geomorfologie 2013-2017 Eindrapport.Waternet, Arens BSDO RAP2018.05

van den Berg, C. (2021). Nationaal Park Zuid-Kennemerland [Photo]. Retrieved from https://www.droomplekken.nl/nederland/nationaal-park-zuid-kennemerland/

Brilha, J., Gray, M., Pereira, D. I., & Pereira, P. (2018). Geodiversity: An integrative review as a contribution to the sustainable management of the whole of nature. Environmental Science & Policy, 86, 19-28.’

Chakraborty, A., & Gray, M. (2020). A call for mainstreaming geodiversity in nature conservation research and praxis. Journal for Nature Conservation, 56, 125862. doi: 10.1016/j.jnc.2020.125862 Comer, P. J., Pressey, R. L., Hunter, M. L., Schloss, C. A., Buttrick, S. C., Heller, N. E., ...Shaffer, M. L. (2015). Incorporating geodiversity into conservation decisions. Conserv. Biol., 29(3), 692–701. doi: 10.1111/cobi.12508

Delgado-Fernandez, I., O'Keeffe, N., & Davidson-Arnott, R. G. (2019). Natural and human controls on dune vegetation cover and disturbance. Science of The Total Environment, 672, 643-656. Erikstad, L. (2008). History of geoconservation in Europe. Geological Society, London, Special Publications, 300(1), 249-256.

Erikstad, L. (2013). Geoheritage and geodiversity management–the questions for tomorrow. Proceedings of the Geologists' Association, 124(4), 713-719.

Esri Inc. (2020). ArcGIS Pro (Version 2.5). Esri Inc. https://www.esri.com/en-us/arcgis/products/arcgis-pro/overview.

Gonçalves, J., Mansur, K., Santos, D., Henriques, R., & Pereira, P. (2020). A Discussion on the Quantification and Classification of Geodiversity Indices Based on GIS Methodological Tests. Geoheritage, 12(2), 1–20. doi: 10.1007/s12371-020-00458-3

(23)

Gray, M. (2008). Geoheritage 1. Geodiversity: a new paradigm for valuing and conserving geoheritage. Geoscience Canada.

Hengl T. (2006). Finding the right pixel size. Computers & Geosciences 32: 1283–1298. DOI: 10.1016/j.cageo.2005.11.008.

Hesp, P. (2002). Foredunes and blowouts: initiation, geomorphology and dynamics. Geomorphology, 48(1), 245–268. doi: 10.1016/S0169-555X(02)00184-8

Hjort, J., & Luoto, M. (2010). Geodiversity of high-latitude landscapes in northern Finland. Geomorphology, 115(1), 109–116. doi: 10.1016/j.geomorph.2009.09.039

Hjort, J., & Luoto, M. (2012). Can geodiversity be predicted from space?. Geomorphology, 153, 74-80.

Hjort, J., Gordon, J. E., Gray, M., & Hunter Jr, M. L. (2015). Why geodiversity matters in valuing nature's stage. Conservation Biology, 29(3), 630-639.

Jenks, G. F. (1967). The Data Model Concept in Statistical Mapping. International Yearbook of Cartography, 7, 186-190.

KNMI, & Huiskamp, A. (2018, September 3). KNMI - Zomer 2018 (juni, juli, augustus) [Press release]. Retrieved June 1, 2021, from https://www.knmi.nl/nederland-nu/klimatologie/maand-en-seizoensoverzichten/2018/zomer

Kooijman, A., van der Hagen, H., & Noordijk, E. (2012). Stikstof depositie in de duinen: alles in beeld? Landschap, 29(3), 147-154.

http://www.landschap.nl/images/stories/landschap/2012/editie3/2012-3_147-154.pdf Koster, E. A. (2009). The “European Aeolian Sand Belt”: geoconservation of drift sand landscapes. Geoheritage, 1(2-4), 93-110

Kuipers, M., Arens, S. M., & Ruessink, B. G. (2016). Grootschalig herstel van stuivende duinen. De Levende Natuur, 117(003), 90-93.

Łabuz, T. A. (2016). A review of field methods to survey coastal dunes—experience based on research from South Baltic coast. Journal of Coastal Conservation, 20(2), 175-190.

Lawler, J. J., Ackerly, D. D., Albano, C. M., Anderson, M. G., Dobrowski, S. Z., Gill, J. L., ...Weiss, S. B. (2015). The theory behind, and the challenges of, conserving nature's stage in a time of rapid change. Conserv. Biol., 29(3), 618–629. doi: 10.1111/cobi.12505

Martinez, M. L., Gallego-Fernandez, J. B., Garcia-Franco, J. G., Moctezuma, C., & Jimenez, C. D. (2006). Assessment of coastal dune vulnerability to natural and anthropogenic disturbances along the Gulf of Mexico. Environmental Conservation, 109-117

Maul, G. A. (1993). Climatic Change in the Intra-Americas Sea: Implications of Future Climate on the Ecosystems and Socio-economic Structure in the Marine and Coastal Regions of the Caribbean Sea, Gulf of Mexico, Bahamas and the Northeast Coast of South America. Edward Arnold.

Muñoz-Vallés S., Cambrollé J. (2016) Dunes. In: Harff J., Meschede M., Petersen S., Thiede J. (eds) Encyclopedia of Marine Geosciences. Encyclopedia of Earth Sciences Series. Springer, Dordrecht. https://doi.org/10.1007/978-94-007-6238-1_160

(24)

Nicholls, R. J., & Cazenave, A. (2010). Sea-level rise and its impact on coastal zones. science, 328(5985), 1517-1520.

Planbureau voor de Leefomgeving. (2015). Basiskaart Aquatisch: de Watertypenkaart 2013 | Data overheid [De basiskaart aquatisch watertypen is een gedetailleerde dataset met de ligging van het Nederlandse oppervlaktewater.]. Retrieved from https://data.overheid.nl/dataset/13653-basiskaart-aquatisch--de-watertypenkaart-2013

Provincie Noord-Holland. (2018, October). N2000 beheerplan Kennemerland-Zuid 2018–2024. Retrieved from

https://www.noordholland.nl/Onderwerpen/Natuur/Projecten/Stikstof_en_Natura_2000/Natura_2000_ beheerplannen/Beheerplannen/N2000_beheerplan_Kennemerland_Zuid_2018_2024

dos Santos, D. S., Mansur, K. L., de Arruda Jr, E. R., Dantas, M. E., & Shinzato, E. (2019). Geodiversity mapping and relationship with vegetation: A regional-scale application in SE Brazil. Geoheritage, 11(2), 399-415.

Seijmonsbergen, A. C., Guldenaar, J., & Rijsdijk, K. F. (2018). Exploring Hawaiian long-term insular geodiversity dynamics. Landform Analysis, 35, 31-43.

https://doi.org/10.12657/landfana.035.007

Serrano, E., & Ruiz-Flaño, P. (2007). Geodiversity: a theoretical and applied concept. Geographica Helvetica, 62(3), 140-147.

Stallins, J. A., & Parker, A. J. (2003). The Influence of Complex Systems Interactions on Barrier Island Dune Vegetation Pattern and Process. Ann. Assoc. Am. Geogr., 93(1), 13–29. doi: 10.1111/1467-8306.93102

Stuurgroep AHN & Rijkswaterstaat. (2018). Actueel Hoogtebestand Nederland 3 (AHN3) [Het Actueel Hoogtebestand Nederland 3 (AHN3) zijn hoogtebestanden en hulpbestanden vervaardigd uit laseraltimetrie in de periode 2013–2018.]. Retrieved from

https://www.nationaalgeoregister.nl/geonetwork/srv/dut/catalog.search#/metadata/41daef8b-155e-4608-b49c-c87ea45d931c?tab=general

Tukiainen, H., Bailey, J. J., Field, R., Kangas, K., & Hjort, J. (2017). Combining geodiversity with climate and topography to account for threatened species richness. Conservation Biology, 31(2), 364-375.

U.S. Geological Survey & EarthExplorer. (2010). Landsat 5 TM [Landsat 5 TM raster bands 3 & 4]. Retrieved from https://earthexplorer.usgs.gov/

U.S. Geological Survey & EarthExplorer. (2018). Landsat 8 OLI [Landsat 8 OLI raster bands 4 & 5]. Retrieved from https://earthexplorer.usgs.gov/

van der Valk, B., & Arens, S. M. (2012). Stuivende kustduinen : illusie of werkelijkheid? |

Groenekennis. Holland’s Duinen: Informatie over Het Duinonderzoek in Berkheide, Meijendel En Solleveld, (2012)60. Retrieved from https://library.wur.nl

Wageningen Environmental Research - GeoDesk. (2014–2018). BRO - Bodemkaart van Nederland 1:50.000 [Deze dataset toont verschillende aspecten van de geactualiseerde Bodemkaart van Nederland, schaal 1 : 50 000]. Retrieved from

(25)

https://www.nationaalgeoregister.nl/geonetwork/srv/api/records/ed960299-a147-4c1a-bc57-41ff83a2264f

Wageningen Environmental Research (WENR) - GeoDesk. (2019). BRO - Geomorfologische Kaart van Nederland 1: 50:000 [De Geomorfologische Kaart van Nederland 1 : 50 000 beschrijft de in Nederland voorkomende geomorfologische eenheden en hun ontstaanswijze.]. Retrieved from https://www.nationaalgeoregister.nl/geonetwork/sitemap/api/records/459231d0-7379-4f26-a444-7616e1d888f0

Waterschap Hollandse Delta & EEA. (2019). Natura 2000 locaties [Nederland kent ruim 160 Natura 2000-gebieden. Dit Natura 2000-netwerk bestaat uit gebieden die zijn aangewezen onder de

Vogelrichtlijn en onder de Habitatrichtlijn.]. Retrieved from

https://bio.discomap.eea.europa.eu/arcgis/rest/services/ProtectedSites/Natura2000Sites/MapServer Yousefi Lalimi, F., Silvestri, S., Moore, L. J., & Marani, M. (2017). Coupled topographic and vegetation patterns in coastal dunes: Remote sensing observations and ecomorphodynamic implications. Journal of Geophysical Research: Biogeosciences, 122(1), 119-130.

Zarnetske, P. L., Read, Q. D., Record, S., Gaddis, K. D., Pau, S., Hobi, M. L., ... & Finley, A. O. (2019). Towards connecting biodiversity and geodiversity across scales with satellite remote sensing. Global Ecology and Biogeography, 28(5), 548-556.

Zhu, Z., & Woodcock, C. E. (2014). Continuous change detection and classification of land cover using all available Landsat data. Remote sensing of Environment, 144, 152-171.

Zwoliñski, Z., Najwer, A., Giardino, M., (2018) Methods for assessing geodiversity. In: Reynard E, Brilha J (eds) Geoheritage: assessment, protection, and management. Elsevier, pp 27–47

(26)

Appendices

Appendix A: Input maps

(27)

(28)

(29)

(30)

Appendix B: Diversity maps of the sub-indices & NDVI maps

(31)

(32)

(33)

(34)

(35)

Appendix C: Models used in ArcGIS Pro

Appendix C1: Pre-processing model for fishnet

(36)

Appendix C3: Analysis & reclassification of the sub-indices

(37)

Unravelling the relationship between vegetation change and geodiversity in the Dutch dunes: a case study