Land–Ocean–Atmosphere Influences on Groundwater Variability in the South Atlantic–Gulf Region

Abstract

Climate association between Groundwater Storage (GWS) and sea level changes have been missing from the Intergovernmental Panel on Climate Change, demanding a requisite study of their linkage and responses. Variability in the Hydrologic Unit Code—03 region, i.e., one of the major U.S. watersheds in the southeast caused by Sea Surface Temperature (SST) variability in the Pacific and Atlantic Ocean, was identified. Furthermore, the SST regions were identified to assess its relationship with GWS, sea level, precipitation, and terrestrial water storage. Temporal and spatial variability were obtained utilizing the singular value decomposition statistical method. A gridded GWS anomaly from the Gravity Recovery and Climate Experiment (GRACE) was used to understand the relationship with sea level and SST. The negative pockets of SST were negatively linked with GWS. The identification of teleconnections with groundwater may substantiate temporal patterns of groundwater variability. The results confirmed that the SST regions exhibited El Niño Southern Oscillation patterns, resulting in GWS changes. Moreover, a positive correlation between GWS and sea level was observed on the east coast in contrast to the southwestern United States. The findings highlight the importance of climate-driven changes in groundwater attributing changes in sea level. Therefore, SST could be a good predictor, possibly utilized for prior assessment of variabilities plus groundwater forecasting.

1. Introduction

Climate variability exhibits strong signals, mainly resulting from Sea Surface Temperature (SST) variability in hydroclimatic variables, including precipitation, temperature, groundwater, Terrestrial Water Storage (TWS), snow water equivalent, and streamflow. Hydroclimatic variables are essential elements of climate-related studies, natural disasters [1], and play a key role in the hydrological cycle [2,3]. These are documented to have a close association with low-frequency climate variability—the El Niño Southern Oscillation (ENSO) that occurs in the equatorial Pacific Ocean [4,5,6,7,8,9]. Ocean–atmosphere oscillations are naturally occurring cycles that affect SST and precipitation [10]. In addition, precipitation and sea level are the driving factors that balance terrestrial water (groundwater, soil moisture, rivers, lakes, and snow). The resulting volume of groundwater along the coast may have a close association with the Pacific and Atlantic atmospheric phenomena directly. This may be correlated with climatic response to ENSO, which drives groundwater dynamics.

Climatic variation leads to short-term fluctuations of which ENSO is the most common change in the climate. It is one of the most important climate phenomena, which can influence global temperature and precipitation by altering the global atmospheric circulation. It is also an important factor that influences groundwater dynamics. Groundwater dynamics are driven by climatic factors, including precipitation, evaporation, and temperature, as well as human factors (irrigation, extraction, and reservoir operation). Climate indices, such as Atlantic Multi-decadal Oscillation (AMO), Pacific Decadal Oscillations (PDO), ENSO, North Atlantic Oscillation (NAO), Niño Regions 3 and 3.4 SST Indices, are a good measure of the deviation from 30-year climatological normal/historical records and are predefined indices of climate variability. These deviations are aggravated by climate-induced effects of natural variability, such as ocean currents, volcanic eruptions, duration of sunshine, and earth orbital changes. For instance, adjacent lands near the ocean are influenced by the ocean currents; warm currents raise the temperatures of the coastal areas, whereas cold currents lower them. All these natural variabilities govern the global climate system and result in hydroclimatic variability. The climate phenomenon with global and regional influence is evaluated by different indices like ENSO, PDO, AMO, and NAO. Different regions of the globe have exhibited different responses in groundwater as an effect of ENSO [11] and PDO signals [12]. However, utilizing SST for the entire Pacific and Atlantic region eliminates the spatial bias impacting hydrology, further aiding in forecasting the hydroclimatic variables [8,13].

Comprehending the importance of linkage between climatic and natural variabilities that governs the fluctuations in Groundwater Storage (GWS) is a challenging task. Comparatively, studies on precipitation variability due to SST [14,15,16,17,18] and soil moisture variability [19,20] have high consideration, whereas the analysis of the relationship between climate and groundwater variability is a focused research domain. For instance, Velasco et al. [21] used groundwater level datasets to identify and quantify quasi-periodic signals in precipitation and groundwater time series that are potentially related to interannual to multidecadal modes of climate variability. Additionally, the GWS anomaly in Rajasthan, India, calculated using the Gravity Recovery and Climate Experiment (GRACE) estimates, showed an increasing trend in annual recharge [22]. GRACE is an advanced remote sensing product that provides monthly TWS measurements and has proved to be an effective tool for hydrological studies [23,24]. GRACE is being successively utilized for climate studies, quantification of hydrologic processes, and large-scale GWS changes due to limited spatial and temporal in situ groundwater data. Opie et al. [23] studied the correlation between GRACE GWS and precipitation for 37 large aquifer systems of the world. Identifying the dominant hydroclimatic variables affecting GRACE-derived groundwater variability using singular value decomposition (SVD) is a novel and the main motivation of this study. Finding prevailing factors to understand long-term groundwater variability is crucial whilst temporal fluctuation due to natural variability could not be neglected.

Estimates from the Intergovernmental Panel on Climate Change (IPCC) consider the effects of thermal expansion, ice sheets melting, and anthropogenic cause of changes in land-water storage as contributing to sea level rise. However, groundwater variability was neglected while estimating sea level budgets. According to Wada et al. [25], groundwater depletion contributes to sea level rise whereas Bjerklie et al. [26] assessed the effect of sea level rise on groundwater level under two hypothetical scenarios (0.91 m sea level rise during 100 years and 0.91 m sea level rise plus a 12-percent increase in recharge) and corroborated that water table rises with a rise in sea level and increase in recharge could exacerbate high water table issues in the New Haven coastal area. Li et al. [27] corroborated the interconnection between groundwater level and sea level changes up to the million-year scale sea level oscillations. Groundwater level change depends on the precipitation-induced infiltration and the soil-water condition, plus both are influenced by climate variability. In addition, sea level variability has not been a major area of concern regarding groundwater rise or fall. Consequently, this study aims to focus on assessing the relationship between groundwater, precipitation, and sea level over 15 years from 2002 to 2017. Additionally, the study seeks to identify the influence of climate variability on these variables to address a missing link between sea level and groundwater.

