Study Area Overview

The study area is located in the coastal plain of Rizhao City, with an average elevation below 50 m. Characterized by a temperate monsoon climate, the region experiences four distinct seasons with high humidity and frequent typhoon events. The annual mean temperature is 12.7ºC, with an average precipitation of 874 mm. The coastal plain has many rivers flowing into the sea. In the north are the Tianshui, Chaobai and Qinggang rivers. The central area has the Shadun, Gu, Futuan, and Chuanzi rivers. The south contains the Jufeng, Longwang, and Xiuzhen rivers. The mainland coastline of Rizhao City extends from the Baima River estuary in the north to the Xiuzhen River estuary in the south, with a total length of 226.35 km, including 32.45 km of natural coastline (14.33%), 167.45 km of artificial coastline (73.98%), and 26.45 km of other coastline (including estuaries and restored shoreline, accounting for 11.69%). The long and meandering coastline has complex effects on groundwater dynamics. The coastal districts of Donggang and Lanshan within the study area are densely populated coastal tourism cities and port cities. The anthropogenic influences on groundwater dynamics from nearby coastal agricultural irrigation zones and other human activities cannot be overlooked.

The groundwater primarily consists of Quaternary porous phreatic water, distributed across the mid-lower alluvial plains of river systems [21]. Additionally, bedrock fissure water is distributed in hilly and mountainous areas. The majority of the groundwater studied consists of porous phreatic water, which provides a prerequisite for selecting regional porous water monitoring wells. The dynamics of porous water are primarily influenced by atmospheric precipitation [22]. The intensive groundwater exploitation in the region has induced complex variations in both groundwater depth and aquatic environments [23]. The hydrogeological conditions of the study area are illustrated in Fig. 1.

Data Acquisition

Groundwater depth data were collected on-site from 20 monitoring wells in the coastal plain of Rizhao between 2018 and 2020. The wells were distributed across the study area to capture spatial variability. The depth was aligned with the Quaternary porous phreatic aquifer to ensure consistency. Data were recorded monthly to analyze seasonal trends.

Fig. 1.

Hydrogeological map.

Empirical Orthogonal Function (EOF) Analysis

The Empirical Orthogonal Function (EOF) method was initially proposed by statistician Pearson in the early 20^th century and was first applied to atmospheric science by Lorenz in 1950. The EOF method analyzes spatiotemporal datasets by decomposing an element field into a linear combination of spatial functions (independent of temporal variation) and temporal functions (dependent solely on time variation), thereby elucidating the spatiotemporal structure of the field. This approach has been extensively applied in geoscience research [24-26].

The EOF method constructs a spatiotemporal matrix from time-series data observed at different locations, then decomposes this matrix into the product of temporal and spatial matrices. By extracting the leading components of both spatial and temporal functions, it forms a new matrix that captures the majority of information from the original spatiotemporal field.

In this paper, the monthly groundwater depth data of 20 monitoring wells in 2018~2020 are calculated X = (x_it) (i = 1, 2, ..., m, t = 1, 2, ... n).. The m and n are the number of monitoring wells and the number of months, respectively.

Prior to decomposition, data preprocessing is essential: Temporal centering (removing station-wise means) and Optional standardization.

where x¯i represents the average water level at all moments in well i; x~it centered data (i.e., outliers that deviate from the mean); σ_i the standard deviation of the time series of water depths in well i; and x^it the relative fluctuation intensity: how many standard deviations from the mean is the current depth of burial.

The matrix X can be decomposed into a spatial function matrix A and a temporal matrix B.

a_j = (a_j1, a_j2, ..., a_jm)^T. This corresponds to the j-th empirical orthogonal mode (EOF mode), which is a spatial-dependent function exclusively determined by geographic coordinates.

Spatial Analysis

Spatial interpolation was performed using ArcGIS to map the spatial distribution of different characteristic vectors. The mean center analysis was conducted by weighting based on groundwater depth values to calculate the mean center position of the groundwater table.

Here, i = 1, 2, 3, ..., n represents the number of monitoring wells, X̅ _w is the x-coordinate of the mean center of groundwater depth, Y̅ _w is the y-coordinate of the mean center, and w_i is the weight at location i, calculated based on groundwater level depth.

Correlation Analysis

