Population and Landslide Risk Evolution in Long Time Series. Case Study of the Valencian Community (1920–2021)

1. Introduction

In recent years, society has seen an increase in damage caused by natural hazards. Floods stand out in particular for their visibility, extent, intensity and frequency. All signs indicate that climate change is responsible for an increased incidence of these events and their severity is accentuated by a lack of adequate planning. In the case of landslides, with unequal influence but marked by unprecedented meteorological events (Chen et al., 2024; Gariano & Guzzetti, 2016; Haque et al., 2019), one important factor is the increased exposure to risk of susceptible elements: population, buildings, infrastructure (Moore & McInnes, 2016).

The assessment, management, and mitigation of this type of risk has been extensively discussed in the literature over the last few decades. A multitude of papers on these issues can be cited, the most important being the synthesis proposed by Dai et al., 2002 with a critical review of the landslide research and the strategies for reducing losses, as well as relevant publications by Lee & Jones, 2004 and Bell et al., 2006 with a multidisciplinary perspective on landslide management. Corominas et al.’s comprehensive review of quantitative risk assessment methods should also be mentioned (2014).

One of the most important aspects of these types of events is the prior assessment of the potential damage that may be caused taking into consideration the vulnerability of the elements at risk. Different dimensions of this vulnerability are often taken into account and there can be of many types of them, although the social, personal and structural are those generally addressed (Fuchs, 2009; Kienberger et al., 2009) and linked together (Fuchs, 2009; Kappes et al., 2012; Papathoma-Köhle et al., 2011). Integrated approaches to evaluation are increasingly being addressed as one would expect (Fuchs et al., 2011; Karagiorgos et al., 2016) although they require a highly detailed quantitative evaluation. Indeed, the variety of elements that are subject to damage and loss together and their different characteristics (e.g., buildings, public works, people) require complex analysis, so such studies are not abundant (Keiler, 2004; Michael-Leiba et al., 2003; Promper & Glade, 2016).

It is common for only two areas of work to be taken into account, namely physical and social. Attempts have been made to establish empirical relationships with regard to property damage between the intensity of the process and the type and quantity of exposed elements to value the estimated loss, usually in monetary units (e.g., Fuchs et al., 2015; Guillard-Goncalves et al., 2016; Kappes et al., 2012; Petrucci & Gullà, 2010; Promper et al., 2015; Silva & Pereira, 2014; Uzielli et al., 2015b, 2015a; Winter et al., 2014).

In the other area of work, on social damage, research has focused on how affected societies adapt to these events, and how to cope with their consequences (e.g., De Oliveira Mendes, 2009; Di Napoli et al., 2023; Wood et al., 2010). The study by (Bukhari et al., 2023) on long-term social impacts in mountain communities is of note. Work has also been done to assess the relationship between the occurrence of a process and the damage caused to people by defining curves between the two variables (e.g., Guzzetti, 2000; Vranken et al., 2013). Human activity itself, by causing changes in soil characteristics, is clearly a trigger for these processes in the most inhabited areas, leading to a greater impact on the population (Fidan et al., 2024).

All of these studies are generally based on historical data on risky events that affected the population and collected in specific inventories (e.g., Gariano et al., 2015; Guzzetti et al., 2005). Logically, these historical inventories are often insufficient and incomplete, which means that the probabilistic calculation has often been based on expert analysis (Michael-Leiba et al., 2003). Consequently, very different approaches need to be adequately addressed to tackle these complex problems (Fuchs et al., 2011). A significant number of data types with different ranges and units (material, personal, social) must also be taken into account to estimate the value of the damage caused if quantitative risk analysis is used (Corominas et al., 2014; Remondo et al., 2008; Schwendtner et al., 2013; Zêzere et al., 2007, 2008).

An accurate assessment of the number of people present at a given time and place is a necessary first step, as it is crucial to develop a consistent emergency plan. Knowing as precisely as possible the location of the potentially exposed population is essential to ensure emergency plan efficiency and to reduce associated costs related to rescue and community recovery (e.g., Aubrecht et al., 2013; Freire & Aubrecht, 2012; Opdyke & Fatima, 2024; Su et al., 2010). According to Bhaduri et al. (2002), locating the population at risk should be the first step to mitigating personal damage and reducing casualties.

Their results must be accompanied by an analysis of the local situation, developing a system for assessing the available social, economic, and health care fabric Wen et al. (2024). Urban resilience is generally particularly dependent on emergency response along with its potential for immediate recovery and reconstruction after any catastrophic event (Meerow et al., 2016). Civil society must be able to prepare, anticipate, preserve, absorb, respond, resist, recover, mitigate, learn, and adapt to these situations (Tian & Lan, 2023).