Dimension reduction is a very useful technique to manage complex data that represents the overall dataset variance of each representative cluster. Maximum information can be extracted from multiple datasets while still accomplishing a manageable analysis. For a detailed analysis, hypothesis testing can be done using statistical tools like N-way ANOVA, analysis of variance through scales, multiple regression, multiple logistic regression, Principal Component Analysis (PCA), canonical correlation analysis, and SVD. Several studies utilizing SVD to examine the relationship between hydrologic variability and SST [13,28,29], as well as precipitation and groundwater level [30], have been documented. Identifying the dominant variables affecting groundwater variability using SVD is novel and the main motivation of this study. This is a comprehensive study using SVD in GWS for the first time to assess the relationship between land–ocean and atmosphere variability.

SVD is a robust statistical tool [31] to examine the teleconnections among spatiotemporal climate variables, precipitation, and groundwater volume, to identify the linkage with oceanic–atmospheric variabilities. It is one of the common methods of PCA in studies related to hydrological variables. It uses orthogonal transformation and decomposes one or more variables into linear combinations or components into spatial and/or temporal patterns. All patterns are independent of one another. It also uses eigenvectors and evaluates the cross-covariance matrix at different spatial and equal temporal scales. SVD is a useful technique that can be used to generate alternative indices instead of relying on the existing climate indices—ENSO and PDO. Utilizing the results from an SVD analysis in this way can remove regional bias from predefined indices, as mentioned earlier [8,13]. SVD assumes linearity and it is a domain-specific method.

The current study investigated the complex issues of groundwater response to Pacific and Atlantic Ocean atmospheric circulation, like ENSO in the southeast U.S., concentrating on the HUC-03 region, which also drives variability in precipitation and sea level. This study will address the following key questions: (1) Does precipitation affect GWS variability? (2) Is variability in groundwater related to Pacific and Atlantic Ocean variability? (3) How is sea level variability related to GWS variability? These aforementioned research questions will assess the association and interactions among groundwater variability, land water storage, and the ocean. The specific objectives are (1) to evaluate the effect of precipitation on GWS in HUC-03; (2) to investigate the groundwater variability as a response of Pacific and Atlantic climate variability; and (3) to evaluate the relationship between sea level and groundwater variability across the contiguous United States.

The subsequent results from this study will help to determine interconnection amidst precipitation, sea level, and climate indices, which can be useful for short-term groundwater management decisions. Evaluating the correlation between the climate indices, precipitation, GWS, TWS, and sea level will further help in the prediction and estimation of groundwater.

2. Materials and Methods

2.1. Study Area

The study area in the South Atlantic–Gulf region, i.e., the HUC-03 watershed, was considered to evaluate the influence of precipitation and SST on groundwater response as it extends in the climatic response in the Southeast U.S., Gulf of Mexico in the west part of the HUC-03 region, and the mountain region in the northwest portion. The area covers Florida, Alabama, Georgia, North Carolina, South Carolina, a few areas of Mississippi, and Virginia. The climate of the HUC-03 watershed is subtropical; it is humid with hot summers and relatively warm but crisp winters in Virginia. The annual mean temperature ranges from 14.3 °C to 19.78 °C and the average annual precipitation ranges from 94 to 142 cm. The average annual recharge ranges from 13 to 38 cm [32].

The study region (shown in Figure 1) consists of most of the Floridian aquifer systems and is one of the world’s most productive aquifers [33]. The presence of varieties of texture, thickness, and types of carbonated rocks has the corresponding variation in regional units and zones. The Floridian aquifer system spans Southeast Georgia, Northeast and South Florida, as well as the Florida panhandle, whilst few parts of Central and Southwest Peninsular Florida is overlaid by an upper confining unit and underlain by a lower confining unit throughout the study region. The middle confining unit is a discontinuous layer that separates the Upper and Lower Floridian Aquifer. It contains confining, leaky, or semi-confining units (rocks having similar hydraulic properties within both or the adjacent aquifers) [33]. The water exchange phenomenon takes place in these units. The variability in climate, wide range of watershed sizes, and its proximity to the ocean were the basic selection criteria. These variables potentially help represent an influencing role in groundwater variability and to evaluate the hydrologic response to climatic variations and sea level variability within the study area.

2.2. Data Acquisition and Pre-Processing

The dataset used in the analysis was comprised of TWS from GRACE [34] data, Canopy Water Storage (CWS), Soil Moisture (SM), Snow-Water Equivalent (SWE) from the Global Land Data Assimilation System (GLDAS) [35], Sea Surface Temperature (SST) from the National Oceanic and Atmospheric Administration (NOAA) Physical Sciences Division [36], and monthly gridded precipitation data from the NOAA Climate Prediction Center (CPC) [37]. A continuous time series was used in this study (i.e., without missing values).

2.2.1. GRACE Dataset

GRACE, a joint U.S.–German satellite mission, measures TWS from three laboratories (i.e., Center for Space Research (CSR), Helmholtz-Zentrum Potsdam (GFZ), and NASA Jet Propulsion Laboratory (JPL)), and has a spatial resolution of 1° × 1°. The monthly GRACE–TWS data were downloaded from the GRACE FTP website [38] in NetCDF formats. Each product was multiplied by the Community Land Model version 4 (CLM4) scaling factor to obtain the average TWS. TWS constitutes a vertically integrated water storage system, including GWS, SM, SWE, and surface water.

The monthly TWS anomalies from 2002 to 2017 were extracted, which includes CWS, Ground Water Storage (GWS), SWE, and SM. A terrestrial water balance approach is shown in Equation (1) below:

$$\mathrm{TWS}=\mathrm{GWS}+\mathrm{SM}+\mathrm{SWE}+\mathrm{CWS}$$

(1)

2.2.2. GLDAS Dataset

GLDAS is a product of the Hydrological Sciences Laboratory, NASA. It provides a global dataset at temporal resolution from 1948 and a good spatial resolution of 1° × 1° and 0.25° × 0.25° and is very useful where data are limited [39,40]. It utilizes data assimilation to integrate observed and satellite data in four Land Surface Models (LSMs)—Variable Infiltration Capacity (VIC), Noah (NOAH), Mosaic (MOS), and Community Land Model (CLM). It is a widely used hydroclimatological dataset for climate and water resources research [39,40]. In addition, past studies, including Syed et al. [41], Nie et al. [42], and Moghim [43], confirmed that GLDAS TWS is consistent with GRACE TWS and is valid where it lacks a hydroclimatological dataset. In this study, the monthly SM, SWE, and CWS gridded spatial resolution of 1° × 1° grids from 2002 to 2017 were obtained from the four aforementioned LSMs determined by GLDAS.

The soil moisture obtained from NOAH and VIC was extrapolated to calculate the soil moisture depths to 3.4 m using simple linear extrapolation. The vertical layering structure of each model was model-specific, as shown in

