Delineation of Potential Groundwater Zones and Assessment of Their Vulnerability to Pollution from Cemeteries Using GIS and AHP Approaches Based on the DRASTIC Index and Speci ﬁ c DRASTIC

: The risk of aquifer contamination is determined by the interaction between the pollutant load and the vulnerability of an aquifer. Owing to the decomposition of bodies and degradation of artefacts, cemeteries may have a negative impact on groundwater quality and suitability for use due to the leaching of organic compounds (e.g., biodegradable organics, pharmaceuticals, and formal-dehyde), inorganic compounds (e.g., nitrate and heavy metals), pathogenic bacteria, and viruses. Factors such as burial and soil type, rainfall amount, and groundwater depth may increase aquifer vulnerability to pollutants generated in cemeteries. The potential for groundwater contamination was investigated in two cemeteries of the Soure region in Portugal (Samuel–UC9 and Vinha da Rainha–UC10), using the classic DRASTIC model, followed by some adjustments, depending on the particularities of the locations, resulting in a Final Classi ﬁ cation considered as Speci ﬁ c DRASTIC. By combining Remote Sensing (RS), Geographic Information System (GIS), and Analytical Hierarchy Process (AHP), groundwater potential zones (GWPZs) were identi ﬁ ed, and aquifer vulnerability was assessed, which included the elaboration of thematic maps using GIS operation tools. The maps allowed for the identi ﬁ cation of areas with di ﬀ erent susceptibilities to contamination: from “Low” to “Very high” for the DRASTIC index and from “Very Low” to “Very high” for the Speci ﬁ c DRASTIC index. Although the di ﬀ erence between the UC9 and UC10 cemeteries is negligible, UC10 is more vulnerable because of its proximity to the community and critically important mineral water resources (such as Bicanho Medical Spa). The Speci ﬁ c model seems be tt er-suited for describing vulnerability to cemeteries. Although there is limited groundwater quality data for the area, the development of vulnerability maps can identify areas that can be sensitive spots for groundwater contamination and establish procedures for pollution prevention.


