A quantitative assessment of geodiversity and human influence in geoparks

(1)

1

A quantitative assessment of

geodiversity and human influence

in geoparks

MSc thesis Emma Polman - 10799478 9/11/2020

Supervisor dr. Harry Seijmonsbergen Co-assessor dr. rer. nat. Daniel Kissling

(2)

2

Abstract

UNESCO global geoparks aim at protecting globally significant geoheritage and geodiversity using a management strategy that relies upon sustainable development for the local population. However, both geodiversity and human influence in geoparks have never been quantified and assessed in a global context and the extent to which geodiversity and human influences are represented in geoparks is therefore unknown. Here, geodiversity, human influence and human influence change between 1993 and 2009 are quantified in 147 geoparks and compared to global Asian and European geodiversity, human influence and human influence change. In addition, the representation of soil and lithology types in geoparks was assessed. The results show that the total geodiversity and lithological and topographic diversity were in general significantly higher in geoparks, while soil diversity was not significantly higher and hydrological diversity was lower. This is most likely caused by the emphasis on geology and geomorphology in the UNESCO application procedure and the large share of geoparks located in mountainous areas, where the lithological and topographic diversity and subsequently the total geodiversity are high. Igneous and volcanic rocks and andosols were represented best in geoparks, especially in Europe. Of all soil and lithology types 22 and 65% were not represented in any geopark, respectively, likely due to the uneven global distribution of geoparks. The human influence indices were higher in geoparks indicating that there is a substantial amount of human activities present in geopark areas. Human influence change was not significantly different in geoparks, compared to global, Asian and European human influence change. It is recommended that future geoparks put more emphasis on soil and hydrological diversity and on the soil and lithology types that are not or hardly represented in geoparks to better reflect global geodiversity. Human activities in geoparks should be monitored in order to assess if damage is caused to geodiversity.

List of abbreviations

GI geodiversity index GGI global geodiversity index HF human footprint

(3)

3

1. Introduction

The conservation of geological heritage has become of increased importance over the past years. Awareness has risen that geological features hold valuable information about past earth processes and that their loss is irreversible given the large time scales associated with their formation (Gray, 2008; Brilha, 2016). Related to this increasing valuation of the earth's abiotic features is the concept of geodiversity. The term geodiversity was introduced in the 1990s (Gray, 2008) and the now most widely accepted definition is the "natural range (diversity) of geological (rocks, minerals, fossils), geomorphological (landforms, topography, physical processes), soil and hydrological features. It includes their assemblages, structures, systems and contributions to landscapes" (Gray, 2013). Geodiversity is now increasingly used within the scientific community (Ibáñez et al., 2018; Zwoliński et al., 2018).

An attempt to conserve globally significant geoheritage and geodiversity is through the establishment of UNESCO Global Geoparks. From 2001-2019, UNESCO has awarded the status of "Global Geopark" to 147 areas. Geoparks are defined as "single unified geographical regions where sites and landscapes of international geological significance are managed with a holistic concept of protection, education and sustainable development" (UNESCO, 2015). Geodiversity is one of the assessment criteria during the geopark application procedure (UNESCO, 2016), but most geodiversity related research in geoparks concerns qualitative analyses about the relation between geodiversity and touristic values (Boothroyd & Henry, 2019). The term "geodiversity" is only used in the official description of 41 geoparks, and many geoparks consider geodiversity as a synonym for geoheritage or the presence of unique geological features (Ruban & Yashalova, 2018), while these are only a part of geodiversity (Crofts & Gordon, 2015; Brilha, 2016). Up to this date, no systematic research has been conducted on quantifying geodiversity in geoparks and assessing its relation to global geodiversity.

The majority of geoparks are located in Europe and China (Appendix I fig 1), making it impossible to capture the entirety of different geofeatures present on earth (Ibáñez et al., 2018). Geodiversity in the geopark application procedure is only assessed as the number of different rock types, geomorphological features and geological time periods represented in the park (UNESCO, 2016), and the main focus is thus on lithological and geomorphological diversity. Other abiotic features such as soils and hydrology are not taken into account, and neither are the spatial distribution and assemblage of different features, while these are often considered in geodiversity indices (e.g. Cañadas & Ruiz-Flaño, 2007; Hjort & Luoto, 2010; Seijmonsbergen et al., 2018). Geoparks also focus on protecting geosites, sites with a high scientific, cultural educational or aesthetic value. Though Gray (2008) states that geosites aim at being representative for their abiotic surroundings, Hjort and Luoto (2010) found that geosites alone are often not capturing the total geodiversity of an area. It is therefore unlikely that geoparks capture and conserve global geodiversity.

In order to assess geodiversity, various qualitative and especially quantitative methods have been proposed in the last years (Zwoliński et al., 2018). Many studies quantify geodiversity by calculating a geodiversity index (GI) which is the sum of multiple sub indices, such as geological, hydrological, soil, topographic, geomorphological and paleontological diversity (Cañadas & Ruiz-Flaño, 2007; Pereira et al., 2013; Araujo & Pereira, 2018; Seijmonsbergen et al.,. 2018; da Silva et al.,

(5)

5 2019; Gonçalves et al., 2020; dos Santos et al., 2020). Sub indices are often calculated using a grid-based approach. Thematic maps of the research area are overlain with a grid and the diversity is determined for each grid cell by counting the number of unique geo features within the cell (e.g. Hjort & Luoto, 2010; Santos et al., 2017; Fernández et al., 2020). However, the different methods and sub indices used to calculate the GIs and the local extent of the studies make it impossible to compare GI values across areas (Ibáñez & Brevik, 2019). A global geodiversity index (GGI) based on topographic, hydrologic, lithology and soil datasets, as used by (Muellner‐Riehl et al., 2019) does allow for a transparent analysis of geodiversity. However, a GI only takes into account the variety of geo features within a grid cell, but ignores the spatial distribution of specific types of geo features, e.g. geology and soil types. Additional assessment of the types of geo features in an area is necessary to quantify their global representation. There is at the moment no standardised method for assessing geodiversity.

Geoparks are not legally protected areas, and parks encompass roads, railways, airports, urban and industrial areas, and agricultural lands. Within the parks geodiversity is also used to support sustainable development while protecting the landscape (UNESCO, 2015), but certain human activities such as urbanisation, mining and land use changes can be a potential threat to geodiversity as these can destruct landforms, disrupt geomorphological processes and cause soil erosion (Crofts & Gordon, 2015; Hjort et al., 2015). Geotourism is a main component of the sustainable development in parks (Farsani et al., 2011), but it can cause an increase in infrastructure (Shui & Xu, 2016), which can eventually lead to visitors damaging sites (Wang et al., 2015), erosion (Sumanapala & Wolf, 2020) and fragmentation of sites (Crofts & Gordon, 2015; Hjort et al., 2015). Many geoparks are in Europe and East Asia, where human influence is high (Venter et al., 2016b), but the management bodies in geoparks should prevent that human activities cause harm to geodiversity (UNESCO, 2015). There is however no global study quantifying the status of human influence and how human influence changes over time in geoparks in relation to global human influence and global human influence change.