Table 1. Vertical layering structure of the four Land Surface Models (LSMs).

Model	Layers	Depths
CLM 2.0	10	0–0.018, 0.018–0.045, 0.045–0.091, 0.091–0.166, 0.166–0.289, 0.289–0.493, 0.493–0.829, 0.829–1.383, 1.383–2.296, and 2.296–3.433 m.
MOS	3	0–0.02,0.02–1.50, and 1.5–3.50 m
NOAH	4	0–0.1, 0.1–0.4, 0.4–1.0, and 1.0–2.0 m.
VIC	3	0–0.1, 0.1–1.6, and 1.6–1.9 m.

. The average soil moisture is the depth-averaged amount of water in the soil depth of all model soil layers. For robust analysis, as suggested by Skaskevych [44], Moore and Fisher [45], and Xiao et al. [46], all four models were utilized to obtain the average of SM, SWE, and CWS.

2.2.3. Anomaly Calculation to Obtain Groundwater Storage Anomaly (GWSA)

Monthly anomalies were obtained by differencing the January 2004 to December 2009 climatological mean annual cycle, as shown in Equation (2).

$$\mathsf{\Delta }{\mathrm{P}}_{\mathrm{t}}\mathrm{A}={\mathrm{P}}_{\mathrm{t}}-\overline{{\mathrm{P}}_{2004-2009}}$$

(2)

where P_t represents GWS/SM/SWE/CWS for a particular time, t, and $\overline{{\mathrm{P}}_{2004-2009}}$ represents the corresponding average P from January 2004 to December 2009. The monthly anomalies of TWS, SM, SWE, and CWS were obtained using Equation (2). The GWSA, as estimated by Rodell and Famiglietti [47], includes the Terrestrial Water Storage Anomaly (TWSA), Snow Water Equivalent Anomaly (SWEA), Soil Moisture Anomaly (SMA), and Canopy Water Storage Anomaly (CWSA). The mass balance approach is shown in Equation (3).

$$\mathrm{GWSA}=\mathrm{TWSA}-\mathrm{SMA}-\mathrm{SWEA}-\mathrm{CWSA}$$

(3)

For the Atlantic and Pacific Analysis, the anomalies of SST, TWSA, GWSA, sea level, and precipitation were estimated as the deviation of the GRACE period (2004–2009) averages. The standardized anomalies of each dataset were used for SVD analysis by normalizing it by the standard deviation.

2.2.4. Precipitation Dataset

Precipitation data were extracted from a gridded daily (0.5° × 0.5°) precipitation database from CPC in the NetCDF format and processed into a 1° × 1° spatial resolution using a spatial averaging approach for comparison with the GRACE and GLDAS output on 1° × 1° grids. The average monthly precipitation was computed for the study period of 2002–2017. The obtained dataset was standardized before the analysis.

2.2.5. Sea Surface Temperature Dataset

The monthly SST dataset was extracted from NOAA available in 2° × 2° grid cells. The selected area in the study extended from 100° E to 80° W longitude and 30° S to 70° N latitude over the Pacific Ocean and 80° W to 20° W longitude and 30° S to 70° N latitude to represent the Atlantic Ocean. The climatology was used based on the period 2002 to 2017. The obtained dataset was standardized before the analysis.

2.2.6. Sea Level Anomaly

Monthly Sea Level Anomaly (SLA) from NOAA [48] for each month along the US coast and monthly time series was extracted. The data were extracted from 2002 to 2017. A total of 59 stations were used based on the available continuous data along the US coast. Each station has more than 30 years of sea level records. The obtained dataset was standardized before the analysis.

2.3. Methods

In this current study, each monthly precipitation (PPT), TWSA, GWSA, SLA value in the HUC-03 region, and monthly SST data of the Pacific and Atlantic Oceans were extracted to investigate the lag–lead relationship of the SST and precipitation with the other selected variables. TWSA, GWSA, and SLA will be mentioned as the HUC-03 variables hereafter. A separate analysis using SST and PPT with HUC-03 variables was carried out to investigate the most influencing variable between SST and PPT. Moreover, an SVD analysis was also utilized to study the heterogeneous relationship between the aforementioned variables at a 90% confidence level.

Climate indices based on non-atmospheric parameters, monthly SST, raw monthly mean PPT, monthly tide gauge data, monthly TWS, and monthly GWS for extended periods from 2002 to 2017 were used to investigate the association between large-scale climate effects and PPT with groundwater variability in the HUC-03 region. PPT, TWSA, GWSA, and sea-level responses to the SST anomalies were evaluated using SVD. The heterogeneous correlations were computed in hurricane months using SVD. Through these spatiotemporal analyses, the land–ocean–atmosphere interaction can be effectively studied to recognize the oscillation in climate variables. Spatial teleconnection indices were helpful to investigate the connection among groundwater with PPT, sea level, and SST.

2.3.1. Lead–Lag Relationship

The lag–lead relationship of SST with PPT, TWSA, GWSA, and SLA at each 1° × 1° grid cell was calculated over the obtained time series (2002-2017) to evaluate the monthly lag using cross-correlation analysis taking detrended SST or PPT as a reference signal and the HUC-03 variables grid cells as a 3D spatiotemporal dataset. The correlation values ranged from −1 to +1. The positive and negative correlation was extracted from a lag–lead correlation. A positive correlation is attained when the SST increases are teleconnected with increases in PPT (or TWSA, GWSA, and SLA). This can indicate that the atmosphere is forced by SST [49]. A negative correlation demonstrated that the SST is inversely teleconnected to the PPT (or TWSA, GWSA, and SLA). The magnitude of the minimum (maximum) correlation signifies the intensity of the SST response (forcing). The corresponding lag (lead) time denotes the response of SST (atmosphere) to atmospheric (SST) forcing [50]. The lag (months) and corresponding correlation to lead/lag time were evaluated in each 1° × 1° grid cell of each HUC-03 region to determine the number of lag/lead months. The lag/lead months were obtained to further compute the SVD analysis between the HUC-03 variables and SST/PPT. If a lag of 1 month was obtained, then the SST values of June were used to obtain the co-variability with the PPT/GWSA/TWSA/SLA values of July. Similarly, for the lag of two and five months, the June SST was used to obtain the co-variability with August and November, respectively. Before computing the cross-correlation analysis, the time series were detrended to remove the autocorrelation in the dataset.

2.3.2. Singular Value Decomposition

SVD is a robust statistical tool to analyze multivariate data that is commonly used to examine spatiotemporal variability between multiple variables. SVD is also used to evaluate the lagged relationship of primary variabilities at different spatial and equal temporal scales. It also uses eigenvectors similar to a PCA, and evaluates a cross-covariance matrix using Equation (4) between two fields (S and Q) and decomposes the matrix by SVD (Equation (5)) and determines the correlations at different modes.