Population distribution data for large areas and in broad historical series are essential first to learn how the risk has evolved, the places where it is most intensely concentrated and also to provide adequate criteria to optimise territorial management and prevent damage (Haque et al., 2019). Significant differences in population density are expected, especially in larger study areas where there may be different types of use (urban concentration or rural or holiday dispersion). This study’s goal is to analyse the historical series of geolocalised population affected by these events to establish variation patterns and define areas of preferential application of territorial policies that minimise possible damage caused by landslides.

The scope of the study is a large swathe of Spain, the Valencian Community. This region covers 22,000 km² and has a population of five million, is located on the east of the Iberian Peninsula and bordered by the Mediterranean Sea. It is a suitable example, offering a complex terrain and different types of urban layouts, both traditional and holiday. The following sections explain the procedure followed to determine possible causes and variation patterns in a major part of Spain that is exposed to these dangerous events and that can be generalised to other areas.

2. General Methodology

2.1. Objectives

The aim of this paper is to define the preferred areas of action to prevent social damage, to the population in particular, caused by landslides. Finding variation patterns in the historical series and their possible causes is one of the possible ways to achieve the goals set. This requires a sufficiently accurate geo-referenced population, which must be complemented by other data including detailed demographic information.

There are a plethora of reasons why collecting this type of data is arduous, starting with the fact that public bodies are only able to provide the most recent figures, and not always at the appropriate resolution. Indeed, it is rare to find detailed information provided by national demographic agencies that goes beyond a regularised 1 km grid size for large areas and only for recent time periods.

There are studies that, based on dasymetric techniques applied to data provided by national land registries, offer a satisfactory resolution (Garcia et al., 2016), although this study is applied to a reduced scope and without a time frame. Fortunately, in Spain there is a very recent paper that offers a population grid from the beginning of the last century to the present day, with a resolution of 100x100m (Goerlich, 2025). Although it is not the most appropriate resolution for events that are not characterised by their enormity, it is the starting point for a first step towards identifying populations affected by natural hazards.

2.2. Status and Trend Indices

The main basis of the proposed methodology is based on determining a series of common explanatory characteristics in each local administrative unit (hereinafter referred to as Urban Administrative Divisions, UADs) that make it possible to determine the situation and evolution over time of the population affected by landslide risk.

The most appropriate way to parameterise risk in a given area is to define a set of dimensionless indices that are easily compared to each other and are not dependent on the size of the area under study so that such comparisons can be made easily. The main data used is population, preferably referring to the area where it exists. Generally, the population density of a UAD is determined by the population over its entire area, which is not wise. The territorial extension of a UAD does not provide a direct relationship with the population. A more specific size value must be used, such as Gross Floor Area (GFA).

Two types of population density will be defined in terms of built-up area. For a given UAD, the following are calculated:

$$TPD={\displaystyle \frac{popT\left(UAD\right)}{GFA\left(UAD\right)}}$$

(1)

$$RPD={\displaystyle \frac{popR\left(UAD\right)}{GF{A}_R\left(UAD\right)}}$$

(2)

TPD is “Total Population Density” and RPD “Risk Population Density”. The popT is the total population of the UAD and popR is the population located in the area at risk of landslides. GFA_R is the built-up area in the area at risk.

The RPD value gives an idea of the absolute concentration of risk in a given area. Assuming the population is evenly distributed over the risk area, a direct relationship could be established with the spatial probability of the extent of the landslide (Fell et al., 2008; Günther et al., 2014). In other words, the higher the RPD, the greater the likelihood of harm to the affected population. See examples in Figure 1.

Three main types of status indices or ratios are defined below. A first dimensionless indicator is the amount of risk in an area as the ratio between the built-up area at risk (GFA_R) and the total built-up area (GFA_T) for each UAD, without taking population into account. This is the “Risk Surface Ratio” or RSR.

$$RSR={\displaystyle \frac{{GFA}_R}{{GFA}_T}}\le 1$$

(3)

This index reflects the amount of risk over the total built-up area in each administrative unit. In some cases this value can be very high, especially in small inland towns in mountainous areas with medium–high susceptibility. In these cases, many of the dwellings are located in historically at-risk areas. Without good options for residential expansion outside the risk zone, this value may be stable over time.

Secondly, another dimensionless index defined as RPR or “Risk Population Ratio” can be calculated for each UAD.

$$RPR={\displaystyle \frac{popR}{popT}}\le 1$$

(4)

This is the ratio of the population affected or at risk to the total population in a given area. This index alone does not define the situation of risk of a population well, because two areas with the same RPR may have very different situations.