This research aims at determining to what extent global geodiversity and human influence are represented in UNESCO global geoparks by conducting an objective quantification and assessment. This assessment will be done on the global, Asian and European extent, as the majority of geoparks are located in Asia and Europe. This research aim leads to the following four research questions: 1) how does geodiversity in geoparks compare to global, Asian and European geodiversity? 2) how do the soil and lithology types present in geoparks compare to the soil and lithology types globally, in Asia and in Europe? 3) how does human influence in geoparks compare to global, Asian and European human influence? and 4) how does human influence change in geoparks compare to global, Asian and European human influence change? It is hypothesised that 1) the lithological and topographic diversity and subsequently the total geodiversity will be higher in geoparks, while soil and hydrological diversity will be lower, because these are not explicitly assessed by UNESCO and associated less with geodiversity in the geoparks context (UNESCO, 2016; Ruban, 2017) 2) not all soil and lithology types are represented in geoparks because of the uneven spatial distribution of the parks (Ibáñez et al., 2018) 3) human influence is equally high or higher in geoparks compared to the world, Asia and Europe, since parks are established in areas with high human influence (Venter et al., 2016b; Kennedy et al., 2019) 4) human influence in geoparks is stable or decreasing in geoparks, because of the management requirements set by UNESCO.

(6)

6

2. Methods

2.1 General Framework

In order to compare geodiversity, human influence and human influence change between geoparks, the world, Asia and Europe data layers were collected, sampled and (statistically) compared (Fig 1). Additionally, a correlation analysis was part of the GI analysis to gain insight in the contribution of each GI sub layer to the total GI (GDsum) (Fig 1, H1).

Figure 1. Workflow diagram showing the three components (A-C) used to answer the four research questions and hypotheses (H1-H4). All input data layers were collected (row A) and GGI (H1), human influence (H3) and human influence change (H4) data layers were sampled in geoparks and random samples in the geoprocessing phase (row B). The soil and lithology data layers were sampled in geoparks, the world, Asia and Europe (H2). Finally, the correlations between the GDsum and GGI sub layers were calculated (H1) and the differences between geoparks and random samples were statistically tested for H1, H3 and H4 (row C). For H2 the percentage area of each soil and lithology type in geoparks relative to the global, Asian and European area was calculated.

Geopark data

A database containing all geoparks and data attributes was built using the lists of geoparks available at the Global Geoparks Network website (2019) and the UNESCO website (2019). The polygons of the geoparks were requested to all individual geoparks by mail, since these data are not publicly available. In addition, several geopark polygons were found online. This resulted in a total of 45 geopark polygons. The remaining 102 geopark polygons areas were estimated using the point location coordinates and surface area data made available by UNESCO (2019). The estimated geopark polygons were created by defining a circular buffer around each point coordinate with a radius corresponding to the geopark area. Some of these estimated polygons overlapped with the nearby sea surface, and to avoid bias since the analysis and data only concern the terrestrial surface,

(7)

7 maps of these geoparks were searched online and digitized by hand. A list of geoparks and the type of polygon used in the analysis is provided in Appendix II and the geoparks shapefile in Appendix III.

Geodiversity index data

For the GI analysis the GGI data described in the research by Muellner‐Riehl et al. (2019) were used (table 1). The data were obtained in person from Hannes Versteegh, who contributed as a co-author in that research. This dataset consists of six diversity raster layers and one raster with the summed geodiversity (GDsum). The resolution of the GGI and sub layers is 10x10 km. The six sub layers used are the number of soil types in a grid cell (SoilDiv) (Hengl et al., 2017), number of lithology types in a grid cell (LithoDiv) (Hartmann & Moosdorf, 2012), the slope range (SlopeRange) and slope standard deviation (SlopeSTD) (Yamazaki et al., 2017) in a grid cell, the total river length in a grid cell (RiverLength) (Lehner et al., 2008) and the total lake area (Lake Area) (Messager et al., 2016) in a grid cell. These diversity layers were reclassified into four sub-indices (soil, lithology, hydrology and topography) using Jenks natural breaks. The GDsum layer is the sum of these four sub indices.

Soil and lithology type data

For the analysis of lithology the Global Lithological Map (GLiM) was used (Hartmann & Moosdorf, 2012) (table 1). This polygon dataset distinguishes 15 unique main lithological types at the first legend level and 437 unique rock types at the combined first, second and third legend level. The SoilGrids (Hengl et al., 2017) (table 1) data were used for assessing the unique soil types. This raster dataset contains 118 unique soil types based on the World Reference Base (WRB) classification (IUSS Working Group, 2006). This classification includes 30 of the 32 main WRB soil types and the most common prefix and suffix qualifiers.

Human influence data

The human influence analysis included two open source datasets: the human footprint (HF), available for 1993 (HF1993) (Sanderson et al., 2002) and 2009 (HF2009) (Venter et al., 2016b) and the Global Human modification of 2016 (gHM) (Kennedy et al., 2019) (table 1). The HF maps are based on population density, electrical power infrastructure, human accessibility and land transformation, which are considered to be explicit proxies of human population and infrastructure (Venter et al., 2016a). The gHM uses similar drivers as the HF, but added higher coverage of transportation infrastructure and mining and energy activities as stressors.

Human Footprint Change

The human footprint change (HFchange) was calculated by subtracting the HF2009 raster from the HF1993 raster. The HF increased mainly in tropical areas and areas with high biodiversity. Small decreases were observed in countries with relatively high wealth and strong legislative enforcement (Venter et al., 2016b).

Continents and terrestrial outline

The continent outlines used in the research were made by merging the Esri world countries dataset (Esri, 2015) on the continent attribute. This maintained the country-level detail to the borders of the continent. Since no human influence and geodiversity data were available for Antarctica it was excluded from the analysis.

(8)

8 Table 1. Descriptions and metadata of the geodiversity and human influence datasets used in the research. Data set Data

type Cell Size (m)

Projection GCS Extent Description of data attributes

Value Range

Source

Geodiversity

Soil diversity raster 10,000 Mollweide WGS84 Global* Number of soil types in a grid cell

1-27 Muellner‐ Riehl et al. (2019)

Lithology diversity

raster 10,000 Mollweide WGS84 Global* Number of lithology types in a grid cell 1- 11 Muellner‐ Riehl et al. (2019) Slope standard deviation

raster 10,000 Mollweide WGS84 Global* Standard deviation of slope in a grid cell

0-22 Muellner‐ Riehl et al. (2019)

Slope range raster 10,000 Mollweide WGS84 Global* Slope range in a grid cell

0-78.9 Muellner‐ Riehl et al. (2019)

River Length raster 10,000 Mollweide WGS84 180 W - 180 O 90 S - 56.6 N* Total length (km) of all rivers in a grid cell 0-103 Muellner‐ Riehl et al. (2019)

Lake Area raster 10,000 Mollweide WGS84 Global* Total lake area (km2 ) in a grid cell 0-100 Muellner‐ Riehl et al. (2019) Geodiversity Sum Raster 10,000 Mollweide WGS84 180 W - 180 O 90 S - 56.6 N

Sum of the sub indices

2-20 Muellner‐ Riehl et al. (2019) In person

Soil grids raster 250 Not projected

WGS84 Global* Soil type code 118 unique types

Hengl et al. (2017)

link

GLiM vector - World Eckert IV

WGS84 Global Lithology types (first level and third level classification) 437 unique types Hartmann & Moosdorf (2012) link Human Influence Human Footprint 1993 & 2009

raster 1000 Mollweide WGS84 Global* Human footprint score 0-50 Venter et al. (2016a), 2018 release link Human Modification 2016

raster 1000 Mollweide WGS84 Global* Human

modification score