$$\mathrm{A}=\frac{\mathrm{S}\ast {\mathrm{Q}}^{\mathrm{T}}}{\mathrm{N}}$$

(4)

Different modes of variation are obtained from the Squared Covariance Fraction (SCF), which is calculated using Equation (6).

$${\mathrm{SVD}}_{\mathrm{A}}={\mathrm{USV}}^{\mathrm{T}}$$

(5)

where U^TU and V^TV are equal to 1 and U and V are standardized matrices and S is a nonnegative singular value matrix.

$${\mathrm{SCF}}_{\mathrm{i}}=\frac{{\mathrm{C}}_{\mathrm{i}}^2}{{{\displaystyle \sum }}_{\mathrm{j}=1}^{\mathrm{r}}{\mathrm{C}}_{\mathrm{i}}^2}$$

(6)

where $\mathrm{C}$ is the singular value of the ith mode. Further, a linear correlation between two fields in all grid points is explained by the Normalized Squared Covariance (NSC) that is calculated by using Equation (7):

$$\mathrm{NSC}=\frac{C^2}{N_s\ast N_i}$$

(7)

where C² is the squared singular values sum, $N_s$ is the number of grids, and $N_i$ is the total number of independent variable stations.

2.3.3. Statistical Approach

SVD was utilized to assess the spatio-temporal relationships between the HUC-03 variables and Pacific plus Atlantic SST (or PPT). SVD was utilized to decompose the SST (or PPT) and monthly HUC-03 variables. The significant regions at a 90% confidence level were examined by utilizing the correlation values and the heterogeneous correlation maps were generated. Monthly SST and HUC-03 variables were used from 2002 to 2017 for each hurricane month. SVD assumes linearity and is a domain-specific method.

A matrix of standardized SST anomalies and standardized anomalies of each HUC-03 variable were developed. Then, a cross-covariance was generated by computing for two spatio-temporal matrices and SVD was applied to the cross-covariance matrix. Please note that an individual matrix was created for each hurricane month. The SVD analysis of the cross-covariance matrix resulted in two orthogonal matrices with singular vectors, referred to as the left and right matrix, and another one matrix of singular values. The singular value matrix was ordered such that the first mode (i.e. first singular value) was greater than the second mode, and so on. The modes corresponding to SCF greater than 10% were considered and retained. The teleconnection relationship between the two fields was obtained. SCF is a helpful measurement to compare the comparative significance of modes in the decomposition [31].

The significant spatial relationships between the left singular vector (SST) matrix and right singular vector (HUC-03 variables) matrix were evaluated by correlating the first mode or column of the left matrix (right matrix) with a standardized HUC-03 variables (SST) anomalies matrix. This resulted in the first temporal expansion series on the left and the right series. The left (right) heterogeneous correlation figures for the first mode that explained the maximum variability were determined by correlating the SST (HUC-03 variables) values of the left (right) matrix with the first temporal expansion series of right (left) field. Statistically significant correlation values (i.e. at 90% confidence level) in the heterogeneous matrices were then mapped. The temporal expansion series could represent a new index of SST variability that has a physical meaning but is excluded in existing predefined indices. This may be useful for groundwater forecasting by using high correlation values with the temporal expansion series.

The capability of SVD to correlate spatio-temporal datasets has already been confirmed and tested in surface water studies for regional assessment [13,28,29,51]. These studies were able to provide the estimated climate signals. Bretherton et al. [31] compared four different techniques (PCA, CCA, SVD, and single-field PCA) and Chitsaz et al. [52] compared SVD with PCA, and both concluded that SVD was the easiest method to implement and superior to the other methods. It is easy to implement as it untangles data into independent components that are easy to analyze. In addition, Wallace et al. [51] concluded that SVD segregates the important modes of variability. In contrast to PCA, SVD can be applied to non-square matrices, which is its main advantage. It is also an effective tool as it depends on maximizing the covariance structure of the data and not only on the correlation between two variables. It can also work well with two different spatiotemporal fields of different dimensions. SVD works well if the unit of measurement is different and differences in covariance are important.

While SVD is a powerful and well-accepted statistical approach of two spatial-temporal fields, some limitations do exist [53]. Therefore, caution must be taken while explaining the physical relationship and it should be substantiated by the literature. It may not give reasonable results for a strongly non-linear dataset, it is easy to misuse and/or misinterpret, and it has the potential to produce spurious patterns and correlations. However, according to Newman and Sardeshmukh [53], SVD can be applied if the first, second, and third mode explains the significant variance of the two fields to evaluate the strength of the coupled variability.

3. Results

The SVD analysis was utilized to investigate the association of climatic effects on the HUC-03 variables. The section is organized into two sub-sections: (3.1) Lag–Lead Correlation Analysis and (3.2) Co-Variability and Correlation Analysis using the SVD Technique.

3.1. Lag–Lead Correlation Analysis

The SST response time to PPT, TWSA, GWSA, and SLA was identified in the HUC-03 region. The SST response time evaluation using 3D spatiotemporal cross-correlation analysis showed lag/lead months in both the Pacific and the Atlantic Ocean. The lag/lead correlation analysis was fundamental for further assessment of the heterogeneous correlation using SVD analysis. The highest spatial average of the forcing time of SST was found to lead TWSA by 5 months. The SST forcing with the highest correlation of 0.5 was found to be near the southern part of Florida with the largest lead time of 5 months. The SST forcing was found to have a longer lag correlation in the western Atlantic Ocean and the northern Gulf of Mexico than in the Pacific Ocean near north of Mexico while investigating with TWSA, whereas the SST response had a shorter lag correlation in the northern Atlantic Ocean near west of South America than in the Pacific Ocean near north of Mexico. The highest spatial average was obtained where the SST led PPT by 2 months. The minimum correlation between SST and PPT was observed in the west of the Gulf of Mexico. Further, SST showed a strong response on GWSA and SLA with a lead by 1 month, as observed from the highest spatial average of the grid cells within the study area, as shown in Figure 2. The maximum correlation shifts southward in GWSA. The minimum correlation region was observed in the west of the Gulf of Mexico in SLA.

Similarly, the PPT response time to GWSA, TWSA, and SLA was also identified using lag–lead correlation. PPT lead TWSA by 5 months, whereas a PPT showed a lag of one month with GWSA and SLA. In addition, cross-correlation analysis was also carried out between GWSA and SLA, which was found to have a lag of 1 month between them.