A third RDR or “Risk Density Ratio” has therefore been defined.

$$RDR={\displaystyle \frac{RPD}{TPD}}={\displaystyle \frac{RPR}{RSR}}$$

(5)

The RDR index is a ratio of the population density in the at-risk area to the overall population density, considering the built-up area as the distribution area (GFA). It is interpreted as an indicator of the concentration of the risk-prone population to the overall population density reference. The higher the RDR, the higher the risk concentration and the smaller the difference with the overall population density, which is undesirable. Indeed, a low RDR indicates a population not at higher risk, confirming a conservative approach to risk or that affects the population less.

Obviously, the RDR index may exceed one when RPD > TPD; in reality it is very rare that there is more population at risk. In other words, the population density in the risk zone should be as low as possible in relation to the overall density. The RDR value allows comparing localities with very different population values but similar relative risk. The RDR index is also highly correlated with the spatial probability of a damaging event occurring although by definition RPD is the variable that best estimates this probability.

To illustrate this, Figure 1 shows an example with three different UADs covering a rectangular enclosure with a built-up area of 100 area units. Each circle represents a group of inhabitants (for the RPD calculation, 1 inhabitant). The three defined indices are calculated together with RPD under different scenarios: a) change in population distribution; b) change in risk area size; c) change in number of inhabitants.

In these examples, the situation of the UAD3 is worth highlighting. With a higher population density, it has a lower RDR value than UAD1. This is explained by its better behaviour, as it offers greater containment of population occupation in the risk zone. Similarly (see UAD2 and UAD3), an increase in the population at risk can lead to a decrease in RDR (even with an increase in RPD), which can be explained by the same reasoning as above. In short, a large population may have a lower RDR than a small population with the same RPD because the population is better spread out avoiding the at-risk zone.

High RSR values are paired with high RPR values and therefore high RSR values. These cases occur in small municipalities which are at high risk and which have virtually identical population density throughout the area. In these cases, the RDR value will be close to 1. In large populations, with greater discrepancy in population densities, the value of RSR and RPR is usually low and, consequently, so is RDR. This is consistent with the explanation in the paragraph above.

Determining the evolution of the RDR index over time is an essential aspect to indicate the risk variation trend type. Indeed, determining whether there is an upward trend in risk is crucial as it is signalling an increase in the population at risk, a situation that should be avoided and controlled. Therefore, a new time-type index is defined as the slope of the regression line of the RDR (or mRDR) index for a given time period. This slope determines the value and direction of the trend of the RDR time series.

Based on the results obtained in this paper, a period of 30 or 40 years in ten-year periods can be considered to be enough to discover the stability, growth, or decrease of risk in a given area. These decennial series are more stable than the annual series and categorically determine the direction of the historical evolution of risk in a given UAD.

As the RDR is defined, if for a given period of time and in the same area this index increases along with the population, this can only be because there is an increase in the population at risk being affected greater than what is experienced by the overall population. Conversely, population density can also decrease, as is the case with depopulation in small municipalities located in poorly connected mountain valleys. In this case, if the decline is more pronounced in total population density than in risk concentration, the RDR value will increase. This situation is unusual, as normally in a municipality suffering from depopulation the population density and the population affected by risk is usually similar, especially if the municipality has a large area at risk (high RSR).

Figure 2 shows an example of two very different situations in which the RDR value increases. In 2.a., the UAD is affected by the phenomenon of depopulation with a loss of inhabitants, but in a more pronounced way in the whole area than in the risk area. In this case, the RDR’s trend is upward (positive mRDR) as the risk concentration decreases at a slower rate than the overall risk concentration. This is the case in Figure 1 with UAD2 and UAD3, and discussed above. The situation is improving since the RPD decreases, but the fact that the population at risk remains more stable than the overall population should be studied.

Figure 2b shows a typical situation of an UAD increasing its risk index. Here, the density of the population at risk increases more than the overall population density, resulting in an increase in the RDR value. In short, this situation is undesirable.

To conclude this section, the Table 1 is included with the indices defined that will be applied in the rest of this paper.

3. Case Study: Valencian Community

The Valencian Community is an autonomous community in Spain located in the east and south-east of the Iberian Peninsula along the coast of the Mediterranean Sea. Covering 23,255 km², it is the eighth largest region in Spain in terms of size and is 4.60% of the national territory. The interior is mountainous, with peaks over 1,800 m high. Its relief is determined by its proximity to the sea, with a fluvial system that is embedded in the mountainous headwaters and expands through alluvial plains towards the Mediterranean.