0-1 Kennedy et al. (2019)

*excluding Antarctica

2.2 Geoprocessing

2.2.1 Sampling in geoparks and random samples

To compare geoparks with the world, Asia and Europe, random samples were taken from the GGI (sub) layers, human influence layers and HFchange layer. For one random sample set, the same amount of random sample areas as geoparks used in the analysis were created, and this procedure was iterated 100 a times (table 2). In the GGI dataset, no RiverLength data was available above 56.6⁰

(9)

9 Northern latitude and subsequently no GDsum data. The analysis of the RiverLength and GDsum layers was therefore constrained to 56.6 deg N, and did not include the geoparks above this latitude. The other GGI sub layers and human influence data have global coverage and included all geoparks (table 2). For the GGI correlation analysis all random samples were constrained to 56.6 degrees northern latitude, since all layers were compared to the GDsum layer. Similar procedures were followed for the Asian and European analysis.

Table 2. The number of geoparks and random samples involved in each of the analyses.

Analysis Extent Ngeoparks Nrandom

Geodiversity:

- SlopeSTD, SlopeRange, LakeArea SoilDiv & LithoDiv comparison.

Human Influence:

- HF change calculation

- HF1993, HF2009, HF change & gHM comparison

Global 147 14,700

Asia 61 6100

Europe 74 7400

Geodiversity:

- GDsum & RiverLength comparison. - Correlation analysis Global 56.6⁰N 139* 13,900 Asia 56.6⁰N 61 6100 Europe 56.6⁰N 66* 6600

* Geoparks excluded from analysis: Shetland (UK), North West Highlands (UK), Magma (Norway), Gea Norvegica (Norway), Trollfjell (Norway), Rokua (Finland), Katla (Iceland), Reykjanes (Iceland).

The locations of the random samples were created using the create random points tool, and the area where the random points could be placed was constrained to the terrestrial surface. The frequency distribution of the random sample areas resembled the frequency distribution of the geopark areas in order to increase the comparability of the data. From a table with radiuses corresponding to the surface area of the included geoparks, a radius was randomly assigned to a random point by performing a join between the two attribute tables. Each point was buffered using its randomly assigned radius, creating a polygon feature class of random sample areas. In the case of polygons overlapping with the sea surface, the overlapping part was removed using the intersect tool, since the data and analysis are constrained to the land surface. The deletion of areas overlapping with the sea surface causes some of the random sample areas to be slightly smaller than the geoparks. For each geodiversity and human influence layer the mean value in the random samples and geoparks was calculated using the Zonal statistics to table tool and stored in geodatabase tables (Appendix III).

2.2.2 Sampling soil and lithology types

For each continent a soil raster was extracted from the SoilGrids (Hengl et al., 2017) raster, using the Extract by Mask tool. These raster layers were projected to the Mollweide projection and converted to polygons. From these polygons the area of each soil type was calculated using Add Geometry Attributes. For the geoparks, the soil grid data was clipped on the geopark polygons and a similar procedure was followed. A spatial join was performed to assign each soil type polygon in geoparks to

(10)

10 the right continent. For lithology types a similar procedure was followed using GliM dataset (Hartmann & Moosdorf, 2012). Since this was vector data, the Clip tool was used instead of Extract by Mask.

2. 3 Statistical Analysis

2.3.1 Comparing geodiversity and human influence between geoparks and random

samples

The GGI, human influence and HFchange scores were compared between geoparks and random samples in order to test if these are significantly higher or different in geoparks (fig. 1, H1, H3 & H4). The normality of the data was examined visually by making histograms and by performing a Lilliefors test. Since the majority of the data was not normally distributed, the non-parametric Mann-Whitney U-test was used to compare the random samples and geoparks. The geopark means were compared to each of the random sample sets, so 100 comparisons were made per layer. For the GGI layers a one sided Mann-Whitney U test was used with:

H0 = the mean geodiversity is the same in geoparks and random samples

Ha = the mean geodiversity is larger in geoparks than in random samples.

For the human influence and HFchange layers a two-sided Mann-Whitney U-test was used, with: H0 = human influence is the same in geoparks and random samples

Ha = human influence is different in geoparks than in random samples.

Since for each layer the geopark values were compared to 100 random sample sets, the p values needed to be corrected for the family-wise error rate (FWER), caused by multiple testing. The FWER is the probability of making at least one type I error and is 99.4% when performing 100 tests (Ranganathan et al., 2016). The Holm-Sidak correction was applied to the p values for both α = 0.05 and alpha = 0.01. This method performs a step-down correction on the p values and has a higher power than other FWER correction methods, such the Bonferroni or Sidak single step methods (Stevens et al., 2017).

2.3.2 Geodiversity index correlation analysis

In order to quantify the contribution of each GGI sub layer to the GDsum (fig 1, H1), spearman's R (rs) was calculated between the GDsum layer and each of the GGI sub layers for both

the random sample sets and geoparks. rs was chosen because it does not assume normality of the

data. Random sample sets that contained null values were removed from the analysis, because the null values changed the order of sample areas, making the results no longer paired observations. For each layer, this resulted in one rs for the geoparks and between 73 and 100 rs for the random

samples. A two-sided test with p<0.025 was used to assess if the geopark rs could be drawn from the

(11)

11

2.3.3 Analysis of distinct soil and lithology types

To assess how well soil and lithology types are represented within geoparks, the total area of each soil and lithology type was calculated on the global, Asian and European extent by summing the surface area of all polygon of a specific type. This calculation was also used for the geopark data. Eventually, this allowed calculating the percentage area of each soil and lithology type that is located in geoparks using:

With

- = the percentage area of a soil/lithology type located in geoparks

- = the total surface area of a soil/lithology type located in geoparks in km2

- = the total surface area of a soil/lithology type in km2

2.4 Software

The pre-processing and a part of the geoprocessing were done both manually and with the model builder in ArcGIS Pro 2.4.2 (Esri, 2019). The random sampling model for geodiversity and human influence layers was written as Python 3.6.8 (Van Rossum & Drake., 2009) scripts, using the ArcPy module (Esri, 2019). Data analysis and visualisation were done in Python as well, using the pandas (McKinney, 2010), SciPy (Virtanen et al., 2020) and Statsmodels (Seabold & Perktold., 2010) packages for the (statistical) analysis and the Matplotlib (Hunter, 2007) and seaborn (Waskom et al., 2020) packages for visualisation. Python scripts for both the geoprocessing and the statistical analysis are available as digital Appendix III via https://github.com/emmamarianina/Geopark-GI-HI.

(12)

12

3. Results

3.1 Geodiversity index in geoparks and random samples

3.1.1 Geodiversity index

The GDsum is significantly higher in geoparks for all comparisons globally and in Europe, but not for comparisons in Asia (fig. 2A; table 3). The LithoDiv, SlopeRange and SlopeSTD are all significantly higher in geoparks than in the random samples, on all extents and on the 99% confidence level (table 3). The median SlopeSTD lies between 1.3-2.6⁰ in random samples, but between 5.1 -6.9⁰ in geoparks, on all extents (fig 2C). This substantial difference is also present in SlopeRange medians (9.7 - 18.4⁰ in random samples, 30.0 - 39⁰ in geoparks; Appendix I fig. 2A). Median LithoDiv ranges from 1.7 - 1.9 in random samples and from 2.5-2.7 types/100 km2 in geoparks (fig. 2D).