Introduction
Cemeteries are places where institutional funeral practices take place and have a special meaning for storing and transforming dead bodies and serving as a collective historical memory [1][2][3][4][5][6].
People and societies have long considered contaminated groundwater near established and unplanned cemeteries to be an urgent concern [7], because it is a slow, chronic, and asymptomatic process [8,9] and should be referred to as decomposition labs [10].
Human cadavers typically contain approximately 35% organic material, 15% bone, and 50% water [11].According to the authors [11], when a person weighing about 70 kg decomposes, 30-40 L of necro leachate is released into the environment between 72 h and 3 years after death [12] and it takes between 15 and 25 years for the person to completely break down into a skeleton [11].The decomposition process releases organic materials, inorganic materials, gases, and trace elements into the groundwater, which harms both the environment and humans [9,12,13].The primary sources of pollution in cemeteries, according to Guttman et al. [1], are materials used to manufacture coffins and embalming fluids.Toxic metals (e.g., Fe, Cu, Ni, Pb, and Zn) are released into the soil by varnishes, sealants, metal handles, and decorations found on wooden coffins [11,[14][15][16][17][18][19][20].
As the world population increases and because of urban land development, cemeteries that were previously located on the outskirts of communities are now located in their centre.Many regions of the world have reported graveyard soil contaminated with extreme physical, chemical, and biological elements [21][22][23].Certain dangerous substances can linger in the atmosphere for extended periods as ultrafine or nanoparticle-sized particles [24][25][26].In the soils of urban cemeteries in Passo Fundo, Brazil, the concentrations of toxic metals were higher than those naturally occurring in control samples [2,27].Dent s [28] research at the Australian Botanical Cemetery indicates that the electrical conductivity (or salinity) around recent burials has increased noticeably.High levels of Cl − , NO3 − , NO2 − , NH₄⁺, PO4 3− , Fe, Na + , K + , and Mg 2+ were found beneath the cemetery [28].Additionally, previous research has shown that a variety of contaminants, including bacteria, viruses, phosphorus, and nitrogen, can contaminate groundwater and pose a health risk to the general public [1,2,5,10,[29][30][31][32][33][34].
Groundwater contamination is significantly influenced by burial practices, including individual or collective graves, grave depth, proximity to water sources, size of cemeteries, and number of burials.Additional factors such as coffin material, soil characteristics (e.g., lithology, mineralogy, grain size distribution, structure, thickness, leaching potential, permeability, plasticity, chemical properties, and presence of porous/fissured zones), topography, land use (e.g., presence of vegetation, agricultural practices, and urbanised areas), climatic characteristics (e.g., precipitation, temperature, and actual evapotranspiration), geological and hydrogeological aspects (e.g., groundwater flow mechanisms), abstraction rates, extension of source protection zones, depth of the water table, and seasonal fluctuations also play a crucial role in groundwater quality [2,11,14,17,31,[35][36][37][38].In cases where a cemetery is situated on permeable and porous soil, such as gravel or sand, leachate from decomposing corpses and coffin seepage can rapidly move and blend with the groundwater below [2].Optimal decomposition requires homogeneous soils with balanced proportions of sand, silt, and clay (roughly 30% of each).During periods of high precipitation, intense runoff, and infiltration, when the water table is close to the soil surface, chemicals and pathogens can swiftly migrate to the groundwater [14].It is advisable to assess cemeteries for potential risks using a comprehensive framework that takes into account risk significance, consequences, magnitude, and hazard identification [35].
Groundwater vulnerability is the tendency or likelihood of contaminants to reach a specific position in the groundwater system after their introduction at a location above the uppermost aquifer.The term first appeared in the 1970s [39] and gained notoriety in the 1980s [40].It involves intrinsic vulnerability, which refers to the characteristics that affect the migration of pollutants towards groundwater [41], and specific vulnerability, which depicts the susceptibility to a specific contaminant or group of contaminants, considering aspects such as biogeochemical attenuation processes [42].Because groundwater vulnerability cannot be measured directly [43], several indicators have been proposed to assess current groundwater quality or predict future scenarios [44,45].Taghavi et al. [46] classified these evaluation methods into four categories: (i) overlay-and index-based methods [40,42,43], (ii) process-based simulation models [43,47], (iii) statistical methods (including orthodox and Bayesian methods) [42,43,48,49], and (iv) hybrid methods [50,51].Other techniques have been put forth to assess the vulnerability of groundwater resources.These include the model of intrinsic groundwater vulnerability and specific vulnerability to pesticide pollution [52,53], techniques for determining karst aquifer vulnerability [54], an approach that incorporates impact modelling, and an index-based approach to determine how vulnerable groundwater resources are to climate change [55].
The DRASTIC index proposed by Aller et al. [40] has already been used for groundwater vulnerability assessment in many studies [56][57][58], and it can be used to assess the risk of groundwater contamination associated with cemeteries.To reduce the subjectivity of the evaluation associated with the original model, modified or updated versions have been developed to identify appropriate ratings and determine weights for the DRASTIC parameters [58][59][60][61].
As cemeteries are sensitivity places with large spatial structures, they need to have a well-designed layout to allow funeral services [62][63][64][65][66]. Generating vulnerability DRASTIC maps involves handling substantial data, and GIS tools have been employed to manipulate hydrogeomorphological, hydrogeological, soil characteristics, and land-use data [67,68].Map algebra calculations facilitate mathematical operations among thematic maps to generate composite spatial maps or charts, such as vulnerability or suitability maps [69].GIS has been previously used to develop DRASTIC-based vulnerability maps, but mainly focused on contamination risks from wastewater facilities, garbage deposits, underground gas or fuel deposits, sanitary landfills, soils contaminated by industrial activities, and agricultural soils contaminated with an excess of fertilizers (specifically nitrate) or pesticides [45,51,57,[70][71][72].For example, Sinan and Razack [73] evaluated the vulnerability of Marrakech s Haouz aquifer to various pollution sources, including Marrakech s industrial park, industrial facilities, cemeteries, and waste deposits near Ourika and Tahanaout.
The main goal of this study was to develop a DRASTIC index-based vulnerability map for assessing the risk of aquifer contamination associated with two cemeteries in the Soure Region (Portugal) using GIS interpolation tools.Using the DRASTIC method [40], a geological analysis of the study area was conducted in two phases: the first phase considered the index independently of the pollutant load, and the second phase was developed based on the locations of particularities, resulting in a Final Classification that was considered the Specific DRASTIC [74].The main innovation of this study is the use of this methodology to create maps that can be useful for defining measures to avoid groundwater pollution from cemeteries, both in existing spaces and new spaces.The DRASTIC approach was selected because it is the easiest to use and fits well with the GIS framework.Moreover, the approach is a computationally efficient model as it eliminates the need for intricate numerical analysis or multi-parameter simulation processes.What is more, though, is that it produces excellent results with little application cost.Due in part to the large number of input data points used, this methodology improves evaluation performance, thereby reducing the impact of errors on the final product.

Design of Cemeteries in Portugal
The construction of new public cemeteries was mandated by a decree dated 9 August 1834 [75].As a result of the high number of deaths in the Portuguese Civil War (1832-1833) and during the cholera epidemic of 1833, the method of burying bodies in the ground was finally regulated [75].The Decree of 21 September 1835 established that municipal authorities should allocate an area of land for the construction of cemeteries in all urban areas (villages, towns, and cities), but located at a safe distance to avoid contamination and health problems.With the new law, bodies had to be buried for 5 years in a pit made of individual soil at least 1.1 m deep and at least 0.33 m apart between graves.In 1962, Decree No. 44,220 appeared, defining the type of soil for burial: siliceous limestone, clayey limestone, or siliceous clayey, and the area must be designed for 50 years.The use of metal (zinc or lead) and solid wood coffins was prohibited to allow the bodies to degrade within five years, with the bones being able to be removed or buried deeper in the ground.Six years later (in 1968), the use of 20 L and 80 L of hydrated lime in wooden and metal coffins, respectively, was authorised to accelerate the decomposition of bodies.Graves and crypts that were unused or unmaintained for 10 years were transferred to the management of local authorities.In 1982, a new law was published (Decree-Law No. 274/82) with instructions on how to bury or cremate mortal remains.
The 1998 law (Decree-Law no.411/98) authorises that only zinc coffins can be buried in a crypt and prohibits the burial of bodies in mass graves, unless the law is revoked for special cases.The ashes of incinerated bodies can be kept in a burial urn, ossuary, or crypt or kept in the care of a family member.The graves must remain closed for at least five years, and the period may be extended if the remains are not degraded.Completely decomposed human remains can be transferred to an ossuary or a family grave or even be cremated at the request of relatives.However, if the bodies are not decomposed at the time of exhumation, they should remain closed until the skeletonization process is complete [75].

Location of Cemeteries and the Study Area
This work involved the study of two cemeteries (Samuel, with identification UC9; and Vinha da Rainha, with identification UC10), both in the municipality of Soure (central region of Portugal) (Figure 1).As no other important sources of anthropogenic pollution are known in the vicinity, these units represent the most serious threat to the quality of groundwater and public health from pollution.Whereas UC9 is situated higher up and farther away from the urban agglomeration, UC10 is situated almost flatly and nearer to urbanisations, where it may initially pose a greater risk of contamination.Despite their approximate linear distances of 2.5 Km (UC10) and 3.3 Km (UC9) from the hydromineral resource, Bicanho Spa, it may still be necessary to investigate the possibility of contamination.

Assessment of Groundwater Vulnerability
An assessment procedure consisting of three steps was devised to evaluate the possible effects of runoff from cemeteries to groundwater in the Soure area.The first step involved delineating the GWPZs.This research aims to evaluate the vulnerability of groundwater pollution using the DRASTIC index model and Specific DRASTIC technique, which was performed during the second and third phases, respectively.

Mapping of GWPZs
In this section of the study, GWPZs were defined based on a variety of geological, hydrogeological, and environmental factors using RS, GIS, and multi-criteria decision analysis (MCDA) [76][77][78][79].Pairwise comparisons can be used to solve complex decisionmaking problems by applying the AHP [79]. Figure 2 shows a flowchart that creates GWPZs using GIS.Ten thematic maps were reclassified (Figure 3): Geology, Slope, Lineament density, Drainage density (Dd), Precipitation, Land-Use/Land-Cover (LULC), Topographic Wetness Index (TWI), Stream Power Index (SPI), Distance to rivers, and Normalised Difference Vegetation Index (NDVI).
The delineation of the GWPZ map is a complex process because different environmental, climatic, and topographical factors are not widely understood [80,81].The development of RS and GIS technologies has facilitated the delineation of large-scale GWPZs [82,83].The different datasets used in the study to compute the GWPZs are detailed in Table 1.Although the spatial resolution of the satellite images provided by ESA (European Space Agency)-Sentinel-2 is generally higher, this implies more clarity and detail but also more data and storage.The Landsat 8 satellite is distinguished by the presence of thermal bands as well as band-8 (panchromatic), which is useful for improving image spectral resolution, and data are distinguished by a high radiometric resolution (16 bits), allowing the measurement of subtle variations in surface conditions.Each raster was normalised using the geometric mean criteria following the evaluation of weights using the AHP method.For every feature class, a rating value between 1 and 5 (meaning "very low", "low", "medium", "high", and "very high") was assumed [84].The rating values represented the suitability of the groundwater potential [84][85][86][87][88][89].Table 2 displays each class s normalised weight and normalised rank for each variable.The geological characteristics play a crucial role in determining groundwater potential because the hydraulic properties of the rock regulate the infiltration and percolation of water [90].The geological map of the study region was converted from vector to raster format, and three categories were created once weight and rank had been assigned (Table 2, Figure 3a): (3) Taveiro sandstones and clays, Boa Viagem sandstones, and Carrascal sandstones; (4) Cabaços limestones and marls, Cabo Mondego limestones and marls and Costa de Arnes crowded limestones; (5) Alluvium and sands and clays with kaolinite.Sedimentary rocks, such as limestones, possess substantial potential for storing groundwater.
Slopes in each area directly affect the rate of infiltration and also surface runoff, which in turn affects the recharge of groundwater, which is impacted by topography and/or slope gradient [91,92].Steep slopes decrease infiltration and groundwater recharge because they allow less water to remain on the surface for longer periods due to rapid runoff.At the same time, because of their high rates of infiltration and low runoff, flat areas are better suited for recharge [93].The slope map (degrees) was produced by using the Digital Elevation Model (DEM) and ArcMap s "Slope Tool".The study area s slope led to the creation of five categories: (1) >30, (2) 15-30, (3) 8-15, (4) 2-8, and (5) 0-2 (Table 2, Figure 3b).
Lineaments, characterised by their straight or nearly straight form, are prominent land features that are accentuated by the permeability of the soil and are widespread across the Earth s surface [94,95].Intrinsic permeability and porosity can be used to broadly characterize underlying fractures, faults, or joints [96,97].The movement and storage of groundwater as well as the facilitation of water infiltration into the subsurface depend on the lineaments [98].Following extraction, lineament discontinuities were examined using Landsat images on ArcGIS, and the "Line Density Tool" was used to create a lineament density map (Km/Km 2 ).Based on natural breaks, the following five categories were established: (1) 0-0.49, (2) 0.49-1.34,(3) 1.34-2.18,(4) 2.18-3.23,and (5) >3.23 (Table 2, Figure 3c).
Drainage is a mechanism that has an important role in controlling the hydrogeological characteristics of soils [85], and drainage density is defined as the surface area of a drained basin divided by the total length of its watercourses [99].The groundwater recharge volume is correlated with the overall length of the drainage densities [84], and a zone with a high drainage density contributes significantly to surface runoff while retaining relatively little groundwater [100].However, the drainage system is affected by several variables, including topography, climate, slope gradient, rainfall, vegetation cover, subsurface features [96], and the type and structure of the bedrock [101].This variable makes it easier to understand and assess data about groundwater infiltration, permeability, runoff potential, and relief by providing a suitable numerical measurement [94].The drainage density (Km/Km 2 ) was determined by using Equation (1) in conjunction with the Stream Network and the Line Density Tool [82].
where () represents the basin area (Km 2 ) and () is the total length of all streams in stream order  (Km).The "Hydrology Tool" in ArcMap, along with the Fill DEM, Flow Direction, Flow Accumulation, Stream Order, and Stream to Feature procedures, was used to create the Stream Network.Based on natural breaks, five categories were established: (1) 0-0.25, (2) 0.25-1.02,(3) 1.02-1.79,(4) 1.79-2.56,and (5) >2.56 (Table 2, Figure 3d).Rainfall is a hydrologic process that restores aquifers, and it is a major factor in determining groundwater potential [86].Although more recent total precipitation data were calculated at the study site, data from 1931 to 1960 [102] were used in the GIS environment because they were available in polygon shapefile format and the most recent data were contained within the polygon.Based on natural breaks, the study area s mean annual rainfall intensity was split into five zones: (1) 0-298, (2) 298-740, (3) 740-1100, (4) 1100-2070, and ( 5) >2070 (Table 2, Figure 3e).
LULC significantly influences how groundwater recharge occurs [103].Plants and trees can store water in their leaves and stems and allow it to enter the earth through their roots and rhizomes, thus contributing to recharging groundwater.This circumstance leads to the demand for groundwater extraction on agricultural and plantation land.However, the increase in the use of concrete in urban areas leads to an increase in surface runoff and a decrease in recharge.The COS2018 chart [104] provided the LULC data, and five categories were created: (1) Urban Area, (2) Bare ground, (3) Water bodies, (4) Vegetation, and (5) Agricultural (Table 2, Figure 3f).
The TWI map was created by Beve and Kirkby [105] and is the most often used map in hydrological studies [102,106].The TWI s upslope area can be used to measure subsurface lateral transmissivity or as a local slope indicator [107,108].Soil moisture content is one of the hydrological parameters that is significantly impacted by TWI in each area [109].Because the zoning and extent of saturated areas affect the occurrence of springs [107], the higher the TWI, the greater the groundwater potential.TWI calculations [110] provide an overview of how foothill, hillslope, and topographic roughness affect lateral groundwater flow.Equation (2) [111] was used to calculate TWI, which measures a cell s propensity to retain water.It also makes it easier to find favourable locations with slow runoff and concentration.
Because local alluvial layers are typically found near river courses and because sites along rivers are best-suited for effective infiltration and subsequent recharge of groundwater, the distance from hydrographic networks is significant in hydrogeological research [112].Rivers contribute to groundwater potential zones within watersheds, which in turn affects them.To begin the distance categories, the Euclidean distance tool from the ArcGIS spatial analyst tools was utilised.Based on natural breaks, the following five categories were established: (1) >858.37,(2) 567.63-858.37,(3) 332.27-567.63,(4) 138.45-33.27,and (5) 0-138.45(Table 2, Figure 3i).
The NDVI layer quantifies vegetation by measuring the difference between near-infrared light, which vegetation strongly reflects, and red light, which vegetation absorbs.It was created using ArcMap and Landsat 8 images, with water being the most likely result given the negative values.Conversely, there is a strong likelihood that it has dense green leaves if the NDVI value is near +1.On the other hand, a region with an NDVI close to zero is probably urbanised and lacks vegetation.Based on natural breaks, the following five categories were established: (1) −1-(−0.02),(2) −0.02-0.09,(3) 0.09-0.22,(4) 0.22-0.31,and (5) >0.31 (Table 2, Figure 3j).
The AHP approach was used to determine the weight of various layers.The first step was to create a Pairwise Comparison Matrix (PCM) (Table 3) using Saaty s (1-9) relative importance scale (Table 4) [113].2).To minimise the related subjectivity, the normalised weights were computed in the second step of this procedure.Equation (3) [114] was also used to calculate the sum of the values in each column, which is shown in Table 5.
where ( ) represents the variable used in the analysis and ( ) represents the PCM s total column value.
To create the Normalised Pairwise Comparison Matrix (NPCM), each column value was divided by the sum of the column values [87,114].Each variable s normalised weight (NWt) was calculated by averaging all the values in the associated row of the NPCM (Table 5) [83,93].Each normalised weight multiplied by all is equal to 1.Because the AHP method is dependent on subjective or individual judgements, its application may lead to some inconsistencies [96].The Consistency Ratio (CR) was computed to assess the accuracy.First, each PCM column was multiplied by the variable weight.The weighted sum value was then obtained by adding the values of each row.A division between the variable s weight and the weighted sum value was then performed, yielding a λ value [87].Equation (4) [115] can be used to determine the maximum eigenvalue ( ): The () values are (1) through (), and the number of criteria is ().A ( max ) value of 11.283 was found in this study.
The value of the Consistency Index (CI) was then calculated using Equation ( 5) [115]: where ( ) is the judgement matrix s maximum eigenvalue and () is the number of criteria.This study yielded a CI value of 0.143.Finally, Equation ( 6) was used to calculate the CR [115]: where, according to [115], () denotes the Random Consistency Index and () stands for CI (Table 6).A consistency ratio value of 1.51 was found in this investigation.If the CR is less than 0.10, the inconsistency is acceptable; if the CR is greater than 0.10, the judgements must be updated.The value 0.094 was determined to be a valid CR value in this investigation.
The GWPZ map was created by Equation (7), which integrates all parameters in order of significance using the Groundwater Potential Index (GWPI) [87].
(Wj) is the normalised weight of the (j ) variable, and (Xi) is the normalised weight of the variable s (i ) class.The Raster Calculator Tool in ArcGIS was used to perform the corresponding integration.Table 7 summarises the assigned normalised weights and ranks of thematic layers found for each cemetery.The cartography created for the GWPZs will be used to map the areas where aquifer recharge is favoured, as well as to define the various indices in the R parameter used in the DRASTIC index.

Mapping of DRASTIC Index Vulnerability
To map the DRASTIC index vulnerability, seven thematic maps were created: depth to groundwater (D), net recharge (R), aquifer material typology (A), soil type (S), topography (T), impact of the vadose zone (I), and hydraulic conductivity (C) [40].Each parameter was further separated into representative classes, each of which was assigned an index (i), as presented in Table 7, to correlate with the local hydrogeological characteristics (Equation ( 8)).
The weight () of each DRASTIC index parameter represents its relative importance to other attributes.The vulnerability of the aquifer to pollution increases with increasing DRASTIC index.The procedures for creating the DRASTIC-based vulnerability map are presented in Figure 4.The adopted weights and indices were proposed by Aller et al. [40] and have already been successfully validated in other works [60,69,70,[116][117][118][119][120][121][122].The seven hydrogeological layers were overlayed to produce the DRASTIC vulnerability index map using the ArcGIS raster calculator.Table 9 displays the quantitative and qualitative classifications of aquifer pollution susceptibility, which are categories modified from the values presented in Hamza et al. [57] and LNEC.In Portuguese studies, this division is the most prevalent.Equation (10) presents the computational procedure used to generate the DRASTICbased vulnerability map.It was adapted from [68] and involved arithmetic operations of maps according to Equation ( 9) and the values in Table 8 to overlap the seven thematic maps.Equation ( 10) was used to generate the value of every cell in the vulnerability map by performing an arithmetic operation.Values were stored in every cell of each thematic map.To create the vulnerability map, Equation ( 8) was added to the raster calculator function.
where M is the vector of cell values from each thematic map in line () and row (), () and () are the dimensions of the thematic grid map, () is the thematic map, () is the number of thematic maps, and () is the vector of values associated with each cell.
where S is the vector of cell values for the suitability map in lines () and () and () and () are the dimensions of the suitability grid map.

Mapping of Specific DRASTIC Vulnerability
Changes were then introduced to the specific DRASTIC [124], considering the current groundwater abstractions in the territory, in two phases: (i) first, the lithological units were classified according to the classic DRASTIC index (DI), with the same values that define potential vulnerability, degree of vulnerability, and qualitative vulnerability class (Table 10); (ii) second, the various units were reclassified according to the location of water catchments and springs; factors considered: (a) presence of geological singularities (OGSs), such as lithological contacts, veins, faults, and fractures, with real or potential connection to water catchments and aquifers; (b) location of the geological unit (LGU) in relation to water catchments and springs.Depending on the circumstances listed below, detailed reclassification may result in a higher or lower degree of classification; the final classification was identified as Specific DRASTIC to differentiate it from the typical situation.Several scenarios were considered for the OGS Factor: (i) the unit maintains the vulnerability class under the general DRASTIC index (DI) if there were no discontinuities with an actual or potential connection to the water abstractions and aquifers; (ii) if there were discontinuities or springs connected to the aquifer, these locations were classified as "very high" to "extremely high" vulnerability (G = 7 to 8); (iii) if there were discontinuities or locations with the potential for springs, these areas were classified as "medium" to "high" vulnerability (G = 5 to 6, Table 10).The following scenarios were considered when it came to the LGU Factor: (i) if the unit was located upstream of areas that either currently or potentially discharge water (groundwater abstractions and springs), it must maintain its vulnerability class under the general DRASTIC index (DI); (ii) if the unit was located downstream of established or prospective areas of natural groundwater discharge, the unit s class must be lower than the overall DRASTIC index (DI).In certain cases, the unit class may even go from high vulnerability to low vulnerability, depending on the details of each case.Units with a very high vulnerability index may not be accessible to other types of occupation and, as long as the quality of underground resources is maintained, it is often necessary to use the entire area.
When applying the methodology to the case study, the main objective was to maintain good water quality at the different extraction points around the cemeteries.

Development of Maps Depicting Site Characteristics
Parameter D affects the extent and degree of physical and chemical attenuation and degradation, as well as the degree of interaction between subsurface constituents and percolating pollutants.Parameter D was estimated using lithology and correlation with studies carried out in that area [124,125].
Parameter R signifies the volume of water that infiltrates through the ground surface and reaches the water table within a specified land area The expected recharge rate was computed by the Thornthwaite method [126] and, additionally, it will be connected to the GWPZ map.
Geotechnical properties are intricately connected to parameter A, which denotes the attenuation potential based on the lithology within the saturated zone.The geological map of Portugal, sourced from LNEG at a scale of 1:500,000 [127], supplied the necessary information for computing the partial indices A, I, and C. The importance of parameter I in determining vulnerability is due to its impact on the residence time of pollutants in the unsaturated zone, and consequently, the probability of attenuation.The capacity of the aquifer to transport water, as indicated by parameter C, affects both the hydraulic gradient and the groundwater flow.High conductivity readings indicate a high risk of contamination.Singhal and Gupta s abacus [128] was used to calculate this parameter.
A geological map of the study area is shown in Figure 5.With a geological history spanning approximately 180 million years, the Figueira da Foz region is distinguished by both the more recent geodynamic setting of the Cenozoic deposits and the numerous Mesozoic evolutionary stages of the Lusitanian Basin [129].Stratigraphic units within the study region are arranged in a substantial column extending from the Mesozoic (Upper Triassic) to the present, positioned discordantly over Precambrian and Paleozoic metasediments [129].Cemetery UC9 is located within the Cabaços limestone and marl units, categorised within the middle to upper Oxfordian period.The brown or black finegrained flint nodules are connected to the micro-sparitic and clayey micritic limestone decimetre scales deposited in freshwater limnic environments [130].This unit has a thickness of about 250 m, attitudes of N20°W, 10-15°W; the limestone expresses itself more towards the base; this unit has poor aquifer suitability despite being quite fractured.Cemetery UC10 is in the Costa de Arnes crowded limestone unit, which consists of marly limestones, limestone sandstones, and marls with a lapped surface that is concreted or piled [130].This unit has a thickness of 50 to 60 m, with a semi-parallel attitude to the previous unit; the base is composed of marls with detrital components, acquiring aquifer-aquitard characteristics.Clayish soils with a high specific surface area and cation-exchange capacity (CEC) are the most common because they maximise the retention of fluids and metals [131].For a few reasons, limestone aquifers are especially susceptible to pollution namely due to the karst morphology that they are prone to develop.Sinkholes and sinking streams are excellent ways for pollutants to seep from the topsoil into an underlying aquifer.
Figure 5. Geological map of the study area (adapted from [132]).
The cemeteries are in the Mondego River Basin.One of the Mondego River s tributaries on the left bank, the Pranto River, forms the western boundary of the municipality of Soure.It is distinguished by a low-lying area with elevations below 50 m until it joins the Mondego River at a height of roughly 2 m, known as the Vale do Pranto [124,132].The river basin receives water from multiple streams and flows over alluvium that has been deposited on clays, marls, and limestones [124,129].The extent of the fluctuations in water levels between the wet and dry seasons, along with the significant alterations in the flow of the springs through which they discharge, suggests that the self-regulating capacity of the karst aquifer systems is limited.There could be a 50-60% infiltration rate [132].UC10 and the Bicanho Medical Spa are situated in the same lithological unit.Cemetery UC9 is in the Orla Ocidental Indiferenciada da Bacia do Mondego (OOIBM) aquifer system (Figure 6).Cemetery UC10 is in the Figueira da Foz-Gesteira aquifer system (Meso-Cenozoic) (Figure 7).UC10 and the Bicanho Medical Spa are both located in the same aquifer system, as are UC9, numerous water wells, and some springs.The conceptual flow model of the Figueira da Foz-Gesteira aquifer system (a) is essentially a geological volume composed primarily of porous detrital sediments that exhibit a diverse array of textures and lenticular structures.The system appears multi-layered because the clayey layers divide the multiple aquifer units.Owing to the wide range of granulometric compositions, hydraulic properties can vary significantly between sites.Karstification also affects the transmissive and storage capacities.These aquifers are especially susceptible to pollution owing to infiltration and rapid flow through karst structures.They also have a very low capacity for self-cleaning and the rapid spread of pathogens.A conceptual model (Figure 7) proposed by Portugal Ferreira [125] for the study area suggests that recharge occurs to the NE and at higher elevations, particularly in the Cabo Mondego limestone and marl units, and then evolves to the SW until the Pranto Fault.The cemeteries are located near the Atlantic coast and nested within the climate region of the west coast.Their Köppen Geiger classification is Csb, meaning that their climate is mesothermal (humid temperate) with a long and hot dry season in July.This climate is typical of the Mediterranean region owing to the influence of the ocean [133].The coastal climate of the Mondego Basin is classified as type C2 B 2 according to the Thornthwaite climate classification [126] and it becomes wetter as the height of the basin increases.Using information gathered from the Portuguese Climate website [133], the Thornthwaite method (Figure 8) was used to determine the region s actual annual evapotranspiration.The air temperature in the study area ranges from 14.2 °C to 29.1 °C, with an average annual precipitation of 852.4 mm and actual evapotranspiration of 587.9 mm.The hydrological balance results led to the following conclusions: there is a dry period and a wet period.The first, known as the wet period, is represented by the water deficit (DH), which lasts from May/June to September, and the second, known as the water surplus (SH), extends from November to May.The water surplus is divided into two components: surface runoff (R) and underground runoff (G), resulting in SH = R + G = 264.5 L/m 2 , which is initially very modest due to the contribution from underground recharge.Parameter S assesses soil characteristics in the upper weathered zone to minimize the risk of pollution.Despite the endogenous features of cadavers influencing decomposition (e.g., age at death, cause of death, and fat content), concerns about burial soil and its impact on human taphonomy remain significant.The soil comprises an aqueous phase with dissolved elements, a gas phase, as well as biological and solid components, encompassing both organic and inorganic materials.Ferreira Gomes [124] looked closely at the data regarding each unit s soil type.The soil of both cemeteries is clay loam.This soil has a high potential for surface runoff when fully saturated.Permeability, or the ability of water to pass through soil, is either low or very low.Usually comprising less than 50% sand and more than 40% clay, they have a clayey texture.In some areas, they might also have a high potential for contraction and expansion.All soils that are less than 50 cm deep to a restrictive layer and all soils that have a groundwater table within the first 60 cm of depth were included in this group [132].Soils with a range of intermediate characteristics, such as clayey sand and sandy clay [134][135][136], are best-suited for cemetery locations.
Changes in slope affect the T parameter as they influence drainage patterns.Flatter regions have become vulnerable to contaminant flows that can reach aquifers.Using topographic data from the USGS [137] and a previously prepared DTM, the slope map (Figure 9) was developed, which delineates zones suitable for aquifer recharge and prone to infiltration of pollutants.According to the slope map, cemetery UC9 has a higher slope (6-12%) than cemetery UC10 (<2%).

Development of the Thematic Maps and the DRASTIC-Based Vulnerability Map
Based on the data presented in Table 11, seven thematic maps (Figure 10) were created and reclassified for each DRASTIC parameter.Parameter D assumes a rating of nine for both cemeteries because the unsaturated zone depth is between 1.5 and 4.6 m.The index is relatively high because possible pollutants can enter the aquifer due to the water table being relatively close to the surface (Table 11 and Figure 10a).
The hydrological balance computation and the GWPZ chart indicate that both cemeteries have an index of eight for parameter R. The two cemeteries were evaluated with a recharge rate of 178-254 mm/year because they are situated in lithological units of karst limestones, which are favoured for aquifer recharge (Table 11 and Figure 10b).
Pollutant dilution and dispersion were significantly impacted by saturation zone water content and net recharge.Parameter A is determined by the lithological material in the saturated zone; lithological unit III-Costa de Arnes crowded limestones (Upper Cretaceous, free aquifer), and lithological unit VI-Cabaços limestones and Marls (Upper Jurassic, free to confined/semi-confined aquifer), were assigned an index of 10, the maximum because they were classified as Karstified limestones, units that are extremely permeable and allow for much faster water flow within the saturated zone (Table 11 and Figure 10c).As a pollutant reservoir and filter, the soil solid phase s capacity to retain and release hazardous chemical species and microbes is essential.Both cemeteries received an index of three due to their location on clay loam soil.Clay, which is composed of smaller particles, has a greater surface area and can retain more water (Table 10 and Figure 10d).
The UC9 cemetery is situated in a convex area with medium slopes (6-12%), earning a rating of five on the T parameter.In contrast, the UC10 cemetery is situated in a flat area with a 2% slope, earning a rating of 10.Contaminants in UC10 can linger long enough on the surface to penetrate (Table 11 and Figure 10e).
In terms of parameter I, both the UC9 and UC10 cemeteries are in lithological units with a lithological index of 10, the highest level, because they are Karstified limestones with a very short contact time with the pollutant (Table 11 and Figure 10f).
Lastly, an index of eight was given to the UC9 and UC10 cemeteries concerning parameter C. For the respective units where the cemeteries are located, a hydraulic conductivity of approximately 40.7 to 81.5 m/day was estimated (Table 11 and Figure 10g).The vulnerability map (Figure 11) was then produced using the raster calculator function and Equation ( 8), weights from Figure 4, ratings from Table 8, and GIS matrix operations through Equations ( 9) and (10).The study area s values ranged from 105 (extremely low vulnerability) to 197 (very high vulnerability).The cemeteries at UC9 and UC19 are in areas with values of 192 and 197, respectively, indicating extremely high vulnerability (Figure 11a).It is required that an environmental monitoring programme is implemented at cemetery UC10, akin to that described in Directive 1999/31/EC [138], for groundwater uses (e.g., wells, holes, springs, and hot springs).With less clay present, weaker, less purifying soils at the surface, and a location in a flat area where pollutants are more likely to seep into the aquifer, UC10 is more vulnerable than UC9.The following factors contributed to vulnerability: The occupation rate map (Figure 11b) was constructed from the cemetery surface area (TSC), grave surface area (SAB) (typically 2.6 m × 1.5 m, length × width), and occupancy rate (SAB/TSC ratio) observed for 2014 in [139].A flow direction map was also created to safeguard water quality (Figure 11c).Pedrosa et al. [139] note that UC9 (14.0%) has a marginally lower occupancy rate than UC10 (15.4%) for each cemetery (Figure 11b).
According to the flow direction tool, surface waters in UC9 flow to the north and in UC10 to the west (Figure 11c, blue arrows)Although the 500 m buffer applied to all georeferenced water points is still quite far from the cemeteries in question, it is important to remember that the shared aquifer units are quite close (Figure 11c).However, it is established that in this case, there is no reason for the cemeteries under study to be concerned about any water hole becoming contaminated.It is critical to remember that many homes have water extraction points that are not listed in any database.
Every lithological unit s unique DRASTIC analysis is displayed in degrees in Table 12, and the representation is shown in Figure 12.Because Unit I have no OGS and does not affect the locations of cemeteries or georeferenced water sources, it will move from degree 5 to 4. Unit II behaves similarly to Unit I, but because UC10 has a superficial flow to the west that is draining in that direction, it will keep the degree at 4, not drop it to 3. Considering that Unit III is a part of the same aquifer system as the mineral water resource at Bicanho Medical Spa, there are numerous particularities to consider.The OGS reports that there are no discontinuities or faults near the UC10.Meanwhile, the LGU notes that although the cemetery is situated downstream of the georeferenced water points, the grade of 7 will be preserved because of the significance of these same resources.Georeferenced water points are present in Unit IV, but they are all upstream of Unit 10 itself.Even though there are a lot of discontinuities in the unit in question, none of them seem to be dangerous for moving potential contaminants from UC10.The grade of five will be upheld for security concerns.There are large discontinuities and some faults cross-unit V, but they do not affect where the cemeteries are located.Upstream, there is only one georeferenced water point, so the grade will drop from 4 to 3. Unit VI contains UC9 as well as a water point that provides a public supply upstream.The degree will persist because UC9 is in an area with northerly surface runoff.Because of how far away UC9 is from the water point, it has not been raised to a higher level.Unit VII features a georeferenced water point, discontinuity, and a few faults; however, because it is situated upstream of the cemeteries, the grade will drop from 7 to 6. Unit VIII is finally distinguished from the other units by a density of lineament and two clearly defined faults crossing it.Nevertheless, they do not affect the cemeteries susceptibility to pollution.If the grade was dropped from 4 to 3, georeferenced water points were also protected.

Final Considerations
Despite the large number of studies that have evaluated the quality of groundwater under the influence of cemeteries [140][141][142], relatively few have looked at soil and unsaturated zone characteristics in their evaluations.The only study that examined cemeteries was conducted by Razack and Sinan [73], who discovered that the observed DRASTIC indices ranged from 71 to 204.By applying the alternative GOD method to create vulnerability maps in four cemeteries in Santa Maria (Brazil), Kemerich et al. [143] determined that these cemeteries were the cause of bacterial contamination in the groundwater.
Owing to urban sprawl, a growing population, and ongoing conflicts between different land uses, the number of deaths is currently rising, while the amount of available land is decreasing.According to several sources [35,56,141], cemeteries should be located 250-500 m away from sources of potable groundwater and 30 m away from water courses or springs to reduce the risk of groundwater contamination; sands underlain by impermeable layers, for example, are not suitable as a burial substrate due to their high permeability; It is advantageous to have a thick aeration layer and a deep underground water table; the ground between graves and tombs must be made watertight; they must not be situated in sloped terrain or areas susceptible to landslides; there are no water-filled graves.According to previous studies, cemeteries have a high potential for pollution, especially if improperly built [144].Cemeteries should have their surrounding groundwater and surface water quality investigated.In the absence of specific guidelines, monitoring should adhere to the Landfill Directive s best practices for water monitoring near landfill sites (Directive 1999/31/EC) [53].
Human decomposition can contaminate groundwater in the vicinity of cemeteries, but not because of any specific toxicity, but rather because it raises naturally occurring organic and inorganic substance levels to a point where the groundwater becomes unsuitable for any use [145,146].Cemetery and burial ground risk management has been researched [28,127,[147][148][149][150][151][152].Increased nutrient concentrations, particularly nitrate compounds [7,13], have been found, and groundwater has been identified as the primary cemetery pollutant receptor [28,144,[153][154][155].
The impact of numerous anthropogenic sources of pollution is the driving force behind most studies on groundwater vulnerability assessment; however, because of the location of the equipment and the lack of other notable nearby sources of pollution, this study particularly focused on pollution from cemeteries.People who use contaminated groundwater as their household water supply are at risk of spreading regional epidemics.As a result, because it contains important information, management organisations for cemeteries as well as entities in charge of environmental vulnerability and public health vigilance should replicate this study.A risk-based decision-making framework proposed by Pollard et al. [155] has been widely adopted in the UK and other European countries.
Future societal challenges will encourage the construction of technologically advanced cemeteries with digital systems (humidity, temperature, pH, and physical-chemical parameters sensors) that will allow the state of degradation of bodies to be accurately assessed without the need to open the graves.Currently, "green funeral" practices are increasing, where a tree is planted next to the buried body, and there are already forest cemeteries, eco-cemeteries, and natural memorial reserves.In the future, new contaminants will emerge related to the development of medical, industrial, and agricultural practices (so-called emerging contaminants), which will generate concern in the management of municipal services, such as cemeteries, as they can affect the water cycle.

Conclusions
This investigation made it possible to identify areas at risk of groundwater contamination from surface runoff from two cemeteries in the Soure region (Portugal), through the construction of a vulnerability map based on the DRASTIC and DRASTIC-specific indices and applying GIS tools and operations.Cemeteries can be a significant source of water contamination, particularly in vulnerable areas where this practice is the main source of pollution.The vulnerability map allowed for the identification of areas with different susceptibilities to contamination (ranging from "Low" to "Very high" for the DRASTIC index and from "Very Low" to "Very high" for the Specific DRASTIC).
Both cemeteries are in an area of high vulnerability to aquifer contamination, though UC10 is slightly more vulnerable in quantitative terms.Its location in an area with a lower slope (2%), which promotes infiltration, along with a higher drainage density and more favourable soil occupation and use-along with a higher TWI-all contribute to this.UC9 is in an area that has higher line density, SPI value, distance to rivers, and NDVI, but the environment is not as favourable for aquifer recharge and infiltration as UC10.The two cemeteries are situated in nearly identical lithological units in terms of hydraulics, which justifies the same vulnerability in terms of quality.Because the UC10 cemetery is close to the mineral resource of the Bicanho Medical Spa, within the same aquifer unit, and is a highly unique and sensitive resource, it must be closely monitored.
Hydrogeological cartography and groundwater vulnerability maps are excellent resources for helping the description, analysis, modelling, and communication of groundwater resource management.The production of maps from hydrogeological models like the DRASTIC index is made possible by the high potential of GIS for processing and analysing complex geo-referenced data.Particularly now that space is an issue in densely populated areas, GIS has shown to be an effective cemetery management tool.

Figure 1 .
Figure 1.Geographic locations of two cemeteries in the Soure region.

Figure 2 .
Figure 2. Flowchart that uses GIS to create a GWPZ map.

Figure 4 .
Figure 4. Flowchart that uses the GIS DRASTIC index to create a groundwater vulnerability map.

Figure 6 .
Figure 6.Hydrogeological map of the study area.

Figure 8 .
Figure 8. Monthly hydrological balance in the study area.

Figure 10 .
Figure 10.(a) Parameter D map; (b) Parameter R map; (c) Parameter A map; (d) Parameter S map; (e) Parameter T map; (f) Parameter I map; (g) Parameter C map of the studied area for the DRASTIC model.

Figure 11 .
Figure 11.(a) Vulnerability map to pollution-Drastic index; (b) Occupation rate map; (c) Flow direction and buffer zones to water abstraction points of the studied area.

Figure 12 .
Figure 12.(a) Vulnerability map to pollution-General Drastic index in degrees; (b) Vulnerability map to pollution-Specific Drastic index in degrees.

Table 1 .
Data used for creating GWPZ input data.

Table 2 .
Values taken for normalised weights and thematic layer classifications.

Table 3 .
Matrix for pairwise comparison of variables in the AHP method.

Table 4 .
Saaty s scale of relative importance.

Table 5 .
Normalised matrix for pairwise comparison of variables in the AHP method.

Table 6 .
Random Consistency Index (RI) values for n variables.

Table 7 .
Assigned normalised weights of thematic layers.

Table 8 .
Partial indices (Ip) used for calculating the DRASTIC index (DI) according to the various parameters and classes.

Table 11 .
Characteristics of lithological units for the development of thematic DRASTIC maps.

Table 12 .
Drastic index and Specific Drastic index in degrees.