This territory, which is narrow and stretches along the Greenwich meridian, is located between the provinces of Tarragona to the north and Murcia to the south, bordered on the east by the Mediterranean. It comprises the provinces of Castellón (CST), Valencia (VLC), and Alicante (ALC), which make up a total of 542 municipalities or equivalent UAD in Spain. It is organised into three well-defined units, namely the coastal strip, the intermediate strip, and the mountainous interior. The latter is where landslide events occur most frequently and it suffers from a more constant and marked rural exodus. The coastal strip, which in Alicante and Castellón also has coastal mountain ranges, is where the population and recognised tourist destinations are mainly located.

The size and population for 2023, according to the Spanish National Statistics Institute, is as follows: Alicante covers 5,820 km2 and has 1,955,268 inhabitants, Castellón covers 6,637 km2 and has 604,086 inhabitants, and Valencia covers 10,810 km2 and has 2,656,841 inhabitants.

The maps below (Figure 3) shows the distribution of the study area on the terrain and the RPD population density in inhabitants per 100 m² of built-up residential area (GFA).

3.1. Data Used

3.1.1. Susceptibility Mapping

The calculation of the danger level was based on the landslide susceptibility map (LSM) developed in a previous study (Cantarino et al., 2019). Its characteristics are 25×25 m pixels as the surface area unit and spatial multi-criteria evaluation (SCME) to weight the factors to obtain the susceptibility values. The three significant factors used were slope, lithology, and land cover.

The susceptibility class thresholds defined in the above-mentioned previous study were used, obtained by means of an objective and thorough classification based on a Receiver Operating Characterisation (ROC) analysis, which uses the intrinsic variability of the data. From this map, the three highest susceptibility levels (Landslide Susceptibility Index, LSI) were selected to calculate the affected cadastral areas. These were: LSI3 or L3 (medium), LSI4 or L4 (high), and LSI5 or L5 (very high).

3.1.2. Population

The Spanish National Statistics Institute (INE) is the official body that provides population data but it does not have complete data for series prior to 1981. This body offers very complete data from 1981 onwards and every 10 years thereafter, since that is when population censuses for the whole of Spain are carried out. But for previous years there is no geolocalised population data by municipality.

As a result, resort was made to a series of raster distributions of population from outside of official bodies. An example of the progress being made in this field on a global scale, through processes of downscaling the municipal population to the grid, is the Global Human Settlement Layer (GHSL) grid population data version R2023A(Joint Research Centre Data Catalogue - GHS-POP R2023A - GHS Population Grid Multitemporal... - European Commission, n.d.) , with a resolution of 1km x 1km, which provides results for the years 1975, 1980, 1990, 2000, 2010, and 2020. The GHSL is a worldwide product that had to be masked by Spain’s administrative boundaries, however, its resolution and temporal distribution is inadequate for the proposed objectives.

Fortunately, some very recent work is available that does fit this study’s approach. Indeed, Goerlich (2025) offers a grid with 100x100m resolution and decennial periodicity between 1900 and 2021 (HIPGDAC-ES), which is based on previous studies by the author and on the extensive disaggregation of the Spanish cadastre by (Uhl et al., 2023). Population was allocated to each cell by applying dasymetric techniques based on the data provided in the above-mentioned publications.

This population grid avoids population dispersion and agglomerates it in the most populated areas unlike the GHSL which disperses the population more widely, and with a much larger cell size. GHSL uses raster layers of built-up areas derived from satellite imagery, categorised according to their functional classification (residential versus non-residential) and incorporating an estimate of the building as a volume obtained by building height. The Goerlich method favours aggregation over dispersion, so many possible second homes in isolated areas are not assigned population. In contrast, GHSL offers many more inhabited cells than HIPGDAC-ES.

The concentration of the population offered by HIPGDAC-ES is useful for studies on landslides, for instance, which generally occur in peri-urban areas and which would otherwise overestimate the population affected by landslides, noting that many of them can occur at second homes.

3.1.3. Cadastre

The cadastral plot data were obtained from the Cadastral Mapping Services according to Inspire offered by the General Directorate of the Cadastre. The cadastral information adapted to the European directive Inspire is offered through interoperable services (WMS and WFS) and downloaded by municipality.

There are fields within this information necessary for this paper, such as the built-up area (GFA), the year of construction, and the type of use. The functional and residential cadastral plots have been selected, removing all those whose construction date was prior to 1800 for the initial calculation.

This date filter serves as a starting point, and has been considered adequate as a reference to start processing the cadastral information and to perform the cumulative calculation of the built-up area. Logically, it should be noted that there are few buildings dating from before 1800, especially in thinly populated municipalities. The filter in place disregards areas with limited residential buildings.