In contrast, the SoilDiv (fig. 2B) in geoparks is not significantly higher than in most random sample sets. Medians range from 3.7 - 4.1 in random samples and from 4.1 - 4.6 types/100 km2 in geoparks, on all extents. The hydrological diversity is also not significantly higher in geoparks. The median RiverLength is about 3km lower in geoparks than in random samples (fig. 2E) and the difference in LakeArea not consistent between extents and varies from 0.031 (Asia) to 0.29 km2 (Europe) in random samples and from 0.014 (Asia) to 0.17 km2 (Europe) in geoparks (Appendix I fig 2B).

The random samples have in general a larger range and more extreme maximum values than the geoparks (fig 2). This indicates that geoparks capture only a part of the total range of diversity values that are present on earth. For example, the highest mean GDsum in geoparks lays around 14.5, while the highest random sample values are about 17.

Table 3. Percentage of tests where geodiversity is significantly higher in geoparks than in random sample sets at α= 0.05 and α= 0.01. p values are corrected with the Holm-Sidak correction.

Global Asia Europe

α = 0.05 α = 0.01 α = 0.05 α = 0.01 α = 0.05 α = 0.01 GDsum 100 100 77 34 100 100 SoilDiv 2 0 0 0 3 1 LithoDiv 100 100 100 100 100 100 SlopeRange 100 100 100 100 100 100 SlopeSTD 100 100 100 100 100 100 RiverLength 0 0 0 0 0 0 LakeArea 3 3 1 0 0 0

(13)

13 Figure 2. Distribution of GD sum means (A) and geodiversity sub layer means (B - E) in the geoparks and random samples on the global (blue), Asian (red) and European (green) extent. The random sample box plots display the combined data from all 100 random sample sets.

3.1.2 Geodiversity index correlations

In geoparks, the GDsum is most strongly correlated with LithoDiv (rs = 0.67, 0.77, 0.63 globally, Asia and Europe, respectively). These geopark rs fall within the 95% CI of the random sample rs

distributions, and are thus not significantly different from the random sample rs (fig. 3A). SlopeSTD

and SlopeRange are also strongly correlated to the GDsum in the random samples, but moderately correlated in geoparks. While the random sample rs distribution ranges from 0.61-88, the geopark rs

(14)

14 each extent, the difference between SlopeSTD and SlopeRange correlations is smaller than 0.03 (fig. 3D; Appendix I fig. 3A). All rs for LithoDiv, SlopeSTD and SlopeRange are significant at the 95%

confidence level (Appendix I table 1).

Both RiverLength and LakeArea are weakly correlated to the GDsum and the majority of the correlations are not significant (fig. 3C; Appendix I fig. 3B, table 1). The rs of both layers is weakly

positive for geoparks, but the random sample rs are distributed around zero. All geopark correlations

fall within the 95% CI, except for the LakeArea rs in Europe. The GDsum-SoilDiv correlation in

geoparks is higher than the RiverLength and LakeArea correlation, but still weakly positive (rs =0.32,

0.21 and 0.31 global, Asia and Europe respectively) (fig 3B). The geopark rs is lower than the 95% CI

globally and in Europe, while the Asian geopark rs falls on the lower boundary of the 95% CI (fig 3).

Most correlations are significant, which the exception of the Asian geopark correlation (Appendix I table 1).

Figure 3. Blue histograms showing the distributions of Spearman's r (rs) between the GDsum and geodiversity

sub layers (ABCD columns) for the random sample sets for the global (first row), Asian (second row) and European (third row) extent. The shaded area in medium grey shows the 95% confidence interval calculated from the correlation distributions. The geopark rs are displayed as red dotted lines. Geopark and random sample

correlations are considered significantly different from each other when the red dotted line lies outside of the grey shaded 95% confidence interval. The number of random sample sets used is shown in the upper right corner of the plots.

3.2 Soil and lithology types

All of the 30 main WRB soil types are located within the 147 geoparks. On the global extent, the percentage area in geoparks is <1.2% for all soil types. Andosols have the highest percentage area in geoparks, followed by alisols. Only a very small percentage area (<0.01%) of the histosols, plinthisols, solonetz and albeluvisols is located in geoparks (fig. 4A). On the Asian and European extent, the

(15)

15 percentage area in geoparks are slightly higher in Asia and much higher in Europe (fig 4B). Especially andosols are well represented in Europe with over 35% of their area present in geoparks. However, the percentage area in geoparks is still <5% for most soil types on the Asian and European extent. In particular, durisols, solonetz, umbrisols and stagnosols are absent in the Asian geoparks. Ferrasols, gypsisols, planosols, nitisols and solonetz do not occur in European geoparks. Only 7 out of the 118 soil types are represented for >1% in geoparks globally. 22% of these soil types are not present in any geopark on the global extent. In Asia, 44% of 111 unique soil types are not present in geoparks and in Europe this is 30% on a total of 90 unique types.

Figure 4. The percentage of the area of each WRB group soil type from the SoilGrids dataset that is located within geoparks compared to the total global area (A) and total Asian and European area (B) of this type. For each of the 15 first legend level lithology classes in the GLiM dataset less than 1% area is covered by the 147 existing geoparks (fig. 5A), since the total geopark area is small in comparison to the global land surface area. Intermediate volcanic rocks and intermediate plutonic rocks are covered the most by geoparks, evaporites the least. None of the ice and glacier surface is present in geoparks. These percentages are in general slightly higher on the Asian extent, and substantially higher on the European extent (fig. 5B). Especially igneous and volcanic rocks are best represented in geoparks, with a percentage area in geoparks >5% in Europe, while the percentage area in geoparks is lower for sedimentary and metamorphic rocks. The evaporites are the exception, only a very small percentage is located in Asian geoparks and 0% is present in European geoparks. The three legend levels combined yield a total of 437 unique lithology type combinations. For only 31 of these types the

(16)

16 percentage area in geoparks >1% on the global extent. One lithology type, unconsolidated sediments with mixed grain sized and black shales mentioned (sumxbs), is almost entirely (99.97) present in the geoparks. On the other hand, 65% of the 437 unique third level lithology classes is not present in any geopark. Of the 227 unique lithology types in Asia the percentage not present in geoparks is 65% as well. In Europe 45% of 177 unique lithology types does not occur in geoparks.

Figure 5: The percentage of the area of each first legend level lithology type from the GLiM dataset that is located within geoparks compared to the total global area (A) and total Asian and European area (B) of this type.

3.3 Human influence in geoparks and random samples

The HF1993, HF2009 and gHM are all higher in geoparks than in random samples, on all extents (fig. 6A; Appendix I fig. 4). This difference is significant for all tests on the global and Asian extent (table 4). In Europe the difference between geoparks and random samples is significant for the majority of the human footprint comparisons. For the gHM the difference is not significant in the majority of comparisons.

3.4 Human footprint change in geoparks and random samples

The HFchange is slightly positive for geoparks and random samples globally and in Asia, but slightly negative for Europe (fig. 6B). For both geoparks and random samples the distribution is symmetrical around zero. The difference between geopark and random sample HF change is not significant on any extent (table 4). However, the minima and maxima of the random sample HF change are higher than the geopark range.

(17)

17 Table 4. Percentage of tests where HI is significantly different in geoparks than in random sample sets at α= 0.05 and α = 0.01. p values are corrected with the Holm-Sidak method.

Global Asia Europe