Correlation analysis between groundwater depth eigenvectors and precipitation was conducted using SPSS, with significance testing performed. The Pearson correlation coefficient (r) is calculated as follows:

Here, x is the monthly precipitation and y is the monthly groundwater depth eigenvector.

Data

The data used in this study mainly include:

(1) monthly average groundwater depth data from 20 national groundwater monitoring wells in the Rizhao coastal plain during 2018-2020; (2) monthly precipitation data calculated from daily precipitation records at the national meteorological station (Rizhao Station) from 2018 to 2020; (3) annual groundwater withdrawal data of Rizhao City for 2018-2020.

Fig. 2.

Depth map of different observation Wells.

The 20 monitoring wells are relatively uniformly distributed in the study area. Wells W10, W13, W14, and W17 are located on the right bank of the Tianshui River and the left bank of the Chaobai River in the estuary area of the northern boundary. Wells W11 and W16 are situated on the left bank of the Shadun River in the coastal region. Wells W1, W2, W12, and W20 are distributed in the estuarine area on both banks of the Futuan River in the central zone. Wells W3, W4, and W5 are located in the coastal area on both banks of the Chuanzi River. W15 is situated on the right bank of the Jufeng River. Wells W8 and W19 are positioned on both banks of the Longwang River. Wells W18 is located in the plain area between the Longwang River and the Jufeng River. Wells W6, W7, and W9 are situated in the left bank area of the Xiuzhen River at the southern boundary. The meteorological station S1 is located within the premises of the Donggang District Meteorological Bureau in Rizhao City.

Temporal Variation Characteristics of Groundwater Depth

As shown in Fig. 4, the groundwater level depth exhibits both annual and monthly fluctuations. In 2018, the maximum average groundwater depth (3.52 m) in the study area occurred in February, while the minimum (2.47 m) was observed in September. In 2019, the months of extreme values shifted, with the maximum depth ( 3.92 m) recorded in December and the minimum (3.24 m) in August. In 2020, the maximum average groundwater depth (4.46 m) in the study area occurred in April, while the minimum (2.38 m) was observed in August. The seasonal variation characteristics over the three-year period consistently showed smaller interquartile ranges of groundwater depth in summer and autumn compared to winter and spring, with groundwater depths being regionally more clustered during spring and summer seasons.

The mean center results of groundwater depth are shown in Fig. 5. Fig. 5a) displays the annual variation of average groundwater depth, revealing that in 2018, the mean center was distributed in a northeast-southwest orientation, located northwest of the study area’s center. In 2019, the mean center shifted southeastward, with an overall smaller displacement distance compared to 2018. In 2020, the mean center of groundwater depth migrated southwestward, with two distinct clusters emerging: a northeastern cluster overlapping with the 2018 and 2019 clusters, and a newly formed southwestern cluster also exhibiting a northeast-southwest orientation. The linear trend lines fitted to the mean centers of groundwater depth for different years are presented in Table 1.

The temporal coefficient curves obtained from EOF decomposition are shown in Fig. 6. Fig. 6 reveals significant variations in the temporal coefficients across different principal component modes.

The temporal coefficient of the first eigenvector (PC1) reaches its peak during March-April and its trough in August each year. The 2019 trough was nearly five times greater than those in 2018 and 2020. After reaching the minimum value in August, it rebounds sharply, while after peaking in March-April of 2018 and 2019, it declines gradually, contrasting with the abrupt decrease observed in 2020. The temporal coefficient of the second eigenvector (PC2) shows an overall downward trend over the three-year period, with slow increases after August 2018 and August 2020, but begins to fluctuate and decline monthly after August 2019. The minimum value (4.80) occurred in September 2020, while the maximum (2.90) appeared in June 2018. For the third eigenvector (PC3), the temporal coefficient showed monthly fluctuations increasing to a peak (2.80) before September 2019, followed by oscillating declines to a trough (-1.20) in March 2020. Subsequently, a peak (0.84) emerged in August 2020. The temporal coefficient of the fourth eigenvector (PC4) consistently displayed August troughs during 2018-2020 (-1.46, -2.17, and -1.83, respectively), with winter-spring peaks occurring in December 2018 (2.60) and April 2020 (0.77).

Fig. 3.