It has been observed that many cadastral plots in small towns were not registered until 1910 or 1920 even though they are older. This results in very high yet unrepresentative population density values at the beginning of the century, which has made it advisable to disregard those decades and start the calculations in 1920.

3.2. Method Implementation

In accordance with the initial approach, this paper aims to analyse the situation and evolution of the population affected by the risk of landslides over time in relation to the population historically located in susceptible areas within the Valencian Community. The line of work followed is shown in the flow chart below (Figure 4).

As indicated in the methodology section, this study seeks to locate the municipalities in which part of the population is in an area at risk of landslides. The goal is to calculate the distribution of this population within the municipality, together with the rate of area affected with respect to the total and its variation over time. For this purpose, the indices described in section 2 will be used.

An extensive period of reference decades has been selected, as extensive series of cadastral and population data are available. In this period, the 1960s–70s and 1970s–80s stand out as the point of economic growth and the beginning of tourism, as well as 2011 (end of the real estate bubble) and 2021 (the year of the pandemic). As indicated above, the starting point was taken as the year 1800 (origin “0”, previously with no risk and no computable buildings) and the two initial decades were disregarded, the first period being 1800–1920 and the second 1920–1930, already with growth values. The following are calendar decades until 2021. From 1981 onwards, the beginning of each decade coincides with the population census carried out every 10 years by the National Statistics Institute.

The long period analysed allows the evolution of the built-up areas and those at risk to be studied by examining the evolution of the growth curve and finding its singular points. It also shows the affected population according to its location between the medium and very high susceptibility levels of the LSM mapping.

The key strength of this study is that the same procedure was carried out for the large area occupied by the Valencian Community, and the results obtained were perfectly comparable, although the calculation does at times present some uncertainty associated in particular with the resolution of the population grid used.

3.2.1. Calculation of Affected Population

This calculation is based on the aforementioned 100x100m HIPGDAC-ES decennial population grid for the whole of Spain and the 25x25m LSM maps for each province. The procedure applied (see Figure 4) is to determine the population at the highest susceptibility levels and then to count it for each municipality. For this purpose, as a first step, a raster with the same resolution and datum (100x100m and LAEA) was created for the three provinces as a whole with exact cartographic adjustments.

Then the cells of the LSM raster that have a value that represents the joint value of the original cells were selected. They can be grouped into various levels of susceptibility. The affected population will be considered to be populations in the medium to very high levels, or levels 3 to 5. Therefore, the mean of the resulting cell will exactly match the value 3 or higher obtained by the LSM.

All calculations of built-up area, equivalent LSM, and affected population were carried out using scripts programmed in Python with the geoprocessor provided by ESRI’s ArcGIS. As a result, all these data are offered for the 542 municipalities in the Valencian Community, using the same calculation procedure, which can be consulted as an appendix at the end of the article.

With this approach, 361 out of a total of 542 municipalities were identified as having some kind of impact on their population in 2021. It seems a conservative result, on the safe side, as it would represent almost 71,000 affected inhabitants, or 1.5% of the almost five million total inhabitants. To mitigate some inaccuracies in the method caused by the geolocation’s resolution, 78 municipalities with a population of under 10 inhabitants were removed as non-computable residual risk. In short, 283 municipalities were considered to be at significant risk in 2021, and 239 in 1960, in both cases more than half of the Valencian Community’s total. Obviously, it cannot be considered an imminent risk for the population in these municipalities, but the results conclude that there is a group in which the risk and its prognosis are concentrated, and these are the places where the necessary prevention and surveillance should be prioritised.

Special attention has been paid to the dates cadastral plots were registered at the beginning of the last century. Indeed, it has been observed that in quite a number of smaller municipalities there is a notable lack of registered residential area, resulting in excessively high population densities. As shown in the results section, the initial series had to be shortened to avoid the appearance of anomalous data caused by this circumstance.

It is, it should be noted, also difficult to differentiate between second homes and permanent dwellings for the purposes of allocating population. Although this is a problem that cannot be resolved by this study, it has been observed that the population grid does in many cases respect these holiday areas. Indeed, it has been shown that many municipalities may have built-up areas at risk and yet have no population affected. Therefore, the population is not distributed strictly dasymetrically, and the distribution routines, when rounding the population, partially disregard those second homes located on the periphery of the town centre.

3.3.2. Obtaining Indices

Section 2 details a number of indices suitable for parameterising risk in a given area. The dimensionless status indices are RSR, RPR, and RPD, and the calculation of population densities is also included as TPD and RPD. The dynamic index mRDR is shown to be highly explanatory for determining the trend of the population affected by landslide risk in the time series considered.