α = 0.05 α = 0.01 α = 0.05 α = 0.01 α = 0.05 α = 0.01

HF1993 100 100 100 100 81 56

HF2009 100 100 100 100 64 43

HFchange 0 0 0 0 0 0

gHM 100 100 100 100 17 4

Figure 6. Distribution of human influence layer means of the geoparks and random samples on the global (blue), Asian (red) and European (green) extent. The random sample box plots display the combined data from all 100 random sample sets.

(18)

18

4. Discussion

This study compared geodiversity, human influence and human influence change in geoparks with geodiversity, human influence and human influence change in the world, Asia and Europe. These results show that the GDsum, topographic diversity and LithoDiv are significantly higher in geoparks, while soil and hydrological diversity are not. The GDsum is also strongly correlated with the LithoDiv and moderately correlated with the topographic layers. In Europe the highest percentages of lithology and soil types are located in geoparks, but not all types are represented. Human influence is significantly higher in geoparks, except for Europe, and the HFchange is not significantly different from the random samples.

4.1 Comparing the geodiversity index

The GDsum, LithoDiv, SlopeSTD and SlopeRange are significantly higher in geoparks. The strong correlation between the GDsum and LithoDiv indicates that lithological diversity highly contributes to the total geodiversity, in both geoparks and random samples. The finding that the correlation between the GDsum and topographic diversity is only moderate in geoparks, while strong relations occur in random samples, emphasizes the importance of lithological diversity for the total geodiversity in geoparks even more. This is in line with earlier listings of the main geoheritage types in geoparks, which are mostly related to geology and geomorphology (e.g. karst, tectonics, landforms, and stratigraphical geoheritage) (Brilha, 2018). Soil diversity is not significantly higher in geoparks and the correlation between the GDsum and SoilDiv is lower in geoparks than in random samples. This demonstrates that soil diversity is less important in geoparks and is also in agreement with the ill representation of pedological heritage in geoparks reported by Ruban (2017). Though a literature review concluded that hydrological diversity gets significant attention in geoparks (Ruban, 2019), this is not reflected by the lower hydrological diversity in geoparks and weak correlation between the hydrological diversity layers and the GDsum.

A possible explanation for the high lithological and topographical diversity is the emphasis on lithological and geomorphological diversity in the geopark application procedure, which was mentioned in the hypothesis. Aspiring parks get more points assigned for higher numbers of different rock types and distinct geological or geomorphological features (UNESCO, 2016). Hydrology and soils are not explicitly taken into account in this assessment. The higher topographic variability in geoparks is likely also related to high number of geoparks located in mountains. Of all geoparks, 58.5% are located within mountainous areas as defined by the Global Mountain Biodiversity Assessment (GMBA), while only 12.3% of the global land area is considered mountainous (Körner et al., 2017). The areas with the highest GDsum values are also mostly located in mountainous areas (fig. 7). In addition, the many features aesthetic values and the better exposed rocks and tectonic structures make mountainous areas popular for geopark establishment (Brilha, 2018).

The bias towards mountainous areas also offers an additional explanation for the not significantly higher soil and hydrological diversity. Hydrology in mountainous areas consists often of small and fragmented streams that might be too small to appear in the HydroSHEDS dataset that was used for the RiverLength layer (Lehner et al., 2013). Patterns of high soil diversity do also not necessarily coincide with mountainous areas or areas with high lithology diversity (Appendix I fig 5), as soil diversity is the result of the interplay of more factors than relief and parent material alone

(19)

19 (Jenny, 1961). In earlier research lower soil diversity was associated with extreme temperatures and precipitation (Minasny et al., 2010) and increasing elevation (Vacek et al., 2020). In the humid tropics and subtropics soil diversity was highly influenced by time and erosion processes in mountainous and hill areas, while parent material had more influence on diversity in the lower altitudes (Gracheva, 2011).

While the GDsum, SlopeSTD, SlopeRange and LithoDiv are significantly higher in geoparks, the random samples contain higher maximum values (fig 2). Considering that the highest GDsum values are located in mountainous areas (fig. 7), it could be that these areas are too sparsely populated to become geoparks, as an area needs to have at least 7 inhabitants/km2 to become a geopark (UNESCO, 2020). This means that though geoparks capture regions with high geodiversity, they do not capture the areas with the highest geodiversity on earth.

Figure 7. Map showing the areas with the highest GDsum values, overlain with mountain polygons as defined by the GMBA inventory. Areas with high geodiversity highly coincide with mountainous areas. The highest mean GDsum values in geoparks are between 12-15, while the highest mean GDsum in random samples are between 15 and 17.5.

4.2 Comparing soil and lithology types

The percentage area in geoparks are very low on the global and Asian extent, due to the small geopark surface relative to the continent, but considerable percentages of some soil and lithology types were found in Europe. The high percentages of igneous rocks are in line with an earlier study on representation of geodiversity types in geoparks on the national level (Ruban, 2017) and the very high percentage area in geoparks of andosols, pyroclastics and volcanic rocks suggest that relatively many geoparks are in areas shaped by (past) volcanic activity. This endorses earlier statements about the popularity of volcanic geoheritage in geoparks, as concluded from qualitative analyses (Liu et al., 2012; Brilha, 2018). Many of the lowest percentage area in geoparks soil types, such as solonetz, plinthosols, durisols, gypsisols are zonal soils associated with (semi) arid climatic zones. The same holds true for the low percentage area in geoparks for evaporites, which are only preserved in arid climates. This is a result of the absence of geoparks in arid climates, and stresses the importance of

(20)

20 reducing the spatial bias in geopark locations, which was earlier recommended too by Ibáñez et al. (2018) and Ibáñez & Brevik (2019). Geoparks should be present in all climate zones to capture all soil types, as many soil types are zonal and thus restricted to a specific climatic area (Bockheim, 2005).

Despite relatively high percentage area in geoparks in European geoparks, still a considerable amount of soil and lithology types are not present in geoparks when taking into account the most detailed legend levels of the datasets. This percentage is higher globally and in Asia than in Europe, indicating that in Europe geoparks are a better representation of all soil and lithology types present on the continent. This is most likely due to the large amount of geoparks in Europe which are relatively evenly distributed over the continent when compared to Asia or the world (Appendix I, fig. 1). Some soil and lithology types are globally low in abundance and have a small percentage area in geoparks, which means that very little of the already little that is there is protected. This applies for example to histosols, while these soils can have a high scientific value because these can function as paleo archives (Seijmonsbergen et al., 2010; Seijmonsbergen et al., 2019). In future geopark designations attention should be paid to poorly or non-represented soil and lithology types to increase the representativeness of geoparks.

4.3. Comparing human influence

Human influence is higher in geoparks for all three indices, indicating a high presence of human activities in geoparks. In Europe, the difference between geoparks and random samples was not significant for all comparisons, especially not for the gHM. This is most likely because of areas with low human influence, such as Siberia, the Sahara and central Asia (Venter et al., 2016a; Kennedy et al., 2019) that are randomly sampled relatively often because of their large surface area. As such areas are restricted in Europe the difference between geoparks and random samples becomes less pronounced. As mentioned in paragraph 4.1 is having a significant population a prerequisite for becoming a geopark (UNESCO, 2020). Thereby is geoconservation in geoparks combined with the promotion of the local economy (Farsani et al., 2011; UNESCO 2015) and this clearly leads to geoparks being located in areas modified by human activities. These human activities in geoparks do not necessarily have to be a threat to geodiversity when properly managed and some geoheritage types, such as mining, were even created by human activities (Crofts & Gordon, 2015; Brilha, 2018). However, many activities associated with economical development as urbanization, geotourism and land use change can be potential threats to geodiversity (Crofts & Gordon, 2015; Hjort et al., 2015) and the high human influence found in geoparks further emphasizes the need of adequate management in geoparks.