Map of monitoring sites.

Fig. 4.

Monthly statistics of groundwater depth in different years.

Fig. 5.

Annual and seasonal variations of mean central groundwater depth. a) Annual change in depth to groundwater, b) Seasonal variations in depth to groundwater.

Fig. 6.

Time coefficients of different feature vectors and monthly precipitation graphs.

Fig. 6 shows the temporal heterogeneity of precipitation in the bar chart. The annual precipitation in 2018 was 793.3 mm, while that in 2019 was 0.86 times that of 2018. In contrast, 2020 experienced more heavy rainfall events, with precipitation reaching twice the 2018 level. On average, summer rainfall accounted for 66.4% of the annual total, followed by autumn (14.5%), spring (12.7%), and winter (6.4%).

Spatial Distribution Characteristics of Groundwater Depth

Using the EOF method, the spatiotemporal field composed of groundwater depth was decomposed, and the characteristic results of the correlation matrix are shown in Table 2. The ratio represents the variance contribution corresponding to the eigenvalue, while the cumulative value indicates the accumulated variance of the preceding principal components. As shown in Table 2, the eigenvalues of the first four factors are all greater than 1, and the cumulative contribution of the first four principal components reaches 83.3%. Therefore, the first four principal components can effectively capture the majority of the key information from the original spatiotemporal field of groundwater depth. The variance contribution of the first principal component (PC1) is 36.7%, followed by the second principal component (PC2) at 24.2%, the third principal component (PC3) at 13.4%, and the fourth principal component (PC4) at 8.9%.

The first four principal eigenvectors were projected onto the corresponding monitoring well locations on the map, and isoline analysis was performed in ArcGIS to obtain spatial distribution maps of groundwater depth for different principal eigenvectors, as shown in Fig. 7.

The spatial distribution maps of different eigenvectors after EOF decomposition focus on the first four eigenvectors. All values of the first eigenvector are positive, ranging from 0.1 to 0.35. As shown in Fig. 7a), the first eigenvector exhibits three extreme value regions within the study area: a maximum value zone (PC1>0.3) located at the confluence of the Futuan River and Gu River, including monitoring wells W1 and W2; and two minimum value zones (PC1<0.1) one situated along the left bank of the lower reaches of Chuanzi River containing monitoring well W5, and another located on the right bank of Tianshui River at the northern boundary including monitoring well W14. The remaining areas in the study region show alternating distributions of high and low value zones between landward and seaward sides.

The second eigenvector exhibits three minimum value zones (PC2<-0.3): the first extreme zone appears near W14 on the right bank of the Xiuzhen River, the second minimum zone is distributed near W16 in the coastal area of the upper reaches of the Shadun River, and the third minimum zone is located near W8 in the Longwang River basin. Additionally, the second eigenvector contains one maximum value zone (PC2>0.3), which is distributed in the area near W6 on the right bank of the Xiuhuazhen River, as in Fig. 7b).

Fig. 7c) shows that the third eigenvector exhibits three maximum value zones (PC3>0.3): the first maximum zone covers a broad area near W11 in the Rizhao urban district, the second extreme zone is located near W15, and the third extreme zone is situated near W9 in the southern region. The third eigenvector also contains one minimum value zone (PC3<-0.4) located near W5 on the left bank of the Chuanzi River.

The fourth eigenvector shows three maximum value zones (PC4>0.3): the first maximum zone is situated near W14 on the right bank of the Tianshui River, the second maximum zone is located near W11 in the Rizhao urban district, and the third maximum zone is found near W7 on the left bank of the Xiuzhen River. Additionally, the fourth eigenvector contains one minimum value zone (PC4<-0.4) located near W12 on the left bank of the main channel formed by the confluence of the Futuan River and Gu River, see Fig. 7c).