The first step was to differentiate between two working periods because of the long time series and thus to be able to distinguish behaviours. A first period, known as the historical period, was defined as between 1920 and 1960, just at the limit of the demographic and tourism boom that began in 1970 and that marks an important turning point. For the current, “recent” period, the series is defined as 1970 to 2021 (the latter years coinciding with the censuses), which better reflects the current trend.

The status indices (RSR, RPR, RDR, TPD, and RPD) were calculated for the end of the two periods originally considered, i.e., 1960 and 2021. The historical and recent periods were considered to calculate mRDR, and calculated in sexagesimal degrees.

4. Results

All available data were analysed according to the methodology proposed and are summarised in the graphs below. Of particular interest are the data on the evolution of population and built-up area, both restricted to the 239 (historical series) or 283 (recent) municipalities with significant risk (Figure 5). The term “net (n)” refers to the built-up area or population directly affected by the potential for landslides; the term “total (t)” is the overall area of a UAD affected in whole or in part.

Firstly, the evolution of the built-up area is clearly marked by the country’s demographic and economic growth, as well as the beginning of tourism from 1960–70 onwards. The years 1981–2021 were chosen to calculate the trend of the recent period as the previous period of intense growth should be excluded.

The increases in both the overall area and the area at risk are quite similar. However, there are clearly different behaviours by province. In the province of Alicante, the surface area under risk increases more than the overall surface area, while in the other two provinces, construction outside the risk zone is predominant. This seems to be caused by the increased urban pressure on the province as an internationally recognised tourist destination along its coastline. Figure 6 shows the major increase in population at risk in the province of Alicante compared to the stability of Castellón.

Figure 7 offers very significant data regarding the province of Alicante’s different behaviour from Castellón. The drop in population affected right up to the start of the tourist boom in the 1960s is of note. It points to a certain process of depopulation coupled with non-expansive construction in at-risk areas. Alicante shows a notable increase as it is the most benefited by tourism. Castellón, on the other hand, tends to fall as it is not as marked a tourist and holiday destination as Alicante.

The graphs of Figure 8 highlight the high average population density in the at-risk zone at the beginning of the historical period. This can only be explained if the lack of built surface area in the first decades due to the aforementioned problem with cadastral dates is taken into account, a phenomenon which was more accentuated in the provinces of Castellón and Alicante. This situation fully justifies reducing the historical period to the 1920–1960 series. The decline in average density in the recent series is mainly due to a decline in population in the areas concerned which is not accompanied by a significant increase in built-up area. This is even more evident in the graph for the province of Castellón. The curve stabilises with a slight upward trend from 2001 onwards.

The maps in Figure 9 confirm all the indications in the previous graphs. The decrease in the affected population in the province of Castellón is notable, caused by the depopulation of inland villages, while it increases in Alicante, mainly due to the effect of tourism and second homes along the coast.

The slight trend observed in the historical series corresponds to the situation expressed in Figure 2a, but conversely: TPD decreases more than RPD (see Figure 7a and Figure 8a). This signifies a slight tendency for residential construction to avoid risky areas. In the recent series, the slight increase in RDR is due, as can be seen in Figure 7b and Figure 8b, to a change in trend, with construction occupying at-risk zones (see Figure 2a). There is a high RDR value for the province of Castellón (Figure 10), indicating a higher concentration of risk, largely as a result of including municipalities with smaller populations than the other provinces. Small municipalities have higher RDR values and also keep this value stable for a longer period of time throughout the historical series due to their lower construction and population dynamics. Finally, the small variation in RDR values across the series indicates a general continuity or stability of the risk situation.

In the early period before 1900 and in the 1900–10 decade some RSR values > 1.0 (even significantly higher) appear in some municipalities, decreasing afterwards towards normal values. This confirms that in these cases the dates of the cadastral data need to be reviewed and are not entirely reliable. All this confirms the decision taken to exclude this period to exclude Cadastre inaccuracies, selecting the series between 1920–1960 to calculate mRDR60.

The Table 2 and Table 3 below shows the averages of the defined variables and the sum of populations for the two series: historical (239 municipalities) and recent (283).

The tables above highlight the decreasing value of mRDR in the historical series, due to depopulation and lack of continued construction. The fact that the province of Castellón has less population at risk recently than historically is also of note, a finding explained by the higher incidence of risk in inland towns that have not undergone recent urban expansion, whereas in Alicante the opposite is true.

5. Discussion