4.4 Comparing human influence change

In the 1993-2009 period the human footprint decreased in about half of the geoparks and increased in the other half, given the roughly symmetrical distribution around zero. There is no significant difference between geoparks and random samples, indicating that geopark establishment does not lead to a different HFchange trend. It should be noted that UNESCO only begun with establishing geoparks in 2001 and that 87 out of 147 geoparks in the analysis were established after 2009, but as geoparks need to be managed by a recognized legislative body (UNESCO, 2015) it can be assumed that some sort of conservation measures were already implemented. Looking at the spatial distribution of areas with increasing and decreasing trends, the HF decreased the most in

(21)

north-21 western Europe and a part of southeastern China, where many geoparks are situated (Appendix I fig. 6). However, adjacent regions where also many geoparks are situated such as southern Europe, Southeast Asia and central China experienced an increase in HF, which could explain the symmetrical distribution of the HFchange. In Asian geoparks a lacking legislative framework led to geotourism causing an increase in infrastructure, agriculture and urbanization, which eventually lead to erosion of geosites (Sumanapala & Wolf, 2020). Such a study is not available for European geoparks. Decreases in the HF mainly occurred in countries with strong legislative enforcement (Venter et al. 2016b), which further emphasizes the importance of a strong legislative framework and legislative body in geoparks.

4.5 Technical limitations and future directions

This research had to deal with several technical limitations, the most evident being the limited cover of the RiverLength and GDsum layers, causing the northern most geoparks to be excluded from some of the analyses. Due to the incompleteness of the geopark polygon dataset about half of the geopark areas had to be estimated, causing inaccuracies. Thereby, the global datasets used are unavoidably subject to spatially uneven data availability and quality, and previous research has shown that higher diversity in areas can also be the result of the larger scale of the input data in these areas (Minasny et al., 2010; Gerasimova et al., 2020). In addition, due to its global scale the SoilGrids data does not contain some uncommon, but relevant soil types such as paleosoils and anthroposols, while these can contain valuable information about processes that shaped landscapes in the past (Seijmonsbergen et al., 2019). In both the HF and gHM climate change is not included as a pressure, while it is considered to be a threat to geodiversity (Gordon et al., 2012; Crofts & Gordon, 2015). In this research a grid-based GGI was used, but different methods that are qualitative-quantitative based (Zwoliński et al., 2018) or use centroids (Forte et al., 2018; da Silva et al., 2019) are being developed. Topographic diversity was used as a substitute for geomorphological diversity in the GGI, as there is no global geomorphological map available. Hydrological diversity was assumed to increase linearly with increasing river length and lake area, but differences between drainage densities or lake types were not taken into account. Despite these data driven limitations, the data used in this research are amongst the best available datasets that allow an analysis and comparison the global extent. If to increase the accuracy of these results, further refinements and expansion of data coverage are necessary.

These results provide organizations, including UNESCO, as well as policymakers with important information on the extent and representativeness of protected global geodiversity and identify aspects to which more attention should be paid. In order to fully represent global geodiversity, geoparks should focus more on soil and hydrological diversity, in addition to lithological and topographic diversity. Additional research will be needed for this, as soil sites are not sufficiently studied (Brilha & Reynard, 2018). Thereby should future geoparks be targeted at the large number of soil and lithology types that are not yet included in geoparks. This focus will probably also help straightening the uneven distribution of geoparks globally. Establishing monitoring programs in geoparks as proposed by Brilha (2018) and Sumanapala & Wolf (2020) are needed, as these results show that geopark areas are significantly modified by human activities. In addition, the HF and gHM layers only reflected the cumulative human impact on areas. More research into the relative importance of each specific driver of human influence is necessary to make management practices targeted more efficiently.

(22)

22

5. Conclusion

The results show that the geoparks' aim of protecting sites and landscapes with international geological significance leads to the protection of areas with a significantly higher geodiversity, lithological and topographical diversity when compared to global, Asian and European geodiversity. Soil diversity and hydrological diversity, however, are less represented in geoparks and should be focus in future geopark designations. Soil and lithology types that are not or hardly represented in geoparks should be included in future geoparks, paying special attention to those types that are low in abundance and are vulnerable to human impacts. Geoparks are located in areas that are relatively highly modified by human activities and HF change data showed a similar trend in geoparks and random samples. Human activities should be monitored and assessed by a strong legislative body to prevent damage to geodiversity. Despite some shortcomings are the datasets used in this analysis the best available data that allow for a global comparison. Future expansion and refinement of global datasets should further increase the accuracy of these results.

(23)

References

Araujo, A. M., & Pereira, D. Í. (2018). A New Methodological Contribution for the Geodiversity Assessment: Applicability to Ceará State (Brazil). Geoheritage, 10(4), 591–605.

Bockheim, J. G. (2005). Soil endemism and its relation to soil formation theory. Geoderma, 129(3–4), 109–124.

Boothroyd, A., & Henry, M. M. (2019). Old processes, new movements: The inclusion of geodiversity in biological and ecological discourse. Diversity, 11(11).

Brilha, J. (2016). Inventory and Quantitative Assessment of Geosites and Geodiversity Sites: a Review. Geoheritage, 8(2), 119–134.

Brilha, J. (2018). Geoheritage and geoparks. In J. Brilha & E. Reynard (Eds.), Geoheritage: Assessment, Protection, and Management (pp. 323–335). Chennai: Elsevier Inc.

Brilha, J., & Reynard, E. (2018). Geoheritage and geoconservation: The challenges. Geoheritage: Assessment, Protection, and Management. Elsevier Inc.

Cañadas, S., & Ruiz-Flaño, E. (2007). Geodiversity : Concept, Assessment and Territorial Application. Boletin de La A.G.E., (45), 389–394.

Crofts, R., & Gordon, J. E. (2015). Geoconservation in protected areas. In S. F. and I. P. G. L. Worboys, M. Lockwood, A. Kothari (Ed.), Protected Area Governance and Management (pp. 531–568). Canberra: ANU Press.

da Silva, M. L. N., do Nascimento, M. A. L., & Mansur, K. L. (2019). Quantitative Assessments of Geodiversity in the Area of the Seridó Geopark Project, Northeast Brazil: Grid and Centroid Analysis. Geoheritage, 11(3), 1177–1186.

dos Santos, F. M., de La Corte Bacci, D., Saad, A. R., & da Silva Ferreira, A. T. (2020, April 3). Geodiversity index weighted by multivariate statistical analysis. Applied Geomatics. Springer. Esri (2019). ArcGIS Pro: Release 2.4. Redlands, CA: Environmental System Research Institute. Farsani, N. T., Coelho, C., & Costa, C. (2011). Geotourism and geoparks as novel strategies for

socio-economic development in rural areas. International Journal of Tourism Research, 13(1), 68–81. Fernández, A., Fernández, T., Pereira, D. I., & Nieto, L. M. (2020). Assessment of Geodiversity in the