3.2. Covariability and Correlation Analysis Using the SVD Technique

3.2.1. TWSA Variability in the Atlantic and the Pacific Ocean

The first mode heterogeneous correlation maps were represented for each hurricane month (June to November) between SST and TWSA lagged by 5 months. The modes corresponding to SCF greater than 10% were considered. The first mode of SVD was considered, as it explained the maximum variability. The maps depicting blue signify a negative correlation (cold phase), while red signifies a positive correlation (warm phase) above a 90% confidence level (Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10 and Figure 11). The maps also represent the significant regions with “+” signs for increasing and represent the regions with “-” signs for decreasing oceanic-atmospheric variabilities. Each significant region was given a number. The numbers before the “+” and “-” signs (for example, 1(+), 2(-), and 4(+)) represent the first, second, and fourth correlation group numbers, depicting the SST regions. The first region, 1(+) or 1(-), represents the highest correlated region, 2 (+) and 2(-) represent the slightly less correlated regions, and so on, depending on the correlation coefficient. The Squared Covariance Fraction (SCF) in the Pacific Ocean ranged from 83% in Nov-SST and April-TWSA to 93% in June-SST and Nov-TWSA. The NSC was found to range from 11% in Sept-SST and Feb-TWSA to 15% in June-SST and Nov-TWSA, as shown in

Table 2. Squared Covariance Fraction (SCF) and Normalized Squared Covariance (NSC) values corresponding to the SST relationship with Terrestrial Water Storage Anomaly (TWSA), GWSA, Sea Level Anomaly (SLA), and precipitation (PPT) (in %) of the Pacific and Atlantic Ocean for the HUC-03 region.

Climate Variability	SST–TWSA 1	SST–PPT 2	SST–GWSA 3	SST–SLA 4	5-Lag		2-Lag		1-Lag		1-Lag
					TWSA		PPT		GWSA		SLA
					SCF (%)	NSC (%)	SCF (%)	NSC (%)	SCF (%)	NSC (%)	SCF (%)	NSC (%)
Pacific	June–Nov	June–Aug	June–July	June–July	93	15	63	7	55	8	59	7
	July–Dec	July–Sept	July–Aug	July–Aug	92	14	62	6	56	9	40	9
	Aug–Jan	Aug–Oct	Aug–Sept	Aug–Sept	89	11	63	4	55	9	58	10
	Sept–Feb	Sept–Nov	Sept–Oct	Sept–Oct	85	11	76	5	55	11	76	12
	Oct–March	Oct–Dec	Oct–Nov	Oct–Nov	87	11	93	11	65	11	67	9
	Nov–April	Nov–Jan	Nov–Dec	Nov–Dec	83	11	45	5	59	11	72	11
Atlantic	June–Nov	June–Aug	June–July	June–July	92	12	70	6	40	8	50	7
	July–Dec	July–Sept	July–Aug	July–Aug	94	13	62	5	52	8	51	9
	Aug–Jan	Aug–Oct	Aug–Sept	Aug–Sept	88	9	70	6	57	8	57	9
	Sept–Feb	Sept–Nov	Sept–Oct	Sept–Oct	81	8	89	7	59	8	70	10
	Oct–March	Oct–Dec	Oct–Nov	Oct–Nov	90	12	89	10	66	10	60	9
	Nov–April	Nov–Jan	Nov–Dec	Nov–Dec	81	8	54	6	60	9	63	10

¹ June–Nov refers to the SST of June with a lag of 5 months in TWSA (i.e., Nov month) and so on. ² June–Aug refers to the SST of June with a lag of 2 months in PPT (i.e., Aug month) and so on. ³ June–July refers to the SST of June with a lag of 1 month in GWSA (i.e., July month) and so on. ⁴ June–July refers to the SST of June with a lag of 1 month in SLA (i.e., July month) and so on.

.

A distinct ENSO pattern in the Pacific Ocean was observed from June to July SST with a magnitude of correlation values ranging from −0.8 to +0.8, but the ENSO-like pattern reduced after August SST with a magnitude of correlation values ranging from −0.7 to +0.7 in different months. Thus, the ENSO teleconnection pattern was observed and revealed to have a strong correlation with TWSA in the Pacific region, as shown in Figure 3. Similar to the Pacific region, a strong positive SST region was found in the Atlantic Region from June and July in the northern Atlantic Region and northern South America. These regions are significant and positively correlated with TWSA (Figure 4).

From Aug to Nov cooling of SST was observed in almost all the SST regions while one warm SST region was observed in the Pacific Ocean. The climate variable teleconnection resulted in a decrease in TWSA in the HUC-03 region as shown in Figure 4. The SCF was found ranging from 81% in Sept-SST and Feb-TWSA to 94% in July-SST and Dec-TWSA in the Atlantic region. The NSC was found to range from 8% in Sept-SST and Feb-TWSA and Nov-SST and April-TWSA to 13% in July-SST and Dec-TWSA, as shown in Table 2. The Atlantic climatic variables indicated a negative and positive correlation with the HUC-03 region TWSA.

3.2.2. Precipitation Variability in the Atlantic and the Pacific Ocean

The Pacific warming SST regions were detected in all the months from June to Nov with high spatial coverage in Oct and July near tropical zone whereas one SST region was detected in the northern Atlantic region. Furthermore, in Oct, the west coast of the United States, the southern part of Mexico in the Pacific region and area in the Gulf of Mexico were observed to have a significant and positive correlation with PPT. The SCF from the SVD analysis among SST and PPT were ranging from 45% in Nov-SST and Jan-PPT to 93% in Oct-SST and Dec-PPT explaining the first mode of variability as shown in Table 2. The correlation values obtained for a lag of 2 months SST on PPT in the HUC-03 region were from −0.7 to +0.7. The climate variable response that influence PPT had a limited region in Nov, Sept, and Aug.

The decrease in the SST over the Pacific region was observed in Aug-SST and Oct-PPT over the HUC-03 region whereas in May to July-SST and Sept to Oct-SST at a lag of 2 months on PPT (i.e. July to Sept and Nov-Dec PPT respectively). The overall teleconnection showed mixed signals in different months in the HUC-03 region as shown in Figure 5. Comparing the SCF and NSC in the Atlantic Ocean, a higher value than the Pacific was observed. The SST variance of the first mode of SVD ranged from 54% to 89% PPT variability as shown in Table 2. A positive and significant correlation in two regions during Oct-SST and Dec-PPT across northern Atlantic and in the tropical region as shown in Figure 6. Moreover, PPT with a negative and significant correlation in the western Tropical region in July-SST and Sept-PPT. In general, it represents warming in June and Oct-SST with a lag of two months (Aug and Dec, respectively) in PPT signifying increase or decrease in PPT.