With the data processed, it is important to observe how the municipalities studied are organised and what type of association there may be between them, according to the status and trend indicators calculated. To do so, a cluster analysis was run to determine the types of groupings that can be found in the area being studied. This type of analysis is a tool that has amply proven its usefulness in grouping urban areas by means of indicators (Goerlich Gisbert et al., 2017; Huang et al., 2007; Stewart & Janssen, 2014).

In a cluster analysis, when distance is used as a similarity criterion, one should avoid using highly correlated variables because they include redundant information (principle of parsimony). A prior correlation analysis was carried out to select those variables that are most suitable for describing the groups. The criterion was to prioritise the degree of correlation in the recent period and with the variables that should better organise the groups. So mRDR and those relating to the population stand out due to their low correlation with the rest of the variables.

To avoid adding complexity to Table 4, p-values indicating significant correlations when greater than 0.05 have not been included. In general, high correlation values (above approx. 0.20) are 95% statistically significant.

Table 4 shows that the only noteworthy mRDR value is with the at-risk population (0.23). Afterwards, high correlation values are observed for the rest of the variables, and between the historical and recent series in general, which seems to confirm the continuity or stability observed in Figure 10. The RDR value is highly correlated to RPD (0.72), RPR (0.83), and RSR (0.54), as expected, the latter index being the least correlated and most interesting as an indicator of magnitude. RSR, RDR, mRDR, and popT were chosen as explanatory variables for the cluster analysis, excluding the rest of the variables due to their correlation. The municipal population was also included as a variable to indicate the importance of the municipality in absolute value. It should be noted that the total population is significantly correlated to the risk indicators, but with a negative sign, i.e., municipalities with larger populations have lower risk values; however, there is a positive correlation between the total population and the population at risk (0.32).

The Ward method was chosen for the clustering method, applying the Euclidean distance. This method is considered the most appropriate since the information lost from merging elements is less than with other methodologies (Kuiper & Fisher, 1975; Sarstedt & Mooi, 2014). Moreover, this method is not very sensitive to outliers or extreme individuals and tends to form more compact clusters with similar sizes, which is an advantage when studying the behaviour of the groups of municipalities addressed.

According to the correlation matrix, the variables RSR, RDR, mRDR, and popT were selected for the two time periods defined: historical and recent. However, the population value is classically distributed according to a “heavy tail” variable and greatly alters the analysis’ outcome. Indeed, the cluster analysis is very sensitive to extreme values, which is why the three cities with a population of over 200,000 inhabitants in the two series (Valencia, Alicante, Elche) were left out.

Table 5 shows the result of the centroids of the five clusters that were determined as the minimum number to adequately describe the comparison between the two series addressed. Two different analyses were carried out for each series, using the 215 municipalities that meet all the conditions in each case, where all the variables initially used appear. It was found that four of the five clusters have similar centroids in both series, which proves the line of continuity throughout the whole period studied.

Table 6 shows the important reference values to make the groupings, leaving two very representative clusters unpaired. Firstly, the time series cluster that groups together municipalities with a clear decrease in the RDR value (c15). And then, in the recent series, a cluster with the opposite characteristics: municipalities with a marked increase in the value of RDR (c20). Table 6 shows the characteristics of each cluster along with their pairing.

These results show how the clusters of the two periods addressed are distributed, encompassing a total of 215 municipalities, which are the ones that have all the data for both periods.

In Figure 11 a 3D bubble chart shows the distribution of municipalities per cluster, the number of which is proportional to the size of the sphere. The ovals show the municipalities with differential development. In general, the evolution of risk is the natural one according to the clusters obtained, from c11 to c21, c12 to c22, ..., located on the main diagonal, which indicates that there were no major changes in this set of 145 municipalities for a century (almost 70% of the municipalities with significant risk). It should be noted that three municipalities in Alicante evolve from c15 to c20, 16 municipalities from c13 to c20, the rest of the 23 municipalities in c20 had no significant history of risk in the historical period.

Of the 43 municipalities with increasing risk (c20), Millares is notable for its unique location in a mountainous inland area and its small population. The rest are either on the coastal strip or close to large cities or major roads. Millares shows a marked decrease in the RDR in the 1980s, as more than 50% of the residential building in the recent period was registered in that decade (major economic activity due to the location of a short-lived textile factory in the 80s–90s), causing this decrease which results in a slope value that cannot be extrapolated.

Figure 12 shows two maps with the geolocation by municipality of all the clusters in the two periods. The concentration of cluster 20 along the coast of Alicante should be noted. This is the La Marina region, an area that has already been studied by the same authors (Cantarino et al., 2021) and whose results do not differ from those presented in this paper. This study, please note, was carried out by considering only the evolution of the risk assessed in economic units, without including in any way the population affected.