Southern Part of the Central Iberian Zone (Jaén Province): Usefulness for Delimiting and Managing Natural Protected Areas. Geoheritage, 12(1), 1–16.

Forte, J. P., Brilha, J., Pereira, D. I., & Nolasco, M. (2018). Kernel Density Applied to the Quantitative Assessment of Geodiversity. Geoheritage, 10(2), 205–217.

Gerasimova, I. I., Golovleva, I., Konyushkova, M. V., Sorokin, A. S., & Krasilnikov, P. V. (2020). Assessment of soil diversity using soil maps with different scales in Eastern Fennoscandia, Russia. Geoderma Regional, 21, e00274.

Global Geopark Network (2019). Member list. Accessed May, 2019 from

http://www.globalgeopark.org/homepageaux/tupai/6513.htm

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.

Gordon, J. E., Barron, H. F., Hansom, J. D., & Thomas, M. F. (2012). Engaging with geodiversity-why it matters. Proceedings of the Geologists’ Association, 123(1), 1–6.

Gracheva, R. (2011). Formation of soil diversity in the mountainous tropics and subtropics: Rocks, time, and erosion. Geomorphology, 135(3–4), 224–231.

Gray, M. (2008). Geodiversity: developing the paradigm. Proceedings of the Geologists’ Association, 119(3–4), 287–298.

Gray M (2013) Geodiversity: valuing and conserving abiotic nature, 2nd edn. Wiley, Chichester Hartmann, J., & Moosdorf, N. (2012). The new global lithological map database GLiM: A

(24)

24 13(12).

Hengl, T., De Jesus, J. M., Heuvelink, G. B. M., Gonzalez, M. R., Kilibarda, M., Blagotić, A., … Kempen, B. (2017). SoilGrids250m: Global gridded soil information based on machine learning. PLoS ONE, 12(2), e0169748.

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

Hjort, J., & Luoto, M. (2010). Geodiversity of high-latitude landscapes in northern Finland. Geomorphology, 115(1–2), 109–116.

Hunter, J.D. (2007). Matplotlib: A 2D Graphics Environment. Computing in Science & Engineering, 9, 90-95.

Ibáñez, J. J., & Brevik, E. C. (2019, November 1). Divergence in natural diversity studies: The need to standardize methods and goals. Catena. Elsevier B.V.

Ibáñez, J. J., Brevik, E. C., & Cerdà, A. (2018). Geodiversity and geoheritage: Detecting scientific and geographic biases and gaps through a bibliometric study. Science of the Total Environment, 659, 1032–1044.

IUSS Working Group WRB (2006). World reference base for soil resources 2006. 2nd edition. World Soil Resources Reports No. 103. FAO, Rome.

Jenny, H. (1961). Derivation of State Factor Equations of Soils and Ecosystems. Soil Science Society of America Journal, 25(5), 385–388.

Kennedy, C. M., Oakleaf, J. R., Theobald, D. M., Baruch-Mordo, S., & Kiesecker, J. (2019). Managing the middle: A shift in conservation priorities based on the global human modification gradient. Global Change Biology, 25(3), 811–826.

Lehner, B., Verdin, K., & Jarvis, A. (2008). New global hydrography derived from spaceborne elevation data. Eos, Transactions American Geophysical Union, 89(10), 93-94.

Lehner, B., Verdin, K., & Jarvis, A. (2013). HydroSHEDS Technical Documentation Version 1.2. EOS Transactions, 89(10), 26.

Liu, Jiaqi, Liu, Jiali, Chen, X., & Guo, W. (2012). Volcanic Natural Resources and Volcanic Landscape Protection: An Overview. In Updates in Volcanology - New Advances in Understanding Volcanic Systems. InTech.

McKinney (2010). Data structures for statistical computing in Python. Proceedings of the 9th Python in Science Conference, 51-56.

Messager, M. L., Lehner, B., Grill, G., Nedeva, I., & Schmitt, O. (2016). Estimating the volume and age of water stored in global lakes using a geo-statistical approach. Nature communications, 7, 13603.

Minasny, B., McBratney, A. B., & Hartemink, A. E. (2010). Global pedodiversity, taxonomic distance, and the World Reference Base. Geoderma, 155(3–4), 132–139.

Muellner‐Riehl, A. N., Schnitzler, J., Kissling, W. D., Mosbrugger, V., Rijsdijk, K. F., Seijmonsbergen, A. C., … Favre, A. (2019). Origins of global mountain plant biodiversity: Testing the

‘mountain‐geobiodiversity hypothesis.’ Journal of Biogeography, 46(12), 2826–2838. Pereira, D. I., Pereira, P., Brilha, J., & Santos, L. (2013). Geodiversity assessment of Paraná State

(Brazil): An innovative approach. Environmental Management, 52(3), 541–552.

Ranganathan, P., Pramesh, C., & Buyse, M. (2016). Common pitfalls in statistical analysis: The perils of multiple testing. Perspectives in Clinical Research, 7(2), 106.

Ruban, D. A. (2017). Geodiversity as a precious national resource: A note on the role of geoparks. Resources Policy, 53, 103–108.

Ruban, D. A. (2019). Water in descriptions of global geoparks: Not less important than geology? Water (Switzerland), 11(9), 1–12.

Ruban, D. A., & Yashalova, N. N. (2018). Geodiversity meanings in global geoparks: an empirical study. Environmental Earth Sciences, 77(22), 771.

(25)

25 Human Footprint and the Last of the Wild. BioScience, 52(10), 891–904.

Santos, D. S., Mansur, K. L., Gonçalves, J. B., Arruda, E. R., & Manosso, F. C. (2017). Quantitative assessment of geodiversity and urban growth impacts in Armação dos Búzios, Rio de Janeiro, Brazil. Applied Geography, 85, 184–195.

Seabold, S. & Perktold, J. (2010). Statsmodels: Econometric and statistical modeling with python. Proceedings of the 9th Python in Science Conference. 92-96.

Seijmonsbergen, A. C., Sevink, J., Cammeraat, L. H., & Recharte, J. (2010). A potential

geoconservation map of the Las Lagunas area, northern Peru. Environmental Conservation, 37(2), 107–115.

Seijmonsbergen, A. C., De Jong, M. G., Hagendoorn, B., Oostermeijer, J. G., Rijsdijk, K. (2018) Geodiversity Mapping in Alpine Regions In Mountains, Climate and Biodiversity (pp 155-170). Chichester: Wiley Blackwell.

Seijmonsbergen, A. C., van den Ancker, J. A. M., Jungerius, P. D., & Norder, S. J. (2019). Can

geodiversity help to save the soil archives? In J. van Mourik & J. Van der Meer (Eds.), Reading the Soil Archives (pp. 275–298). Elsevier.

Shui, W., & Xu, G. (2016). Analysis of the influential factors for changes to land use in China’s Xingwen Global Geopark against a tourism development background. Geocarto International, 31(1), 22–41.

Stevens, J. R., Al Masud, A., & Suyundikov, A. (2017). A comparison of multiple testing adjustment methods with block-correlation positively-dependent tests. PLOS ONE, 12(4), e0176124. Sumanapala, D., & Wolf, I. D. (2020, September 1). Man-Made Impacts on Emerging Geoparks in the

Asian Region. Geoheritage. Springer.

UNESCO (2015). Statutes of the International Geoscience and Geoparks Programme. Accessed 8-3-2019 from https://unesdoc.Unesco.org/ark:/48223/pf0000260675

