A complex network analysis of groundwater wells in and around the Doñana Natural Space, Spain

Research papers A complex network analysis of groundwater wells in and around the Do˜ nana Natural Space, Spain R. Rodríguez-Alarc´ on a , S. Lozano b,* a Confederaci´ on Hidrogr´ afica del Guadalquivir, Avda. del Brillante, 57, 14012 C´ ordoba, Spain b Escuela Superior de Ingenieros, University of Seville, Isla de la Cartuja, 41092 Seville, Spain ARTICLE INFO This manuscript was handled by Nandita Basu, Editor-in-Chief, with the assistance of Jan Adamowski, Associate Editor Keywords: Groundwater abstraction wells Complex network analysis Scale free Transitivity Assortativity Abstraction-based impact centrality indexes Do˜ nana Natural Space ABSTRACT This paper proposes the use of Complex Network Analysis tools to study the spatial distribution and temporal evolution of groundwater abstraction wells within a region of interest. To that end a weighted undirected network connecting wells within a certain distance was built. The edge weights take into account the distance between the wells as well as their corresponding abstraction volumes. In addition to the conventional complex network metrics, two specific abstraction-based centrality indexes are proposed, based on the total impact of a groundwater abstraction well on its surrounding wells and on the total impact of surrounding wells on a specific location. The proposed approach has been applied to the area of the Do˜ nana Natural Space, one of the most important wetlands in Europe. The topology of the network reveals a large number of small, disconnected components and a few large components. This is reflected also in the degree distribution, which follows a Power Law distribution. The network has a large average clustering coefficient, a small average path length and a low level of centralisation. Global and average local network efficiency are also high. The network is assortative in terms of degree-degree and strength-strength correlations. Eigenvector centrality of the nodes reveals the location of influential wells. 1. Introduction The aim of this paper is to study the spatial distribution and temporal evolution of groundwater wells in a region of interest using complex network analysis (CNA) tools. To that end, a weighted undirected network is built connecting neighbouring wells (i.e., those that are less than a certain distance apart). The edge weights take into account both the separation of the wells and their respective abstraction volumes. The CNA methodology used has been extensively used in hydrology. These applications can be classified in several categories, depending on whether they involve water distribution networks (e.g., Simone et al., 2016, Giustolisi et al., 2019, Sitzenfrei et al., 2020), urban drainage networks and urban stormwater networks (Reyes-Silva et al., 2020, Simone et al., 2022), river networks and water resources at the river basin level (e.g., Porse and Lund, 2016, Rodríguez-Alarc´ on and Lozano, 2019, 2022), streamflow gaging stations (e.g., Sivakumar and Woldemeskel, 2014, Halverson and Fleming, 2015, Deepthi and Sivakumar, 2023), streamflow time series (e.g., Yasmin and Sivakumar, 2018, 2021a,b), rainfall network (e.g., Jha et al., 2015, Jha and Sivakumar, 2017) or extreme precipitation events and drought propagation networks (e.g., Agarwal et al., 2018, Muthuvel and Sivakumar, 2024). More closely related is Bharti et al. (2023), who used CNA to study the groundwater level at a number of wells observed during the period 1990–2009. The locations were clustered and the transfer entropy between the time series of pairs of wells within each cluster was used to compute the connection weights of a directed network. The average clustering coefficient and the distribution of node strengths for each cluster were computed. Another study that is somewhat relevant to this research is Agarwal and Xue (2020), who proposed a model-based clustering of networks based on exponential random graph models and applied it to study sulphate concentrations at 865 sampling sites. The edge weights of the corresponding directed network were defined as the change in sulphate concentration between the destination and source nodes. There is yet another type of CNA applications that are relevant to this research as they involve a certain radius of influence around the location of each node. Examples of these are pond networks (e.g., Fortuna et al., 2006, Thornhill et al., 2018), waterhole networks (e. g., O’Farrill et al., 2014) and wetland networks (e.g., Albanese and Haukos, 2017, Wright, 2010, Kim et al., 2022). How the connections between the different nodes are determined varies but, generally, two * Corresponding author. E-mail addresses: [email protected] (R. Rodríguez-Alarc´ on), [email protected] (S. Lozano). Contents lists available at ScienceDirect Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol https://doi.org/10.1016/j.jhydrol.2024.132079 Received 28 August 2023; Received in revised form 22 June 2024; Accepted 16 September 2024 Journal of Hydrology 644 (2024) 132079 Available online 27 September 2024 0022-1694/© 2024 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/bync-nd/4.0/ ). ponds or wetlands are considered to be connected if their distance (calculated centroid to centroid) is below a certain threshold (called dispersal distance). Typically, the degree distribution and how it changes with the intensity of drought, the sizes of connected components, the average path length, the average clustering coefficient, the network efficiency and betweenness centrality of the nodes are studied. From the literature review carried out it seems that neither CNA nor similar techniques have been applied so far to analyse the spatial distribution and temporal evolution of groundwater abstraction wells within a region. This question becomes of special importance when, as occurs in the case study used to illustrate the proposed approach, the region under study has a high environmental value and, in spite of its official protected status, its medium-term hydrological balance is threatened by this issue (Green et al., 2024). The goal of the study is to explore how CNA can be used to gain insight into the spatial distribution and interactions of the existing wells, to identify the areas where the concentration of wells is higher and the specific wells that are most critical as regards their potential impact on surrounding wells. Thus, the proposed approach uses the geographical location and authorised abstraction volume of registered groundwater wells in the Do˜ nana region to build a network of interconnected wells based on a threshold distance. These undirected connections are weighted to take into account the actual distance between the wells and their respective abstraction volumes, i.e., the closer the wells and the larger their abstraction volumes, the stronger the connection. The degree and strength distributions, the size of the components, the network efficiency and centrality of the nodes can be readily studied and used to characterise the spatial distribution and interactions between groundwater wells. Therefore, apart from the novelty of the proposed CNA-based approach, a major contribution of this paper is to provide a new perspective on the issue of groundwater well proliferation in the area of Do˜ nana. The structure of the paper is as follows. Section 2 presents the CNA-based methods used. Section 3 introduces and describes the region of interest. Section 4 presents and discusses the findings of the study, and, finally, Section 5 summarizes and provides a conclusion. 2. Materials and methods 2.1. Conventional CNA characterisation metrics Let (V,E)be an undirected network composed of a set of nodes V and a set of undirected edges E. Let |V| = n the number of nodes, |E| = m the number of edges and A=(Aij)nxn the adjacency matrix whose components are binary and indicate whether the edge corresponding to each pair of nodes (i,j)exists or not. Let us assume also that the network is weighted, i.e., each edge (i,j)has a weight wij. For an undirected network like the one considered in this paper both the adjacency matrix and the weight matrix are symmetric, i.e., Aij =Aji ,wij =wji ∀i,j∈V. Some basic topological network measures are: Density, i.e., proportion of all possible edges that are actually present ρ =m n⋅(n−1)/2. It measures the sparsity of the network. Node degree, i.e., number of edges connecting a node i: ki=∑jAij = ∑jAji. It measures the number of wells in the vicinity of a given node. Node neighbourhood, i.e., the set of nodes connected with each node i:Ni={j∈V:Aij =1} Node strength, i.e., the sum of weights of the edges connecting a node i: si=∑jwij =∑jwji. It is sometimes called weighted-degree as it measures the number of neighbouring wells but weighting them by the intensity of the corresponding connections. Geodesic distance matrix, i.e., the length of the shortest path between each pair of nodes (i,j); for an undirected network this matrix is symmetric, i.e., dij =dji ∀i,j∈V. These geodesic distances are computed on the unweighted network and represent the minimum number of edges to traverse in order to go from one node to another. Average path length, i.e., average geodesic distance: 〈d〉 = 1 n⋅(n−1)/2∑j>idij. This is the average of the entries of the geodesic distance matrix. It measures the expected separation, in average, between any given pair of nodes. Diameter, i.e., maximum geodesic distance: D=max j>idij. This is the maximum of the entries of the geodesic distance matrix, i.e., the geodesic distance between the two most distant nodes. It represents the overall “extent” of the network. Network efficiency, i.e., the inverse of the harmonic mean geodesic distance: η =1 n⋅(n−1)/2∑j>id−1 ij . The same as the average path length above, only that it uses the harmonic mean instead of the arithmetic average and is inverted, so that higher network efficiency implies shorter geodesic distances and vice versa. Degree centralization index, i.e., the extent to which the network has a centre: CI =∑i(kmax−ki) (n−1)(n−2)(kmax =maxi∈Vki is the maximum degree in the network). It adds the difference between the node degrees and that of the most connected node, using for normalization the corresponding value for a star. Clustering coefficient, i.e., the fraction of neighbours of a node i that are themselves connected: CCi=2 ki⋅(ki−1)∑j>hAijAihAjh =(A3)ii ki⋅(ki−1). It measures how close the neighbourhood of a node is to being a clique, i. e., a complete subgraph with connections between all the nodes. Average clustering coefficient: CCi=1 n∑i∈VCCi. This measures the overall cliquishness of the network, i.e., how meshed it is, on average, at the local level. Transitivity, i.e., the relative frequency of triangles (i.e., triadic closure) in the network: T=3⋅∑h>j>iAij Aih Ajh ∑h>j>i(Aij Aih +AjiAjh+AhiAhj)= 1 2⋅Traza(A3) ∑j>i(A2)ij . Similar to the average clustering coefficient, but computed directly on the whole network instead of as an average of local values. Local efficiency, i.e., the efficiency of the subgraph formed by the neighbours of a node i: η local i=1 ki⋅(ki−1)/2∑j,h∈Ni j>h d−1 jh ( djh are the geodesic distances in the subgraph). Similar to the local clustering coefficient, it measures the density of connections in the neighbourhood of a node. Average local efficiency: η local =1 n∑i∈V η local i. This is similar to the average clustering coefficient. It is also computed as an average of local values. Eigenvector centrality: xi=1 λ∑j∈Vwij xj (λ is a normalization constant). In the case of a weighted network the eigenvector centrality of a node is defined as the weighted sum of the centrality of its neighbours. Eigenvector centrality measures how “well-connected” a node is within the network. By definition, a node has high eigenvector centrality if it has strong (i.e., large-weight) connections with nodes with high eigenvector centralities. Assortativity: This measures whether nodes with high degree are more likely to be connected to nodes with high degree or to nodes with low degree (disassortativity). This can be measured by plotting the average nearest neighbour degree of a node against its degree. The average nearest neighbour degree of a node i is the average degree of its neighbours knn,i=1 ki∑j∈VAij kj=1 ki∑j∈Nikj. If such scatterplot shows a positive correlation, then the network is assortative, and it will be disassortative if the correlation is negative. Degree-degree correlation: The Pearson correlation coefficient between the degree of the two extremes of the edges is another way of measuring assortativity (or disassortativity if the coefficient is negative). Correlation coefficients measuring assortativity with respect to other node attributes can also be computed. For example, in the groundwater wells application considered in this paper the abstraction volume and the age of a well (i.e., the number of years elapsed since the well was officially registered) are node attributes and hence the assortativity with respect to these attributes indicates to what extent any pair of connected R. Rodríguez-Alarc´ on and S. Lozano Journal of Hydrology 644 (2024) 132079 2 Fig. 1. Geographical location of Guadalquivir River Basin, Do˜ nana National Park and GWB of interest. R. Rodríguez-Alarc´ on and S. Lozano Journal of Hydrology 644 (2024) 132079 3 nodes (i.e., nodes that are within a close distance) have similar abstraction volumes or were registered the same year. Power law (PL) distribution function. Mathematically, the density function of the PL distribution function is f(x) = γ−1 xmin (x xmin)−γ x≥xmin, where γ is the PL exponent and xmin is the minimum value of the variable. Some features of this distribution are that the corresponding cumulative distribution function (CDF) or complementary cumulative distribution function (CCDF) are plotted as straight lines on a log–log scale with a negative slope that depends on the PL exponent γ. Moreover, depending on the value of γ three cases can be distinguished: if γ≤2 then neither the mean nor the variance are finite, if 2 <γ≤3 then the mean is finite but the variance is not and for γ>3 both the mean and the variance are finite. Finally, note that the PL distribution function can be observed not only for the node degree distribution but also as for other variables such as the node strength or the size of the network components (i.e., the size of the maximal connected subgraphs of the network). This distribution function occurs very often in real complex networks and several mechanisms have been shown to lead to a PL degree distribution. The main one is preferential attachment, which refers to the preference of new nodes to connect with nodes of higher degree. This produces an accumulation effect that eventually leads to a heavy-tailed distribution in which a few nodes (called hubs) have very high degrees while most nodes have very low degrees. The fact that the PL distribution is scale free, i.e., its functional form does not change with the scale (or the units of measurement) of the variable, has led researchers to label this type of network as scale-free networks. Scale-free networks have been extensively studied in the literature and have a number of interesting features (such as small world character, robustness to random failures, vulnerability to target attacks, etc). 2.2. Specific features of the proposed CNA approach A weight wij is assigned to each edge so that this weight increases with the abstraction volumes of the nodes (i.e., ai and aj) and decreases with their separation dij. The overall average abstraction volume a and the separation threshold R are used to make the weight dimensionless. Specifically, wij = ai+aj 2a dij R Although alternative definitions can be conceived (e.g., using a geometric mean of the abstraction volumes instead of their arithmetic average), the justification of the use of this particular expression is that it satisfies a number of requirements, such as being as simple and parsimonious as possible, being dimensionless, generating symmetric weights and being strictly monotonic in abstraction volumes and in well separation distance. In addition, two specific abstraction-based centrality indexes are proposed for this application. The first one, labelled Ego impact (EIi), aims at measuring the aggregate effect of a well i on its neighbours Ni= {j:dij ≤R}(i.e., on those wells that are located at a distance not larger than R). This index is dimensionless and increases with the abstraction volume of the ego node and with the number and proximity of the alters EIi=ai a ⋅∑ j∈Ni R dij The second index, labelled Alters impact (AIi), aims at measuring the aggregate effect on a well i due to its neighbours. This index is also dimensionless and increases with the number, the abstraction volume and proximity of the alters AIi=∑ j∈Ni aj a⋅R dij The EIi index is related to the overall benefits of closing well i on the cones of depression of its neighbours while the AIi index measures the opposite, i.e., the aggregated effect of the neighbouring wells on the cone of depression of a well i. Moreover, the following relationship between these two indexes and the node strength holds si=EIi+AIi 2 In the case study presented in the next section, the edge length and weight distribution function and their geographical distribution are computed. This will provide information such as the mean distance between neighbouring wells or whether the weight distribution is homogeneous or has a heavy tail. Also, the macroscopic structure of the network given by the number and size of the components will be examined. It is of particular interest to know where the largest components can be found, as these represent large clusters of interconnected nodes, i.e., areas where there is a large concentration of wells. We are also interested in studying the relationship between the size of the components and the abstraction volume of the wells that belong to them and, in particular whether there is an inverse relationship between the size of the components and their corresponding average abstraction volume. This would indicate that those large concentrations of wells involve small wells, not large wells, which would have a larger impact on the corresponding cones of depression. And we want to study the degree and strength distributions, in particular whether they follow a PL distribution. The latter would occur if, for example, new wells tend to be drilled preferably close to wells with many connections, i.e., in areas with many other wells, instead of on areas with fewer wells. This might be due to the perception that it would be safer for a new well to be located in an area where there are already other wells than in a more isolated area. A PL degree distribution would also imply that only a few wells have high degree while most of the wells have very low degree. In addition to their distribution function, the geographical distribution of the node degree and strength are of interest. This will allow the identification of the areas where the nodes with large degrees (hubs) and nodes with large strength are located. The eigenvector centrality of the nodes will also be examined in order to identify where the nodes with the largest centrality indexes lie. However, although the eigenvector centrality is a standard indicator of centrality, we will also compute the two specific centrality indexes presented in section 2.2, namely Ego impact (EIi) and Alters impact (AIi), which are based on the impact that a well produces on its neighbours and the impact a well feels due to its neighbours, respectively. The geographical distribution and the correlation between these two specific centrality indexes will complement the information provided by the conventional eigenvector centrality. Finally, the assortativity of the network as regards the node degree and other node attributes will be studied, as well as the correlation between the degree and the abstraction volume and between the age of the wells and their abstraction volumes. 3. Case study and data The Do˜ nana National Park, in the South-West of Spain (see Fig. 1), was established in 1969 and covers 50,720 ha. In 1994 it was listed as a World Heritage Site and, together with the surrounding Do˜ nana Nature Park (111,645 ha), has been designated to be included in the Montreux Record of Wetlands of International Importance, i.e., Ramsar Sites. Do˜ nana is also a Natura 2000 site under the European Habitats Directive. Do˜ nana is a complex wetlands ecosystem of marshlands, seasonal ponds and lagoons and fixed and mobile dunes, which sustains a high level of biodiversity and, in particular, thousands of breeding birds and millions of migratory birds. Unfortunately, as many experts have warned (e.g., Green et al., 2024, Camacho et al., 2022, Navedo et al., 2022), excessive groundwater abstraction poses an evident risk to Do˜ nana’s groundwater resources. Thus, around 90,000 people live in R. Rodríguez-Alarc´ on and S. Lozano Journal of Hydrology 644 (2024) 132079 4 the vicinity of Do˜ nana, a large part of them engaged in intensive agriculture and tourism activities, both of which require large amounts of groundwater as well as generating pollution. The fact that current abstraction volumes (totalling around 95 hm 3 /year) are above their sustainable limits has forced the responsible government agency, i.e., Confederaci´ on Hidrogr´ afica del Guadalquivir (CHG), to declare several of the groundwater bodies (GWB) in the area as overexploited, i.e., “Bad global status” (BGS), as opposed to Good Global Status (GGS) (Confederaci´ on Hidrogr´ afica del Guadalquivir, 2022). The Almonte-Marismas aquifer arises from the superposition of an aeolian sands aquifer on top of a multilayer, partly-confined alluvial aquifer. The recharge of the aquifer occurs mainly through rainfall infiltration (estimated at 174 hm 3 /year) mainly in the aeolian sands as well as in the unconfined section of the alluvial aquifer. The overall direction of the groundwater flows within the aquifer is NW-SE. The Almonte-Marismas aquifer actually comprises six GWB in the Table 1 Groundwater bodies in the area of interest. Id. Name Area (km 2 ) Population Estimated recharge (hm 3 /y) No. of wells Aggregate abstraction volume (hm 3 ) Global status ES050MSBT000055101 Almonte 621.448 69,897 69.16 3,549 32.514 Bad ES050MSBT000055102 Marismas 515.807 5,894 33.78 142 25.537 Bad ES050MSBT000055103 Marismas de Do˜ nana 421.511 0 12.32 22 0.820 Good ES050MSBT000055104 Manto E´ olico Litoral de Do˜ nana 379.955 2,804 73.29 61 4.899 Good ES050MSBT000055105 La Rocina 493.143 9,030 80.44 1,120 32.838 Bad Total 2,431.863 87,625 268.99 4,894 96.608  Fig. 2. Growth in the number of wells and in the total abstraction volume in the area under study. R. Rodríguez-Alarc´ on and S. Lozano Journal of Hydrology 644 (2024) 132079 5 surroundings of the Do˜ nana Natural Space, of which the five that will be considered in this paper are under the responsibility of the Guadalquivir River Basin Authority (CHG) while a sixth one (namely, Condado) belongs to the Tinto, Odiel and Piedras Hydrographic District (DHTOP). These GWB have very gentle relief, with an altitude that oscillates between 150–190 m above mean sea level in the north and less than 1 m in the marshland near the coast in the south (Guardiola-Albert and Jackson, 2011, Gonz´ alez-Jim´ enez et al., 2023). The sub-humid Mediterranean climate with Atlantic influence predominates. The geological materials are of sedimentary detrital origin (gravel, gravel-sand, silt and clay) and rest on a substrate of blue marls (Guardiola-Albert et al., 2009, Guardiola-Albert and Jackson, 2011). These are detrital aquifers, permeable by primary porosity. In general, they present hydrogeological continuity with the adjacent GWB (i.e., Campo de Tejada and Gerena to the North, Alfaraje Norte, Aljarafe Sur and Sevilla-Carmona to the East). The data for the five GWB under study (see Table 1) have been obtained from the Guadalquivir River Basin Hydrological Plan (https:// www.chguadalquivir.es/segundo-ciclo-guadalquivir), approved by the Spanish Parliament (RD 1/2016) and which is in accordance with the Water Framework Directive (2000/60/CEE). The data from the registered groundwater wells that are the basis of the study come from the spatial data infrastructure provided by the Water Commission of the CHG (https://idechg.chguadalquivir.es/nodo/index.html). The data used are the official groundwater withdrawals legalized under the Spanish Water Law (RDL 1/2001, of July 20th) and included in the Water Registry created by the aforementioned Law. Therefore, the abstraction volume data is the official amount authorized and registered for each well, not the actual volume abstracted, since, at present, the agency does not have verified abstraction data per well. As regards the time span, the data include all wells officially registered by the CHG from 1986 (after the obligation to register wells imposed by the Water Law approved in 1985) up to 2020. It is important to note, however, the existence of a large number of unauthorized wells. Acreman et al. (2022) mentions that there may be more than 1,000 unauthorized wells and 1,700 suspicious irrigation ponds. As regards the extension of illegal farms, i.e., not included in the existing Special Management Plan of the Irrigation Zones, the estimations range between 1,013.6 ha (Acreman et al., 2022) and 1,232.2 ha (Gil Gil et al., 2024). Fig. 2 shows the irregular but incessant growth in the number of registered groundwater abstraction wells and in the corresponding total authorised abstraction volume during the period 1986–2020. Around 25 % of the wells and 50 % of the total abstraction volume already existed in 1990, with approximately another 25 % of the wells and of the total abstraction volume having been authorised in the last eight years of this period. Hence, although the problems of over-exploitation of the aquifer in this area started some time ago, they have been aggravated in recent years and now demand urgent action. Fig. 3 and Table 2 show the evolution of the number of groundwater wells in the last three decades in the different GWB. The thresholds used for the colour scale in Fig. 3 come from the values specified in the Spanish Water Law (Art. 54 of RDL 1/2001, July 20th) and in Ministerial Order ARM/1312/2009, May 20th (Art. 3). Not surprisingly, the large majority of the wells are located in the GWB that have been categorised as BGS. The protected status of the Do˜ nana Natural Space guarantees that practically no wells have been authorised within its limits, although there are a few wells on its border, linked to a coastal tourist resort, and which have been shown to exert a high hydrological pressure on the ponds that fall closest to this pumping area (e.g., Dimitriou et al., 2017, Fig. 3. Groundwater abstraction wells in the area under study. Table 2 Evolution of the number of groundwater wells in the last three decades. Almonte Marismas Marismas de Do˜ nana Manto E´ olico Litoral de Do˜ nana La Rocina Abstraction volume (m 3 /y) 1990 2005 2020 1990 2005 2020 1990 2005 2020 1990 2005 2020 1990 2005 2020 <7,000 596 2,023 2,920 16 25 26 1 18 18 1 33 35 35 82 597 7,000–40,0000 306 488 536 19 19 19 0 2 2 0 4 4 46 70 401 40,000–100,000 30 36 41 9 9 13 0 0 0 0 0 0 1 1 31 100,000–300,000 8 9 45 33 38 59 0 0 0 14 20 20 33 33 64 >300,000 7 7 7 22 22 25 1 2 2 2 2 2 18 19 27 TOTAL 947 2,563 3,549 99 113 142 2 22 22 17 59 61 133 205 1,120 R. Rodríguez-Alarc´ on and S. Lozano Journal of Hydrology 644 (2024) 132079 6 Fern´ andez-Ayuso et al., 2018). There are a large number of wells with small abstraction volumes, often in close proximity to one another, forming clusters that in some cases involve hundreds of wells. In spite of the high local density of nodes, the aggregate abstraction volumes of these large clusters of small wells are not significant. Much greater is the impact of some large wells (fewer in number) that, although more widely separated from each other and not clustered, represent a significant fraction of the total abstraction activity and impact. These highabstraction wells are situated on the border between Almonte and Marismas GWB and within the SE limits of La Rocina GWB. Fig. 4 shows, overlying a geological map, the geographical distribution of wells according to their depth. The thresholds used for the colour scale in Fig. 4 have been selected by the authors, based on the range of values of the observations. Most of the wells, particularly those in the populated areas, are less than 25 m deep. There are a significant number of much deeper wells (50–100 m and more) in the El Rocío and Villamanrique areas and some very deep wells (+100 m depth) in the Matalasca˜ nas area. 4. Results The Do˜ nana network of groundwater abstraction wells comprises 4,894 nodes with an average abstraction volume a=19,140.22 m3. After some initial trials, a separation threshold R=0.5 km was selected as it is a reasonable value that provides the appropriate level of granularity to carry out the study. Thus, a smaller value leads to a network with few edges and many components of very small size, and the opposite is true in the case of a larger value. For R=0.5 km the resulting network has 34,125 undirected edges. The distribution of the edge lengths is almost uniform, with an average of 0.311 km and a standard Fig. 4. Geological map and geographical distribution of well depth. R. Rodríguez-Alarc´ on and S. Lozano Journal of Hydrology 644 (2024) 132079 7 deviation of 0.125 km. Fig. 5 shows the edge weights distribution, while Fig. 6 shows its spatial distribution. The thresholds used for the colour scale in Fig. 6 were determined automatically by the QGIS software with the option “Quantile”, which balances the number of observations in each category. The edge weights seem to follow a PL distribution, with the large majority of the edges having small weights. Note that most of the edges fall within the three GWB with BGS and that the network is not connected, but formed by many scattered clusters of connected nodes. In fact, there are a total of 466 components, the largest of which involves 1,750 nodes. Fig. 7 shows the distribution of the size of these components. It can be seen that it follows a PL distribution with exponent γ=1.8117 (pvalue =0.9861975, K-S statistic =0.02101958). Moreover, as Fig. 8 shows, the largest components all fall in the Almonte GWB, with some intermediate-size components also in the La Rocina GWB. The other three GWB, which include the two with GGS, contain fewer and much smaller components. Note that in Fig. 8 each point corresponds to a well and its colour depends on the size of the component the node belongs to. Note also that the thresholds used for the colour scale in Fig. 8 were determined automatically by the QGIS software with the option “Natural Breaks (Jenks)”, which attempts to find clusters in the data by maximizing the variance between classes and minimizing the variance within each class. Fig. 9 shows the average abstraction volume and average well depth of each component with respect to its size. It can be seen that the larger components have many nodes but of a low abstraction volume. The largest wells, on the contrary, belong to very small components. This is because these large wells tend to be separated by more than 500 m due to the large cone of depression they produce. As regards the average well depth, it tends to be larger for small and for large components, with smaller values for intermediate-size components. The similarity between the two curves suggest that the abstraction volume and the depth of the wells are positively correlated. Fig. 10 shows the scatterplot of both variables. It can be seen that although the wells with the largest abstraction volumes are not very deep and the deepest wells have small abstraction volumes, there is a positive correlation between both variables (correlation coefficient 0.517). It can also be noted there are a large number of wells with low depth and small abstraction volumes. This is because, in this area, the aquifer is less productive, with less depth, less gravel and greater presence of finer material. Table 3 shows some global measures of the network. It can be noted that the network density is very low, the average path length and the Fig. 5. Distribution of edge weights. Fig. 6. Geographical distribution of the edges of the network. R. Rodríguez-Alarc´ on and S. Lozano Journal of Hydrology 644 (2024) 132079 8 Fig. 7. Distribution of component sizes. Fig. 8. Geographical distribution of the components of different sizes. R. Rodríguez-Alarc´ on and S. Lozano Journal of Hydrology 644 (2024) 132079 9 Declaration of competing interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Data availability The data that has been used is confidential. Acknowledgements This research was carried out under the Grant PID2021-123515NB- I00 funded by MCIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe”. The authors are also grateful to two anonymous reviewers for their helpful and constructive comments and suggestions, which have significantly contributed to improve this paper. References Acreman, M., Casier, R., Salathe, T., 2022. Evidence-based risk assessment of ecological damage due to groundwater abstraction; the case of Do˜ nana Natural Space, Spain. Wetlands 42, 63. Agarwal, A., Marwan, N., Maheswaran, R., Merz, B., Kurths, J., 2018. Quantifying the roles of single stations within homogeneous regions using complex network analysis. J. Hydrol. 563, 802–810. Agarwal, A., Xue, L., 2020. Model-based clustering of nonparametric weighted networks with application to water pollution analysis. Technometrics 62, 161–172. Albanese, G., Haukos, D.A., 2017. A network model framework for prioritizing wetland conservation in the Great Plains. Landsc. Ecol. 32, 115130. Bharti, V., Roshni, T., Jha, M.K., Ghorbani, M.A., Ibrahim, O.R.A., 2023. Complex network analysis of groundwater level in Sina Basin, Maharashtra, India,. Environ. Dev. Sustain. https://doi.org/10.1007/s10668-023-03375-x. Camacho, C., Negro, J.J., Elmberg, J., Fox, A.D., Nagy, S., Pain, D.J., Green, A.J., 2022. Groundwater extraction poses extreme threat to Do˜ nana World Heritage Site. Nat. Ecol. Evol. 6, 654–655. Confederaci´ on Hidrogr´ afica del Guadalquivir, Informe de estado de los acuíferos del entorno de Do˜ nana. A˜ no hidrol´ ogico 2020-2021, June 2022 (https://www. chguadalquivir.es/estudios-tecnicos last accessed January 30th 2023). Deepthi, B., Sivakumar, B., 2023. Towards assessing the importance of individual stations in hydrometric networks: application of complex networks. Stoch. Env. Res. Risk A. 37, 1333–1352. Dimitriou, E., Moussoulis, E., Díaz-Paniagua, C., Serrano, L., 2017. Hydrodynamic numerical modelling of the water level decline in four temporary ponds of the Do˜ nana National Park (SW Spain). J. Arid Environ. 147, 90–102. Fern´ andez-Ayuso, A., Rodríguez-Rodríguez, M., Benavente, J., 2018. Assessment of the hydrological status of Do˜ nana dune ponds: a natural World Heritage Site under threat. Hydrol. Sci. J. 63, 2048–2059. Fortuna, M.A., G´ omez-Rodríguez, C., Bascompte, J., 2006. Spatial network structure and amphibian persistence in stochastic environments. Proc. R. Soc. B 273, 14291434. Gil Gil, T., Schmidt, G., et al., 2024. Science to save Do˜ nana. Evidence of its ecological degradation in 2024. WWF Spain. https://wwfes.awsassets.panda.org/downloads /wwf_science-to-save-donana.pdf. Giustolisi, O., Ridolfi, L., Simone, A., 2019. “Tailoring centrality metrics for water distribution networks. Water Resour. Res. 55, 2348–2369. Gonz´ alez-Jim´ enez, M., Guardiola-Albert, C., Ezquerro, P., Aguilera, H., B´ ejar-Pizarro, M., Naranjo-Fern´ andez, N., et al., 2023. Analysis of aquifer-system deformation in the Do˜ nana Natural Space (Spain) using unsupervised cloud-computed InSAR data and wavelet analysis. Water Resour. Res. 59, e2022WR033858. Green, A.J., Guardiola-Albert, C., Bravo-Utrera, M.´ A., et al., 2024. Groundwater abstraction has caused extensive ecological damage to the Do˜ nana World Heritage Site, Spain. Wetlands 44, 20. Guardiola-Albert, C., Garcia Bravo, N., Mediavilla, C., Martin Machuca, M., 2009. Management of groundwater resources in the Donana area with the support of the Almonte-Marismas aquifer mathematical model; Gesti´ on de los recursos hídricos subterr´ aneos en el entorno de Do˜ nana con el apoyo del modelo matem´ atico del acuífero Almonte-Marismas. Boletín Geol´ ogico y Minero De Espa˜ na 120 (3), 361–376. Guardiola-Albert, C., Jackson, C.R., 2011. Potential impacts of climate change on groundwater supplies to the Do˜ nana Wetland, Spain. Wetlands 31, 907–920. Halverson, M.J., Fleming, S.W., 2015. Complex network theory, streamflow, and hydrometric monitoring system design. Hydrol. Earth Syst. Sci. 19, 3301–3318. Jha, S.K., Sivakumar, B., 2017. Complex networks for rainfall modeling: Spatial connections, temporal scale, and network size. J. Hydrol. 554, 482–489. Jha, S.K., Zhao, H., Woldemeskel, F.M., Sivakumar, B., 2015. Network theory and spatial rainfall connections: An interpretation. J. Hydrol. 527, 13–19. Kim, B., Lee, J., Park, J., 2022. Role of small wetlands on the regime shift of ecological network in a wetlandscape. Environ. Res. Commun. 4, 041006. Muthuvel, D., Sivakumar, B., 2024. Spatial propagation of different drought types and their concurrent societal risks: A complex networks-based analysis. J. Hydrol. 636, 131247. Navedo, J.G., Piersma, T., Figuerola, J., Vansteelant, W., 2022. Spain’s Do˜ nana World Heritage Site in danger. Science 376 (6589), 144. O’Farrill, G., Gauthier Schampaert, K., Rayfield, B., Bodin, O., Calme, S., Sengupta, R., Gonzalez, A., 2014. The potential connectivity of waterhole networks and the effectiveness of a protected area under various drought scenarios. PLoS One 9, e95049. Porse, E., Lund, J., 2016. Network analysis and visualizations of water resources infrastructure in California: linking connectivity and resilience. J. Water Resour. Plan. Manag. 142 (1), 04015041. Reyes-Silva, J.D., Zischg, J., Klinkhamer, C., Rao, P.S.C., Sitzenfrei, R., Krebs, P., 2020. Centrality and shortest path length measures for the functional analysis of urban drainage networks. Appl. Network Sci. 5, 1. Rodríguez-Alarc´ on, R., Lozano, S., 2019. A complex network analysis of Spanish river basins. J. Hydrol. 578, 124065. Rodríguez-Alarc´ on, R., Lozano, S., 2022. Complex network modeling of a river basin: an application to the Guadalquivir River in Southern Spain. J. Hydroinf. 24, 559–573. Simone, A., Giustolisi, O., Laucelli, D.B., 2016. “A proposal of optimal sampling design using a modularity strategy. Water Resour. Res. 52, 6171–6185. Simone, A., Cesaro, A., Del Giudice, G., Di Cristo, C., Fecarotta, O., 2022. Potentialities of complex network theory tools for urban drainage networks analysis. Water Resour. Res. 58, e2022WR032277. Sitzenfrei, R., Wang, Q., Kapelan, Z., Savi´ c, D., 2020. Using complex network analysis for optimization of water distribution networks. Water Resour. Res. 56, e2020WR027929. Sivakumar, B., Woldemeskel, F.M., 2014. Complex networks for streamflow dynamics. Hydrol. Earth Syst. Sci. 18 (11), 4565–4578. Thornhill, I., Batty, L., Hewitt, M., Friberg, N.R., Ledger, M.E., 2018. The application of graph theory and percolation analysis for assessing change in the spatial configuration of pond networks. Urban Ecosystems 21, 213–225. Wright, C.K., 2010. Spatiotemporal dynamics of prairie wetland networks: power-law scaling and implications for conservation planning. Ecology 91, 1924–1930. Yasmin, N., Sivakumar, B., 2018. Temporal streamflow analysis: Coupling nonlinear dynamics with complex networks. J. Hydrol. 564, 59–67. Yasmin, N., Sivakumar, B., 2021a. Study of temporal streamflow dynamics with complex networks: network construction and clustering. Stoch. Env. Res. Risk A 35, 579–595. Yasmin, N., Sivakumar, B., 2021b. Spatio-temporal connections in streamflow: a complex networks-based approach. Stoch. Env. Res. Risk A 35, 2375–2390. R. Rodríguez-Alarc´ on and S. Lozano Journal of Hydrology 644 (2024) 132079 16