3.2.3. Groundwater Variability in the Atlantic and the Pacific Ocean

From Table 2, it can be noted that the range of SCF values are quite lower than the other variables, from 55% in June-SST and July-GWSA to 65% in Oct-SST and Nov-GWSA. A lag of 1-month SST, i.e., June, Sept, and Nov, with GWSA, i.e., July, Oct, and Dec: a large spatial extent extending from east to west equatorial Pacific in the central area was observed. A significant negative correlation in the respective regions was identified with a maximum of -0.8 correlation values (Figure 7). The heterogeneous correlation maps showed the decrease in SST is linked with a decrease in GWS, mainly in Florida. From the figure, most parts of Florida were negatively correlated with the SST regions. In other words, the cooling of SST influenced lower GWS in the southeast HUC-03 region. The west coast of the United States showed a teleconnection pattern with GWSA variability. In the Atlantic Region, Northeast US, and Northeast South America showed a significant and high correlation in July-SST and August-GWSA. Similar to the Pacific region, the SCF range was less, i.e., from 40% to 66% explained variability in GWSA by SST in the Atlantic region. Furthermore, in June-SST and July-GWSA, limited SST regions were only observed and may have only little influence. A strong positive correlation pattern in July-SST and Aug-GWSA was observed to increase GWSA due to an increase in SST whereas in other months a decrease in GWSA was observed due to a decrease in SST as illustrated in Figure 8.

3.2.4. Sea Level Variability in the Atlantic and the Pacific Ocean

In the Pacific region, a strong positive correlation is observed, extending from the Eastern Pacific to the Western Pacific and having major spatial coverage in the tropical regions. Throughout the study period, except Nov- and July-SST and Dec- and Aug-SLA, all other months and had a strong warm SST, depicting the warm ENSO phase on the west coast of the United States. A significant strong correlation value of more than +0.8 was identified in the SST regions. Small pockets of positive correlations were observed in July-SST and Aug-SLA, whereas strong negative correlations extending from west to east were identified in Nov-SST and Dec-SLA.

From the SVD analysis, a range from 40% to 76% of the variability was elucidated by the first mode characterizing the El Niño–SST pattern in the eastern equatorial Pacific and a negative correlation in only a few pockets, as shown in Figure 9i. A strong association of SST has been characterized by sea level stations in the East Coast, as shown in Figure 10i. It indicates an increase in sea level with an increase in SST and vice-versa. During the warm phase, the temperature rises that warms the ocean water lead to an increase in the sea level. Figure 9i shows the first mode of SVD of June-Nov SST and sea level anomaly of July-Dec in the Atlantic region. Significant strong correlations were evident in the northern part of South America and the northeast US, except in Oct–Nov, Nov–Dec, and Aug–Sept. Unlike the Pacific, the SCF ranges from 50% to 70%, suggesting sea levels are less closely linked to Atlantic SST. As observed, two to three regions are well correlated with summer months. The cooling SST phase is observed, exhibiting a close relationship to sea level stations over the US coast (Figure 10ii). Areas in this region have had a correlation coefficient greater than ±0.65.

3.2.5. Precipitation Relationship with Terrestrial Water and Groundwater Storage

The SVD analysis was carried out to evaluate the effect of PPT on terrestrial water and groundwater variability. The core area of the study area, i.e., Florida, experiences monsoon season from May to October. A significant positive and negative correlation between SST and TWSA was identified throughout the study period. However, the gridded TWSA showed no significant change due to PPT. This may be due to the higher influence of SST than PPT. The SCF of those months was 85% and 95%, explaining the maximum variability of the first mode, as shown in

Table 3. SCF and NSC values corresponding to the PPT relationship with TWSA, GWSA, and SLA (in %) for the HUC-03 region.

PPT–TWSA 1	PPT–GWSA 2	PPT–SLA 3	5-Lag		1-Lead		1-Lead
			TWSA		GWSA		SLA
			SCF (%)	NSC (%)	SCF (%)	NSC (%)	SCF (%)	NSC (%)
June–Nov	June–May	June–May	85	2	60	5	45	5
July–Dec	July–June	July–June	80	8	77	12	57	6
Aug–Jan	Aug–July	Aug–July	92	5	52	6	73	5
Sept–Feb	Sept–Aug	Sept–Aug	95	8	56	5	70	9
Oct–March	Oct–Sept	Oct–Sept	94	6	82	7	75	5
Nov–April	Nov–Oct	Nov–Oct	88	5	90	7	89	7

¹ June–Nov refers to the PPT of June with a lag of 5 months in TWSA (i.e., Nov month) and so on. ² June–May refers to the PPT of June with a lead of 1 month in GWSA (i.e., May month) and so on. ³ June–May refers to the PPT of June with a lead of 1 month in SLA (i.e., May month) and so on.

. The range resulted in a minimum of 80% in July-PPT and Dec-TWSA and a maximum of 95% in Sept-PPT and Feb-TWSA. A major component that governs GWS is PPT, the correlation values are positive in Aug-July whereas other period showed negative correlation values. The SCF in July-SST and Aug-GWSA was 56% in the Pacific region, which is greater than Aug-PPT and July-GWSA (52%). The NSC is also greater with SST (9%) than PPT (6%). However, in July–June (PPT–GWSA), the NSC is higher with PPT (12%), as shown in Table 3. In all other months, SST has a higher influence on GWSA than that of PPT. The SCF value ranged from 45% June-PPT and May-SLA to 89% in Nov-PPT, and SLA-Oct, respectively. The NSC value ranged from 5% in three months to 9% in Sept-PPT and Aug-SLA.

3.2.6. Relationship between Sea Level and Groundwater Storage within the Contiguous United States

The heterogeneous relationship between the sea level and GWS was also evaluated using the SVD technique. The highest SCF was found to be 79% in Nov–Dec and NSC 12% whereas the maximum NSC resulted in 15% in Sept–Oct. The variability of GWS resulted in SCF values of 71%, 67%, 72%, 76%, 60%, and 79% for the first mode for the result of June–July, July–Aug, Aug–Sept, Sept–Oct, Oct–Nov, and Nov–Dec, as shown in

Table 4. SCF and NSC values corresponding to the GWSA and SLA relationship (in %) for the HUC-03 region.

GWSA–SLA 1	1-Lag
	SCF (%)	NSC (%)