The first four principal feature vectors of groundwater depth derived from EOF decomposition reveal distinct spatial patterns: PC1 (all positive values: 0.1-0.35) exhibits three extreme-value zones – a maximum zone (>0.3) at the Futuan-Gu River confluence (W1, W2), and two minimum zones (<0.1) at the lower Chuanzi River’s left bank (W5) and Tianshui River’s right bank (W14), with remaining areas showing land-sea alternation of high/low values; PC2 forms three minimum-value zones (<–0.3) along the Xiuzhen River (W14), Shadun River coast (W16), and Longwang River (W8), while displaying a maximum zone (>0.3) on the Xiuhuazhen River’s right bank (W6); PC3 contains three maximum-value zones (>0.3) in Rizhao urban area (W11), near W15, and southern region (W9), whereas the Chuanzi River’s left bank (W5) constitutes a minimum zone (<–0.4); PC4 displays three maximum-value zones (>0.3) at Tianshui River (W14), Rizhao urban area (W11), and Xiuzhen River’s left bank (W7), with a minimum zone (<–0.4) on the Futuan-Gu River mainstem’s left bank (W12). These patterns demonstrate PC1’s dominance in global spatial consistency, while higher-order vectors (PC2-PC4) reveal local heterogeneity through extreme-value distributions, their spatial differentiation being intrinsically linked to river networks and monitoring sites.

Fig. 7.

EOF decomposition of different feature vector space distribution. a) First eigenvector, b) Second eigenvector, c) Third eigenvector, d) Fourth eigenvector.

Impact of Precipitation on Principal Eigenvectors of Groundwater Depth

Fig. 6 demonstrates that summer and autumn precipitation collectively account for 80.9% of the annual total. The principal eigenvectors of groundwater depth (PC1 and PC4) exhibit synchronous fluctuation patterns with seasonal precipitation variations. Notably, when precipitation reaches its peak, PC1 and PC4 display a 1-2 months lag before attaining their trough values. This temporal delay primarily results from the time required for precipitation infiltration and runoff generation/concentration processes, which consequently leads to the hysteresis of groundwater level responses relative to precipitation variations [27-29]. PC2 and PC3 exhibit inverse trends in response to seasonal precipitation variations. Specifically, PC3 reaches its peak value with a 2-month lag following the month of maximum rainfall, while PC2 attains its trough value with a 1-2 months delay after the peak precipitation period. To quantitatively examine the relationships between precipitation and these principal eigenvectors, correlation analysis was conducted, with the results presented in Table 3.

As indicated by the correlation analysis between monthly precipitation and principal eigenvectors in Table 3, PC1 and PC4 exhibit negative correlations with monthly precipitation (r = -0.334 and -0.406, respectively), both statistically significant at the 0.05 level (two-tailed p-values: 0.046 and 0.014). Combined with the spatial analysis results from Fig. 7, this demonstrates that PC1 and PC4 represent precipitation-influenced eigenvectors that significantly contribute to groundwater depth dynamics. Consequently, the cumulative contribution of PC1 and PC4 (45.6%) quantitatively reflects precipitation’s influence on regional groundwater systems. The non-significant correlations of PC2 and PC3 (combined contribution: 54.4%) likely result from the combined effects of other natural factors and anthropogenic activities.

Impacts of Socioeconomic Development on Groundwater Resources

The Donggang and Lanshan districts in the study area primarily utilize shallow phreatic water as their water supply source. As shown in Fig. 8, the groundwater supply in Lanshan District, Rizhao City, was relatively limited - maintaining 11.25 million m³ in both 2018 and 2019, before increasing sharply to 38.45 million m³ in 2020 (a 3.42-fold increase). In contrast, Donggang District’s groundwater supply volumes were 39.65 million m³ (2018) and 38.45 million m³ (2019), with a modest increase to 41.05 million m³ in 2020.

As shown in Fig. 5, the average centroid of groundwater depth shifted southwestward in 2020, forming a new concentrated area. Combined with Fig. 8, this phenomenon can be attributed to increased groundwater extraction for water supply in the Lanshan District during 2020 compared to 2018-2019 levels, which resulted in greater groundwater depth and caused the overall southwestward migration of the study area. This demonstrates that the average centroid of groundwater depth is associated with anthropogenic groundwater extraction. To maintain the dynamic equilibrium of regional groundwater levels and protect aquifers from seawater intrusion [29- 31], it is necessary to implement both localized key controls and regional coordinated management. This involves rationally planning and dynamically adjusting the spatial distribution of socio-economic activities, as well as modifying water supply patterns and total volumes. Through policy guidance, sustainable development and utilization of groundwater resources can be achieved.

Fig. 8.

Column chart of groundwater supply in Rizhao in different years.