UNESCO (2016). Evaluation Document A - Self Evaluation. Downloaded 1-4-2019, from

http://www.Unesco.org/new/en/natural-sciences/environment/earth-sciences/Unesco-global-geoparks/application-process/

Unesco (2019). List of geopark members. Accessed May, 2019 from

http://www.unesco.org/new/en/natural-sciences/environment/earth-sciences/unesco-global-geoparks/list-of-unesco-global-geoparks/

UNESCO (2020). Self-Evaluation Checklist, Explanatory Notes for aspiring UNESCO Global Geoparks. Downloaded 29-7-2020 from http://www.unesco.org/new/en/natural-sciences/environment/earth-sciences/unesco-global-geoparks/application-process/

Vacek, O., Vašát, R., & Borůvka, L. (2020). Quantifying the pedodiversity-elevation relations. Geoderma, 373(April), 1–8.

Van Rossum, G. & Drake, F. L. (2009). Python 3 Reference Manual. Scotts Valley, CA: CreateSpace. Venter, O., Sanderson, E. W., & Magrach, A. (2016a). Global terrestrial Human Footprint maps for

1993 and 2009. Scientific Data, 3, 160067.

Venter, O., Sanderson, E. W., Magrach, A., Allan, J. R., Beher, J., Jones, K. R., … Watson, J. E. M. (2016b). Sixteen years of change in the global terrestrial human footprint and implications for biodiversity conservation. Nature Communications, 7(1), 12558.

Virtanen, P., Gommers, R., Oliphant, T.E., Haberland, M., Reddy., T. et al. (2020). SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat Methods 17, 261–272. Wang, Lulin, Tian, M., & Wang, Lei. (2015). Geodiversity, geoconservation and geotourism in Hong

Kong Global Geopark of China. Proceedings of the Geologists’ Association, 126(3), 426–437. Waskom, M., Botvinnik, O., Ostblom, J., Gelbart, M., Lukauskas, S. et al. (2020). mwaskom/seaborn:

v0.10.1 (April 2020) (Version v0.10.1). Zenodo.

Yamazaki, D., Ikeshima, D., Tawatari, R., Yamaguchi, T., O'Loughlin, F. et al. (2017). A high accuracy map of global terrain elevations. Geophysical Research Letters. 44(11), 5845-5853.

(26)

26 Zwoliński, Z., Najwer, A., & Giardino, M. (2018). Methods for Assessing Geodiversity. In E. Reynard &

(27)

27

Appendix I

Additional figures and tables

Figure 1. Distribution of the 147 geoparks. Most geoparks are located in Europe and East Asia, especially in China.

Figure 2. Distribution of SlopeRange (A) and LakeArea means (B ) in the geoparks and random samples on the global (blue), Asian (red) and European (green) extent. The random sample box plots display the combined data from all 100 random sample sets.

(28)

28 Figure 3. Blue histograms showing the distributions of Spearman's r (rs) between the GDsum and geodiversity

sub layers (AB columns) for the random sample sets for the global (first row), Asian (second row) and European (third row) extent. The shaded area shows the 95% confidence interval calculated from the correlation distributions. The geopark rs are displayed as red dotted lines. Geopark and random sample correlations are

considered significantly different from each other when the red dotted line lies outside of the green shaded 95% confidence interval. The number of random sample sets used is shown in the upper right corner of the plots.

Table 1. Correlations between geodiversity sublayers and the GDsum. For the random samples the 95% CI of the rs distribution is given. The percentage of random sample rsthat are significant at alpha = 0.05 is displayed

below the interval.

Geoparks rs Random samples rs [min, max]

Global Asia Europe Global Asia Europe

SoilDiv 0.32 0.21 0.31 [0.46 - 0.70] [0.21 - 0.68] [0.39 - 0.77] p<0.01 n.s. p<0.05 100% 91.4% 100% LithoDiv 0.67 0.77 0.63 [0.62 - 0.78] [0.53 - 0.81] 0.53 - 0.82] p<0.01 p<0.01 p<0.01 100% 100% 100% SlopeSTD 0.47 0.40 0.44 [0.61 - 0.79] [0.67 - 0.88] [0.63 - 0.85] p<0.01 p<0.01 p<0.01 100% 100% 100% SlopeRange 0.46 0.38 0.45 [0.63 - 0.80] [0.67 - 0.88] [0.65 - 0.86] p<0.01 p<0.01 p<0.01 100% 100% 100% RiverLength 0.19 0.23 0.12 [-0.01 - 0.34] [-0.20 - 0.28] [-0.31 - 0.13] p<0.05 n.s. n.s. 55.1% 1% 10.3% LakeArea 0.09 0.11 0.08 [-0.09 - 0.26] [-0.35 - 0.18] [-0.39 - 0.08] n.s. n.s. n.s. 18.2% 7% 25.3%

(29)

29 Figure 4. Distribution of HF1993 (A) and gHM 2016 means (B ) in the geoparks and random samples on the global (blue), Asian (red) and European (green) extent. The random sample box plots display the data from all 100 random sample sets together.

(30)

30 Figure 5. Maps showing the different spatial patterns of soil (A) and lithological (B) diversity. Mountainous areas are marked with black lines. Based on data by Hannes Versteegh from Muellner-Riehl et al. (2019).

(31)

31 Figure 6. Maps showing the areas where the human footprint decreased (A) and increased (B) between 1993 and 2009. Based on data from Venter et al. (2016a)

(32)

32

Appendix II

Geopark database

continent country geopark name

year of design ation Shapfile surface area (km^2)

Africa Morocco mgoun_global 2014 NA 5730

Africa Tanzania ngorongoro_lengai 2018 NA 11886

Asia China hongkong 2011 digitised from online map 150

Asia China leiqiong 2006 digitised from online map 3050

Asia China ningde 2010 digitised from online map 2660.34

Asia China yandangshan 2005 digitised from online map 298.8

Asia China alxa_desert 2009 NA 630.37

Asia China arxan 2017 NA 3653.21

Asia China dali_mount_cangshan 2014 NA 933

Asia China danxiashan 2004 NA 292

Asia China dunhuang 2015 NA 2067

Asia China fangshan 2006 NA 954

Asia China funiushan 2006 NA 5858.52

Asia China guangwushan_nuoshuihe 2018 NA 1818

Asia China hexigten 2005 NA 1750

Asia China huanggang_dabieshan 2018 NA 2625.54

Asia China huangshan 2004 NA 1200

Asia China jingpohu 2006 NA 1400

Asia China jiuhuashan 2019 NA 139.7

Asia China keketuohai 2017 NA 2337.9

Asia China leye_fengshan 2010 NA 930

Asia China longhushan 2007 NA 966.63

Asia China lushan 2004 NA 500

Asia China mount_kunlun 2014 NA 7033

Asia China sanqingshan 2012 NA 433

Asia China shennongjia 2013 NA 1022.72

Asia China songshan 2004 NA 464

Asia China stone_forest 2004 NA 350

Asia China taining 2005 NA 492.5

Asia China taishan 2006 NA 418.36

Asia China tianzhushan 2011 NA 413.14

Asia China wangwushan_daimeishan 2006 NA 986

Asia China wudalianchi 2004 NA 720

Asia China xingwen 2005 NA 156

Asia China yanqing 2013 NA 620.38

Asia China yimengshan 2019 NA 1804.76

A quantitative assessment of geodiversity and human influence in geoparks