June–July	71	11
July–Aug	67	14
Aug–Sept	72	11
Sept–Oct	76	15
Oct–Nov	60	12
Nov–Dec	79	12

¹ June–July refers to the GWSA of June with a lag of 1 month in SLA (i.e., July month) and so on.

. Similarly, NSC values of 11%, 14%, 11%, 15%, 12%, and 12% were obtained. Figure 11 illustrates the spatial and temporal relationship between sea level record and GWS of the contiguous United States. It could be revealed that the East Coast GWS is positively correlated with SLA at the stations in all hurricane months. Most of the stations on the East Coast show a positive correlation with the GWS on the East Coast whereas a negative correlation with the area near California, Mexico, Arizona, and New Mexico. One of the stations in the tide gauge showed a negative correlation in all months whereas few SLA stations showed a non-significant relationship with the GWSA (represented by black dots in Figure 11). In August and September, the tide gauge showed that the increase in sea level record decreases GWS in the western United States. However, the East Coast sea level rise increases GWS. This could be due to other influencing factors, such as climate-induced variability, for instance, ENSO in the western US results in dry conditions, or also due to anthropogenic changes in those areas, such as groundwater pumping and irrigation. The Midwest did not show any significant relationship between sea level and groundwater in a few months.

4. Discussion

Physical climate variability was determined by using datasets such as TWS, sea level, GWS, and PPT. Climate variability has implications on the sea level, groundwater variability, and terrestrial water availability. In Figure 3, it was found that SST variabilities had a strong association with ENSO indices. However, it cannot capture PDO and AMO signals as GRACE has only 15 years of data. During the summer months, TWSA had high SCF and NSC values with SST (Table 2). From the SCF and NSC, it is noteworthy to conclude that SST had a higher influence than PPT, and close association with GWSA also influences sea level variability. Furthermore, a higher correlation was found between TWSA and SST, which may be attributed to the greater influence of SST than PPT in TWSA variability (Table 2 and Table 3). This confirms the discussion by Reager et al. [24] that the variation in land water storage for 12 years from 2002 to 2014 is attributed to climate-driven processes.

In addition, Engström and Waylen [54] found strong AMO and ENSO conditions that increase water shortage in the southern U.S. Hanson [55] investigated the relationship among climate variability, groundwater well levels, streamflow, tree-ring indices, and PPT in the southwestern region of the United States and found that ENSO and PDO-like components exhibited cyclic hydrologic variability. In addition, the Empirical Orthogonal Functions, computed from sea level reconstruction, showed a pattern associated with PDO. However, the study by Hamlington et al. [56] corroborated its association with common forcing rather than PDO and variability dynamics SVD analysis in this study also confirms the teleconnection of the Atlantic and Pacific Region (perhaps like ENSO and PDO) with GWS, sea level, and TWS, explaining the plausible increase in groundwater storage within HUC-03.

Sadeghi et al. [57] found a decrease in streamflow since 2000 due to five (2000, 2007, 2008, 2010, and 2011) La Niña events. A strong El Niño SST pattern with a large positive correlation extending from eastern to western Pacific was observed. Additionally, the warm (cold) ENSO phase showed a direct relationship with positive (negative) TWSA in the HUC-03 region with a lag of five months, as shown in Figure 3, and is also similar to what was observed by Linage et al. [58] for the Amazon Basin. The consistent warming followed by the cooling of SST in the ENSO region can be one of the most important phenomena that drive terrestrial water variability. During the warm phase, high temperature leads to higher evaporation that might cause a decrease in terrestrial water within the HUC-03 region. In Figure 3, warm SST regions were observed in June–July SST and Nov–Dec TWSA, whereas cold SST regions were observed in Aug–Nov SST and Jan–April TWSA within the HUC-03 region. This result supports that a warming Pacific SST is linked to decreasing inland water in the Pacific region.

In this study, TWS from GRACE was used to address the association of TWS with climatic and sea level variability. TWS variability has a substantial effect on sea level variability, as found by Hamlington et al. [56]. PDO further drives TWS variability that is interlinked with sea level changes. Monitoring the relationships among these variables is important as it has a significant role in the water cycle. The studies related to TWS, PPT, and sea level have been done in the past, but this study incorporates a less focused research domain, the oceanic and atmospheric association with groundwater. A close association between the detrended Global Mean Sea Level (GMSL) record from altimetry data and two ENSO events, El Niño and La Niña, is observed [59]. Furthermore, short-term variation in GMSL was associated with global land-water storage that declines throughout the El Niño event with the rise in GMSL, especially in tropical regions [60]. This can be attributed to high PPT in oceans and less PPT in the land throughout the El Niño event and is also supported by Gu and Adler [61]. In addition, La Niña leads to reduced GMSL [62].