The result of the cluster analysis was very revealing as it established groups of municipalities in a large area that have common characteristics both in terms of their situation and their evolution. It is particularly relevant that in four of the five groups established, the characteristics of the historical series are maintained in the recent series.

One of the remaining groups is notable for bringing together dozens of municipalities that show an upward trend in risk, presumably due to their population dynamics and tourism and holiday activity. Indeed, many of them are located on the coastal strip, where areas of high residential value are located. According to the above, construction expansion may lead to the occupation of risk areas, given the lack of conventional land and the need to meet demand, often at high prices.

This is the case of well-known coastal towns in the province of Alicante with a large amount of tourism, although not all of them present the aforementioned increased risk. Factors such as building pressure in certain coastal areas lead to rising demand and consequent increases in land value, making it difficult to access land that is free of risk. Or even often not avoiding risk because the desire for scenic locations or areas closer to the coastline is so great. This problem has already been noted in the study cited above (Cantarino et al., 2021).

At the other extreme is the province of Castellón, with very few towns showing high risk. This is due to a much smaller demographic dynamic, with inland locations, poor connections, and less building pressure. It is also relevant that the coastal strip, which is also mountainous and of great tourism value, has such a different value when compared with the province of Alicante. It should be noted that a large part of this coastal mountainous strip was declared a regional natural park, which means there are major limits on construction that indirectly allow the population to avoid the risk of landslides.

The province of Valencia, to conclude, has no mountains along its coastline, with the exception of a single municipality which, however, is not affected in a major way by the risk of landslides (cluster C24). The Valencian municipalities at more risk are located in mountainous interior regions and close to the main roads.

6. Conclusions

In the scope studied, the increase in risk is mainly due to inappropriate locations of risk-prone elements, i.e., those subject to more exposure. In terms of residential building, this occurs when construction takes place in areas at risk, due to housing or holiday pressure in tourist areas.

To determine risk growth characteristics and trends, in this study a series of indices were defined, which are easy to obtain and replicable in any study of this kind. The RDR index has proven to be in particular very useful as a dimensionless descriptor of risk and adequately characterises municipalities’ risk situation. Its analysis in time series has made it possible to know its evolution, a fundamental aspect to determine the effectiveness of each municipality’s territorial management.

As a result, it was possible to observe a certain stability in the evolution of risk due to the fact that the trend in construction is mainly to build in non-risk areas. This fact is clearly evident in the historical series, a result that can be considered expected, as it is an attenuated expansion that is exclusively a product of the need for housing. In the recent period many municipalities maintain this line of stability (67% of municipalities with significant risk). However, some of them radically change their tendency to accommodate holiday homes whose pace of construction and location requirements (panoramic views, etc.) force them to locate new buildings in risk areas.

All these conclusions cannot be definitive, as the resolution of the population data needs to be improved to obtain more detailed and conclusive results. The problem of the cadastral estates at the beginning of the last century is not relevant, especially when during the first half of the century the risk was stable and allocated to inland municipalities. In terms of resolution improvement, grids of 50x50m, or, even better, 25x25m or 10x10m, would be required. Advances in high detail mapping information may lead to such maps using dasymetric techniques. The major problem lies in cataloguing second homes and being able to allocate inhabitants accordingly. Therefore, bearing in mind that these second homes tend to occupy land at risk, it may be erring on the side of excess to allocate population to second homes at risk that are actually unoccupied for a large part of the year.

This study sets out a series of guidelines together with easily obtainable indicators to determine the status and development of risk. It should be noted that studies that analyse the evolution of risk in the population in long series are not very common, and this study has allowed us to understand the behaviour of risk in previous historical series. The above findings confirm the need to study the evolution of risk over extended periods of time.

This is a pioneering inquiry that needs to be improved and extended, but it is a necessary starting point in the dynamic study of landslide risk in mountainous areas, especially those with a certain degree of holiday interest. These results therefore need to be tested and validated in other settings and with more precise data, and risk data with an economic value incorporated and combined with risk to the population to gain a complete scenario that can better meet the challenge of risk prevention.

Finally, estimating the number of affected inhabitants per municipality and their prognosis is considered essential to establish objective territorial planning procedures and to increase the resilience of the social fabric, as well as to increase the efficiency of the monitoring plans and emergency actions put in place by civil protection agencies. Indeed, by taking into account the potentially affected inhabitants, prioritisation of the identified areas will improve the accuracy of rescue operations. In light of recent disasters, some of them close to us, all protocols aimed at reducing damage to buildings and infrastructures and, of course, saving lives, should be fostered.