Studies suggest that changes in wind forcing are the primary driver for decadal Pacific sea level variability [63] as well as Rossby waves leading to variations near the coastal area (https://oceanservice.noaa.gov/facts/rossby-wave.html). Warming of the ocean increases easterly trade winds, enhancing PPT and sea level rise. Winds are associated with climate modes, such as ENSO, PDO, and AMO. It is essential to account for the ocean circulation to evaluate the spatial relationship with sea level changes. Out of this study, a good correlation between sea level and SST, plus a close association with ENSO and PDO, were revealed. Sea level variability is majorly driven by the temperature above the ocean and atmospheric circulation [64].

Han et al. [65] emphasized understanding natural variability impacts on sea level variation so this study was also focused on identifying regions correlated with SST regions. Reager et al. [24] found a net sea level sink during 2002–2014 that could directly link with the increase in GWS or TWSA. This could further be supported by our study that as GWS increases, the sea level also increases in the East Coast; however, the subsurface catchment properties can act as a particularly strong filter on climate variability [66]. The effects of catchment size on climate–hydrology relationships, the degree of damping of climate signals related to hydrodynamic properties of the aquifer, as well as the influence of human stresses must be taken into consideration. Therefore, climate variability plays an influential role in GWS variability. GWS is further influenced by the changes in PPT driven by climate variability. In the western United States, the groundwater showed cold regions (as shown in Figure 11i); this can be due to ENSO-induced variability in the western states resulting in dry conditions during La Nina events. The IPCC has included thermal expansion and ice melting as a contributing factor of sea level change but has excluded the effects of groundwater variability on the sea level budget [67]. Furthermore, studies also found that the quantity of current GWS is comparable to a sea level differential of ~320 m to 330 ± 41 m. The change is significant if a certain amount of water contributes to sea level change [27,68]. In addition, Wada et al. [69] also found that the groundwater depletion contributes to sea level rise. Furthermore, Sahagian et al. [70] found that the groundwater mining in Arabian and African aquifers can have implications for sea level rise in the 21st century. In this study, higher NSC and SCF values between GWSA and SLA were observed to have higher correlation values. The correlation values also substantiate that the variability of the first mode primarily contributes to a close association between sea level and groundwater. From this, we can conclude that groundwater variability has a substantial effect on sea level variability. Associating the groundwater NSC during the study period with SST and PPT, it was evident that larger NSC values of SST indicate the influential role of SST on groundwater variability rather than PPT within the HUC-03 region. The SCF values of the first mode of June–Nov SST clearly shows that PPT and groundwater is the main regulator of hydroclimatic co-variability in the basin, especially in October. Almanaseer and Sankarasubramanian [71] found that the July to September groundwater level is not statistically significant with Nino 4.5 in the Southeast United States. Furthermore, they found groundwater and streamflow were found to be the dominant factors driving hydroclimatic variability in winter.

The current study examined the response of GWS driven by PPT and SST anomalies during hurricane months over the HUC-03 region. The hydroclimatic co-variability during the fall season (Sept–Nov) was found to be stronger than the summer season (June–Aug), as indicated by the first mode of high SCF values in Atlantic SST–GWSA (Table 2). During the fall season, El Niño events are at their peak, and were found to be an important oscillation affecting hydroclimatic variability. The less variability during summer primarily may be because of high SSTs. A high SST drives jet streams and creates a pressure difference that can change the PPT pattern. In addition, the southeast region is linked with the East Asian Jet stream that leads to increased PPT in the western US [28]. Furthermore, the PPT pattern was found to resemble El Niño events and the study by Hamlington et al. [56] corroborated ocean–atmosphere interactions are related to decadal persistence. Moreover, it emphasized that increased PPT is associated with ENSO, forming a movement of the water cycle from land to ocean. In this study, during La Niña (El Niño), the PPT value tends to increase (decrease). The positive relationship between PPT and TWSA leads to a decrease in terrestrial storage due to a decrease in PPT. The NSC value of PPT–TWSA is low when compared with SST–TWSA, suggesting SST has higher influence in the October and November months. In the Southeastern United States, climate variability tends to have a higher influence on TWSA than PPT from previous studies and is also evident in the current research. The weak PPT season may have a strong influence on climatic effects, further influencing the groundwater variability. The GWSA–PPT relationship was found to be non-significant at a 90% significance level in all the study period. This could be inferred to have a non-linear relationship between them. Non-linearity may prevail uncertainties in climate-related forecast studies. However, SVD is a linear method. Non-linearity may lead to altering short-term and long-term variability, which needs to be studied in the future. In addition, this study does not include anthropogenic effects such as groundwater pumping, withdrawal, irrigation, and cropping patterns. These could also affect sea level. However, a study by Thomas and Famiglietti [72] showed that groundwater use does not have much influence on HUC-03 groundwater storage. A further detailed study could be done for additional understanding of groundwater use and the effects of climate variability. The long-term and short-term variabilities in streamflow, precipitation, and temperature have been done but not yet studied for groundwater. From this, we can conclude that SST has a large influence on the GWSA in the HUC-03 region. Moreover, during La Niña (El Niño), the GWS value tends to increase (decrease) and the negative relationship between PPT and GWSA leads to decreased GWS due to a decrease in PPT. However, in the HUC-03 region, the SST is dominant in June, May, and October. The spatiotemporal correlation maps showed positive and negative values in all the study periods.

Future work will require extended work on longer GWS time series data after it becomes available in the future, which can give a better vision on the effect of the sea level variability record. This will provide an additional understanding of the role of high- and low-frequency teleconnections. The outcome of this study can further be utilized to predict the groundwater using climate forecasts over the study region.

5. Conclusions

The current study focused on the contribution of climate to influence groundwater variability within the HUC-03 utilizing the SVD technique. The analysis was based on monthly PPT and GWS data to comprehend groundwater dynamics driven by PPT and climate processes. The monthly gridded terrestrial water storage obtained from GRACE was correlated with the preceding 5 months of PPT and SST. The highest spatial average of the SST forcing was found to lead PPT by 2 months and SST showed a strong response on GWSA with a lead by 1 month.

The ENSO and PDO are associated with the variability in the Pacific and the Atlantic Ocean, and through the use of correlation and temporal analyses of the first mode of SVD, their signals were typically reflected in GWS, TWS, PPT, and sea level within the HUC-03 region. However, only a few PDO patterns were recognized in the groundwater of the HUC-03 region. The available GRACE dataset was only approximately 1.5 decades so the SST regions, including PDO and AMO, were not well-captured in this study. The SVD analysis noticeably shows that climate variability is the leading variable that influences the hydroclimatic variables in the HUC-03 region. Moreover, regions with high NSC, indicating TWS, GWS, and sea level, are the primary drivers of the hydroclimatic processes. In addition, the SCF during the summer season is lower than the fall season, signifying high temperature decreases the recharge. From this study, the spatial–temporal relationship between sea level and the response of GWS can be apprehended. The responses differ according to location, where the Southeast Coast had groundwater changes with sea level change.

The main conclusions drawn in this article are as follows.

• The precipitation was found to have less influence on groundwater variability than SST.
• The groundwater variability was found to have a significant relationship with the Pacific and Atlantic Ocean variability.
• The GWS in the east, south, and northwest coast of the United States were found to have a positive relationship with sea level variability whereas a negative relationships between GWS and SLA were observed in the western and southwestern United States, near California, Mexico, Arizona, and New Mexico.

Prominent relationships between SST and groundwater may not be linear as assumed in the current study. Moreover, the lagged relationship showed the potential of predicting groundwater prior to the HUC-03 region. By and large, the key contributions of the current study are:

(1)

Comprehensive analysis of the lag and the lead relationship between SST, PPT with groundwater, and sea level anomaly.

(2)

Evaluating the effect of climatic variability on groundwater and sea level.

(3)

Identification of the influencing variable that drives sea level and groundwater variability.

(4)

Categorizing SST regions that drive the groundwater and sea level in different parts of the study region.

(5)

Evaluating the robustness of SVD to assess the relationship between climate variability and groundwater.

In the current study, the hydroclimatic co-variability identification helped to infer the climate variability implications on water resources. Although new regions were obtained, climate change-induced hydrological cycle complexities may have some uncertainties. A similar study can also be done in other HUC regions. Furthermore, the teleconnected regions could be utilized for forecasting. Future work may use a longer time period as it becomes available for wider insight. In addition, groundwater forecasting model development could also be a potential future research area in HUC regions.