Mapping Aquifer Storage Properties Using S-Wave Velocity and InSAR-Derived Surface Displacement in the Kumamoto Area, Southwest Japan

: We present a novel approach to mapping the storage coefﬁcient ( Sk ) from InSAR-derived surface deformation and S-wave velocity ( Vs ). We ﬁrst constructed a 3D Vs model in the Kumamoto area, southwest Japan, by applying 3D empirical Bayesian kriging to the 1D Vs proﬁles estimated by the surface-wave analysis at 676 measured points. We also used the time series of InSAR deformation and groundwater-level data at 13 well sites covering April 2016 and December 2018 and estimated the Sk of the conﬁned aquifer. The Sk estimated from InSAR, and well data ranged from ~0.03 to 2 × 10 − 3 , with an average of 7.23 × 10 − 3 , values typical for semi-conﬁned and conﬁned conditions. We found a clear relationship between the Sk and Vs at well locations, indicating that the compressibility of an aquifer is related to the stiffness or Vs . By applying the relationship to the 3D Vs model, we succeeded in mapping the Sk in an extensive area. Furthermore, the estimated Sk distribution correlates well with the hydrogeological setting: semi-conﬁned conditions are predicted in the Kumamoto alluvial plain with a high Sk . Our approach is thus effective for estimating aquifer storage properties from Vs , even where limited groundwater-level data are available. Furthermore, we can estimate groundwater-level variation from the geodetic data.


Introduction
Groundwater is a vital water resource in arid and remote regions. Given changing climate and human development, there is an increasing need for groundwater exploration, monitoring, and management. Effective groundwater management commonly relies on measurements of hydraulic head and storage properties. Groundwater from deep confined aquifers offers a viable source of fresh water and a buffer against severe drought conditions. Estimating hydrogeological properties is therefore important for developing sustainable, long-term water management strategies. Aquifer characteristics are usually estimated by analyzing pumping or recovery test data, including measuring water-level variations at monitoring wells [1][2][3][4][5][6]. However, the high cost of aquifer testing has restricted the monitoring of groundwater resources, and knowledge of their spatiotemporal evolution remains sparse [7]. their spatial distribution in the Kumamoto area has not been fully investigated because limited groundwater-level data were available.

Study Area
The Kumamoto area lies in central Kyushu Island, southwestern Japan ( Figure 1). Hossain et al. [43,44] and Hosono et al. [45] reported its geological and hydrogeological settings and the geochemical processes controlling its groundwater. The eastern part of the study area is a volcanic area of the Quaternary age covered by pyroclastic-flow deposits erupted from the Aso volcano and alluvial deposits of the Pleistocene age. The western portion (especially coastal area) comprises the Kumamoto alluvial plain, consisting of Holocene Ariake clay [45,46]. The pyroclastic-flow deposits erupted from Mt. Aso over four periods between 90 and 30 ka and are classified into four types: Aso-4, Aso-3, Aso-2, and Aso-1. The pyroxene andesite called Togawa lava lies on the Aso-1, and the upper clacked part connects to the bottom of Aso-3 [46]. Two main faults are distributed in the study area: the Futagawa fault, which cuts the lava plateau and runs along the boundary between Cretaceous rocks and the Kumamoto plain, and the Hinagu fault, which runs from the south next to the alluvial plain toward the north through bedrocks to merge with the Futagawa fault. The 2016 Kumamoto earthquake sequence (Mw 7.0) along these two fault systems [47][48][49] caused widespread damage and disrupted infrastructure (e.g., the work of [50]).
Water supply in the Kumamoto area mainly relies on its abundant groundwater [51]. The Kumamoto aquifer system is mainly recharged from the western rim of Mount Aso, three major highland areas (Ueki, Kikuchi, and Takayubaru), and the midstream area of the Shira River being located in the north next to Takayubaru; its topography leads the groundwater flow direction to the southwest [52]. The Aso volcanic rock (pyroclastic flow and Togawa lava) is highly porous and permeable, making it an excellent aquifer ( Figure 2) [53]. The unconfined aquifer is composed of recent pyroclastic-flow deposits (uppermost Aso-4) and partial marine sediments (<50 m deep). The underlying confined/semi-confined aquifer consists of older pyroclastic deposits and volcanic flow lavas (60-250 m) [45,54].

of 17
Sk, their spatial distribution in the Kumamoto area has not been fully investigated because limited groundwater-level data were available.

Study Area
The Kumamoto area lies in central Kyushu Island, southwestern Japan ( Figure 1). Hossain et al. [43,44] and Hosono et al. [45] reported its geological and hydrogeological settings and the geochemical processes controlling its groundwater. The eastern part of the study area is a volcanic area of the Quaternary age covered by pyroclastic-flow deposits erupted from the Aso volcano and alluvial deposits of the Pleistocene age. The western portion (especially coastal area) comprises the Kumamoto alluvial plain, consisting of Holocene Ariake clay [45,46]. The pyroclastic-flow deposits erupted from Mt. Aso over four periods between 90 and 30 ka and are classified into four types: Aso-4, Aso-3, Aso-2, and Aso-1. The pyroxene andesite called Togawa lava lies on the Aso-1, and the upper clacked part connects to the bottom of Aso-3 [46]. Two main faults are distributed in the study area: the Futagawa fault, which cuts the lava plateau and runs along the boundary between Cretaceous rocks and the Kumamoto plain, and the Hinagu fault, which runs from the south next to the alluvial plain toward the north through bedrocks to merge with the Futagawa fault. The 2016 Kumamoto earthquake sequence (Mw 7.0) along these two fault systems [47][48][49] caused widespread damage and disrupted infrastructure (e.g., the work of [50]).
Water supply in the Kumamoto area mainly relies on its abundant groundwater [51]. The Kumamoto aquifer system is mainly recharged from the western rim of Mount Aso, three major highland areas (Ueki, Kikuchi, and Takayubaru), and the midstream area of the Shira River being located in the north next to Takayubaru; its topography leads the groundwater flow direction to the southwest [52]. The Aso volcanic rock (pyroclastic flow and Togawa lava) is highly porous and permeable, making it an excellent aquifer ( Figure  2) [53]. The unconfined aquifer is composed of recent pyroclastic-flow deposits (uppermost Aso-4) and partial marine sediments (<50 m deep). The underlying confined/semiconfined aquifer consists of older pyroclastic deposits and volcanic flow lavas (60-250 m) [45,54].   [55]. (b) Location map of the study area (red rectangle) created using ArcGIS ® software by Esri. The red rectangle indicates the location of panel a.  [46,52]). The location of cross-section lines is in Figure 1.

InSAR Data
We used InSAR surface displacement data reported by Ishitsuka et al. [42], which was derived from persistent scatterer interferometry (PSI) analysis of ALOS-2/PALSAR-2 synthetic aperture radar (SAR) images. PSI analysis estimates time-series surface displacements using phase differences at coherent pixels called persistent scatters (PSs) [56,57]. In addition, PSI analysis has strategies to mitigate nuisance components of the phase differences attributed to changes in atmospheric states and errors in a digital elevation model. The accurate phase information at PSs and the mitigation of nuisance phase differences (phase differences caused by factors other than surface displacement) enable us to estimate accurate time-series surface displacements.
Ishitsuka et al. [42] used 28 ALOS-2/PALSAR-2 images acquired in a descending track mode orbit between 18 April 2016 and 10 December 2018. The PSI processing flow used in Ishitsuka et al. [42] was based on Kampes [57]. First, they selected interferometric pairs for single reference images. To select PS candidates (PSCs), they used an amplitude dispersion index [56] of 0.30. Following that, differential interferograms were created at the PSCs with a single-look. Next, they removed topographic phase components using a 10 m mesh external digital elevation model (DEM) provided by the Geospatial Information Authority of Japan. The stability of the interferometric phase of PSCs was then assessed based on phase coherence [56] and selected PSs as pixels above a coherence threshold of 0.80. Subsequently, the residual orbital fringes were removed, and the DEM errors and the displacement rates were estimated from a least-squares method. Phase unwrapping was then performed using a minimum cost flow algorithm [58]. Next, they reduced the atmospheric phase component by smoothing the change in the temporal phase because the atmospheric phase contribution is generally characterized by a high temporal frequency (e.g., the work of [56]). They showed that the line-of-sight surface displacements estimated from the PSI analysis are consistent with global navigation satellite system (GNSS) observations in the study area [42]. For characterizing the source of surface  [46,52]). The location of cross-section lines is in Figure 1.

InSAR Data
We used InSAR surface displacement data reported by Ishitsuka et al. [42], which was derived from persistent scatterer interferometry (PSI) analysis of ALOS-2/PALSAR-2 synthetic aperture radar (SAR) images. PSI analysis estimates time-series surface displacements using phase differences at coherent pixels called persistent scatters (PSs) [56,57]. In addition, PSI analysis has strategies to mitigate nuisance components of the phase differences attributed to changes in atmospheric states and errors in a digital elevation model. The accurate phase information at PSs and the mitigation of nuisance phase differences (phase differences caused by factors other than surface displacement) enable us to estimate accurate time-series surface displacements.
Ishitsuka et al. [42] used 28 ALOS-2/PALSAR-2 images acquired in a descending track mode orbit between 18 April 2016 and 10 December 2018. The PSI processing flow used in Ishitsuka et al. [42] was based on Kampes [57]. First, they selected interferometric pairs for single reference images. To select PS candidates (PSCs), they used an amplitude dispersion index [56] of 0.30. Following that, differential interferograms were created at the PSCs with a single-look. Next, they removed topographic phase components using a 10 m mesh external digital elevation model (DEM) provided by the Geospatial Information Authority of Japan. The stability of the interferometric phase of PSCs was then assessed based on phase coherence [56] and selected PSs as pixels above a coherence threshold of 0.80. Subsequently, the residual orbital fringes were removed, and the DEM errors and the displacement rates were estimated from a least-squares method. Phase unwrapping was then performed using a minimum cost flow algorithm [58]. Next, they reduced the atmospheric phase component by smoothing the change in the temporal phase because the atmospheric phase contribution is generally characterized by a high temporal frequency (e.g., the work of [56]). They showed that the line-of-sight surface displacements estimated from the PSI analysis are consistent with global navigation satellite system (GNSS) observations in the study area [42]. For characterizing the source of surface displacement, Ishitsuka et al. [42] created a model of temporal displacement caused by seasonal and long-term factors using sinusoidal and exponential functions, respectively, as follows: where α is the amplitude of the seasonal displacement, and θ represents the phase shift that indicates the beginning of the seasonal displacement. The factor sin(θ) was used to ensure that U(0) = 0. Here, β is a coefficient representing the magnitude of the exponential function, and k is the exponential decay coefficient. From the least-squares method (Equation (1)), four unknown parameters (i.e., α, θ, β, and k) were estimated. To estimate θ and k, a grid search algorithm was used with intervals of 5 • and a search range from −180 to 180 • for θ and an interval of 10 with search range from 20 to 10,000 for k for determining the coefficients α and β. The least-squares method was applied during the grid search. After that, the optimal values of four parameters (α, θ, β, and k) that reduce the error of mean square between the modeled and PSI time-series displacement were determined. The estimated displacement maps show transient and seasonal surface displacements posterior to the 2016 Kumamoto earthquake ( Figure 3). Transient displacements around the central and western part of the study area are attributed to groundwater migration through new coseismic ruptures (black rectangle in Figure 3a) and sediment compaction (the coastal area around the Ariake Sea (red rectangle in Figure 3a). Meanwhile, seasonal surface displacements based on Equation (1) in the northern and central parts are likely related to changes in groundwater levels, as shown in Figure 3b

Three-dimensional S-Wave Velocity Model
To estimate the spatial variation of S-wave velocities in the study area, we used 1D Vs profiles estimated from microtremors using a miniature seismic array [59] at 676 observation points by the National Research Institute for Earth Science and Disaster Resilience (NIED; [60]). The interval between observation points ranged from approximately 100 m to 1000 m. To collect microtremor data, we conducted microtremor surveys combining a triangular array (radius: 0.6 m) with one central station and an irregular triangu-

Piezometric Data
We used the time series of monthly average groundwater levels at 13 well sites (triangles in Figure 3b) acquired between 17 April 2016 and 10 December 2018. All wells were located in the northern part of the study area, where seasonal surface displacements were estimated with a magnitude up to 5 mm (Figure 3b; Ishitsuka et al. [42]), and measured the groundwater level at the confined aquifer. The upper and lower ends of the strainer depth for each well are shown in Table 1. For most wells, the strainer depth ranges from 50 to 170 m (Table 1), corresponding to a single aquifer. Surface displacements (topographic changes) were induced up to a few meters by the 2016 Kumamoto earthquake. Each well top height was re-measured after the earthquake (2017) by the leveling, and post-earthquake groundwater level (m a.s.l.) was adjusted by updated height data.

Three-Dimensional S-Wave Velocity Model
To estimate the spatial variation of S-wave velocities in the study area, we used 1D Vs profiles estimated from microtremors using a miniature seismic array [59] at 676 observation points by the National Research Institute for Earth Science and Disaster Resilience (NIED; [60]). The interval between observation points ranged from approximately 100 m to 1000 m. To collect microtremor data, we conducted microtremor surveys combining a triangular array (radius: 0.6 m) with one central station and an irregular triangular array (side ranges from 4 to 10 m or more) at each observation point with a sampling frequency of 200 Hz (Figure 4). The microtremor surveys were conducted using an integrated, portable housing seismometer, JU410, manufactured by the Hakusan Corporation based on the cloud microtremor observation system [61].
First, the horizontal/vertical (H/V) spectral ratio and the phase velocity dispersion curves of Rayleigh waves were estimated from the observed microtremors. Second, the Vs structure at each observation point was estimated by the joint inversion of the H/V spectral ratio and dispersion curves using a genetic algorithm inversion [62,63]. Finally, geostatistical interpolation by Empirical Bayesian Kriging 3D (EBK3D) was applied to predict the Vs at unsampled points in 3D space. EBK3D is available in ArcGIS Pro software as a geoprocessing tool. Empirical Bayesian kriging (EBK) [64] is different from other kriging methods; the EBK-based approach considers standard prediction errors and automatically estimates the semivariogram parameters using restricted maximum likelihood. EBK involves the following steps: 1. Estimation of a semivariogram model from the original data; 2. Prediction of data at each of the observed points from the semivariogram; 3. Estimation of a new semivariogram model and its weight from the predicted data; 4. Repeating steps 2 and 3 creates a spectrum of the semivariogram models; 5. Prediction of Vs and their standard errors at unmeasured locations using these weights.
Remote Sens. 2021, 13, x FOR PEER REVIEW 7 of 17 1-Estimation of a semivariogram model from the original data; 2-Prediction of data at each of the observed points from the semivariogram; 3-Estimation of a new semivariogram model and its weight from the predicted data; 4-Repeating steps 2 and 3 creates a spectrum of the semivariogram models; 5-Prediction of Vs and their standard errors at unmeasured locations using these weights.

Methodology: Mapping the Skeletal Storage Coefficient
The skeletal storage coefficient, Sk, of a confined aquifer (also called storativity) is a dimensionless measure of the volume of water released from or taken into storage per unit surface area of the aquifer per unit change in the hydraulic head [65]. It is defined as Sk = Ss•b, where Ss is specific storage (L −1 ), and b is the thickness (L) of the aquifer. We estimated the Sk of the confined aquifer system from 2016 to 2018 using the following relationship with surface displacement (Δb, from PSI) and the change in the hydraulic head (Δh) at the 13 well sites [66,67]: This relationship is valid when the deformation is elastic. Theoretically, the coefficient is described by the function of aquifer compressibility (Equation (3)), assuming that the porosity does not change significantly with elastic deformation [68]: where is the pore fluid density, is the acceleration due to gravity, is the matrix compressibility, n is the porosity, is the pore fluid compressibility, and b is the thickness (L) of the aquifer. Thus, the soil's (matrix) compressibility represents its volume change in response to the applied pressure [68].
Before estimating Sk using Equation (2), we eliminated multivariate outliers in the data using Z-score. It is a statistical method to measure the divergence of observed data points from its mean. Z-score can be computed as follows:

Methodology: Mapping the Skeletal Storage Coefficient
The skeletal storage coefficient, Sk, of a confined aquifer (also called storativity) is a dimensionless measure of the volume of water released from or taken into storage per unit surface area of the aquifer per unit change in the hydraulic head [65]. It is defined as Sk = Ss·b, where Ss is specific storage (L −1 ), and b is the thickness (L) of the aquifer. We estimated the Sk of the confined aquifer system from 2016 to 2018 using the following relationship with surface displacement (∆b, from PSI) and the change in the hydraulic head (∆h) at the 13 well sites [66,67]: This relationship is valid when the deformation is elastic. Theoretically, the coefficient is described by the function of aquifer compressibility (Equation (3)), assuming that the porosity does not change significantly with elastic deformation [68]: where ρ is the pore fluid density, g is the acceleration due to gravity, c m is the matrix compressibility, n is the porosity, c w is the pore fluid compressibility, and b is the thickness (L) of the aquifer. Thus, the soil's (matrix) compressibility represents its volume change in response to the applied pressure [68]. Before estimating Sk using Equation (2), we eliminated multivariate outliers in the data using Z-score. It is a statistical method to measure the divergence of observed data points from its mean. Z-score can be computed as follows: where X i is i th sample point, µ refers to the mean of the data set, and σ is the standard deviation. Since we analyzed multivariate data (i.e., surface deformation and hydraulic head), we computed Z-score in 2D space. We define outliers if the Z-score of the data exceeds the threshold ( Figure 5).
between Vs and aquifer compressibility. For example, Cha et al. [38] noted that shear stiffness and compressibility could be estimated using Vs, as they reported relationships between soil compressibility and the small strain parameters used in velocity-stress power relations. These relationships suggested that cemented dense sediments are characterized by low compressibility and high Vs and vice versa. Martin et al. [39] investigated a shear zone formed by an igneous rock to estimate hydraulic conductivity and normal stiffness using borehole data. They reported that low stress and stiffness regions exhibit high storage coefficients and hydraulic conductivity. Li et al. [69] derived a formula for rock compressibility based on soil mechanics concepts [70], indicating the dependence of rock compressibility on the rigidity of rock minerals. Therefore, the Sk can be linked to Vs. We thus constructed the Sk-Vs relationship by plotting the Sk and Vs at the water level depth estimated at each well. Using the Sk-Vs relationship, we then estimated Sk values from only the Vs. This approach enables us to map the Sk distribution where only Vs data were available.  Figure 6 shows maps of Vs at different depths. A low-velocity zone observed in the western coastal area correlated well with the Kumamoto alluvial plain covered by Holocene Ariake clay. We calculated the average Vs for a depth down to 30 m (Vs30) and could identify the faulted/fractured zones as low-Vs zones in the map of Vs30 (Figure 7) that coincided with the Futagawa fault system. We compared surface displacements by InSAR data analysis with the Vs structure to examine whether surface displacements are linked to subsurface properties. The low-velocity regions correlated with surface displacement associated with the 2016 Kumamoto earthquake sequence found by the InSAR technique [41,42,47]. The rock strength is reduced by faulting and fracturing, and such reduction is To predict the Sk values where groundwater-level data (i.e., well data) were not available, we proposed a new approach based on Vs. Previous studies show a relationship between Vs and aquifer compressibility. For example, Cha et al. [38] noted that shear stiffness and compressibility could be estimated using Vs, as they reported relationships between soil compressibility and the small strain parameters used in velocity-stress power relations. These relationships suggested that cemented dense sediments are characterized by low compressibility and high Vs and vice versa. Martin et al. [39] investigated a shear zone formed by an igneous rock to estimate hydraulic conductivity and normal stiffness using borehole data. They reported that low stress and stiffness regions exhibit high storage coefficients and hydraulic conductivity. Li et al. [69] derived a formula for rock compressibility based on soil mechanics concepts [70], indicating the dependence of rock compressibility on the rigidity of rock minerals. Therefore, the Sk can be linked to Vs. We thus constructed the Sk-Vs relationship by plotting the Sk and Vs at the water level depth estimated at each well. Using the Sk-Vs relationship, we then estimated Sk values from only the Vs. This approach enables us to map the Sk distribution where only Vs data were available. Figure 6 shows maps of Vs at different depths. A low-velocity zone observed in the western coastal area correlated well with the Kumamoto alluvial plain covered by Holocene Ariake clay. We calculated the average Vs for a depth down to 30 m (Vs30) and could identify the faulted/fractured zones as low-Vs zones in the map of Vs30 (Figure 7) that coincided with the Futagawa fault system. We compared surface displacements by InSAR data analysis with the Vs structure to examine whether surface displacements are linked to subsurface properties. The low-velocity regions correlated with surface displacement associated with the 2016 Kumamoto earthquake sequence found by the InSAR technique [41,42,47]. The rock strength is reduced by faulting and fracturing, and such reduction is reflected by a reduction in shear coupling and Vs [71]. Therefore, we demonstrate the potential of using microtremors to detect faulted/fractured zones.

Three-Dimensional S-Wave Velocity
Based on the Vs model ( Figure 6), furthermore, we can evaluate the liquefaction area. During liquefaction, the effective stress is close to zero because of the increase in pore pressure. A mechanism of undrained consolidation for soil liquefaction has often been discussed [72,73] based on the theory of poroelasticity [16] that relates effective stress to the pressure of the pore fluids. A reduction in effective stress associated with increasing pore pressure should be reflected by reducing the mechanical strength of sediments that weakens them, so the Vs is a very important parameter to evaluate liquefaction. The soil in the southwestern Kumamoto area was classified as soft or stiff in the Vs30 map according to the classification by the National Earthquake Hazard Reduction Program (NEHRP)and the American Society of Civil Engineers (ASCE) [74,75] (Figure 7). Velocities down to a depth of around 20 m were lower than 200 m/s in the western part (Figure 6), representing an upper limit of the range of velocities in liquefiable soil [76]. The region of low velocity is consistent with the area affected by liquefaction during the 2016 Kumamoto earthquake.
Remote Sens. 2021, 13, x FOR PEER REVIEW 9 of 17 reflected by a reduction in shear coupling and Vs [71]. Therefore, we demonstrate the potential of using microtremors to detect faulted/fractured zones.
Based on the Vs model ( Figure 6), furthermore, we can evaluate the liquefaction area. During liquefaction, the effective stress is close to zero because of the increase in pore pressure. A mechanism of undrained consolidation for soil liquefaction has often been discussed [72,73] based on the theory of poroelasticity [16] that relates effective stress to the pressure of the pore fluids. A reduction in effective stress associated with increasing pore pressure should be reflected by reducing the mechanical strength of sediments that weakens them, so the Vs is a very important parameter to evaluate liquefaction. The soil in the southwestern Kumamoto area was classified as soft or stiff in the Vs30 map according to the classification by the National Earthquake Hazard Reduction Program (NEHRP)and the American Society of Civil Engineers (ASCE) [74,75] (Figure 7). Velocities down to a depth of around 20 m were lower than 200 m/s in the western part (Figure 6), representing an upper limit of the range of velocities in liquefiable soil [76]. The region of low velocity is consistent with the area affected by liquefaction during the 2016 Kumamoto earthquake.  [47], and the green and blue lines refer to Futagawa and Hinagu faults [77].  [47], and th and blue lines refer to Futagawa and Hinagu faults [77]. Background map is from ArcGIS ® s by Esri.

Skeletal Storage Coefficient and S-Waves Velocity Relationship
We used the time series of InSAR-derived surface deformation and ground level data to determine the storage properties of the confined aquifer at 13 well sit ure 3b). Table 1 shows that the storage coefficient values at the 13 sites ranged from to 2 × 10 −3 , with an average of 7.23 × 10 −3 after outlier removal using Z-score. Mos correlations between InSAR displacement and well data are statistically significan on the p-value (<0.05). Since the correlation at SS-17 does not show statistical signi we excluded the estimated Sk in the following analysis. Pearson's correlation coe (R) and the standard errors of the Sk estimates are also shown in Table 1. The est values of Sk correspond to those of semi-confined and confined conditions [78,7 deep aquifer in the study area is semi-confined to confined due to the discontinuit impermeable clay layer of lacustrine sediments [45].
We then observed a negative correlation between the Sk and Vs (Figure 8). Ex tial curve fitting found the following empirical relationship between Sk and Vs: Around 83.3% of data points lie within the ±95% confidence interval of the re ship in Equation (5). The relationship between Sk and Vs can be explained by the re ship between compressibility and stiffness (or Vs), as mentioned in the previous se  [47], and the green and blue lines refer to Futagawa and Hinagu faults [77]. Background map is from ArcGIS ® software by Esri.

Skeletal Storage Coefficient and S-Waves Velocity Relationship
We used the time series of InSAR-derived surface deformation and groundwater-level data to determine the storage properties of the confined aquifer at 13 well sites (Figure 3b). Table 1 shows that the storage coefficient values at the 13 sites ranged from~0.03 to 2 × 10 −3 , with an average of 7.23 × 10 −3 after outlier removal using Z-score. Most of the correlations between InSAR displacement and well data are statistically significant based on the p-value (<0.05). Since the correlation at SS-17 does not show statistical significance, we excluded the estimated Sk in the following analysis. Pearson's correlation coefficient (R) and the standard errors of the Sk estimates are also shown in Table 1. The estimated values of Sk correspond to those of semi-confined and confined conditions [78,79]. The deep aquifer in the study area is semi-confined to confined due to the discontinuity of the impermeable clay layer of lacustrine sediments [45].
We then observed a negative correlation between the Sk and Vs (Figure 8). Exponential curve fitting found the following empirical relationship between Sk and Vs: Around 83.3% of data points lie within the ±95% confidence interval of the relationship in Equation (5). The relationship between Sk and Vs can be explained by the relationship between compressibility and stiffness (or Vs), as mentioned in the previous sections.

Mapping of the Skeletal Storage Coefficient to Monitor Groundwater Level from InSAR Data
To map Sk, we applied the Sk-Vs relationship to areas where only Vs were estimated from the microtremor survey ( Figure 9). The estimated Sk (Figure 9) is correlated with the hydrogeological setting, where the Sk greater than 0.01 in the southwestern part of the study area reflected semi-confined conditions due to discontinuity of the impermeable clay layer of lacustrine sediments underlying the deep aquifer [45]. This area is covered by an alluvial plain and terrace deposits ( Figure 1) composed of large amounts of unconsolidated sediments with high compressibility and low Vs. Wilson and Wöhling [80] found that vertical transmissivity and storage coefficient decrease significantly with depth due to the increasing tortuosity of permeable flow pathways in the deeper parts causing groundwater to flow horizontally through shallow pathways.

Mapping of the Skeletal Storage Coefficient to Monitor Groundwater Level from InSAR Data
To map Sk, we applied the Sk-Vs relationship to areas where only Vs were estimated from the microtremor survey ( Figure 9). The estimated Sk (Figure 9) is correlated with the hydrogeological setting, where the Sk greater than 0.01 in the southwestern part of the study area reflected semi-confined conditions due to discontinuity of the impermeable clay layer of lacustrine sediments underlying the deep aquifer [45]. This area is covered by an alluvial plain and terrace deposits ( Figure 1) composed of large amounts of unconsolidated sediments with high compressibility and low Vs. Wilson and Wöhling [80] found that vertical transmissivity and storage coefficient decrease significantly with depth due to the increasing tortuosity of permeable flow pathways in the deeper parts causing groundwater to flow horizontally through shallow pathways.

Discussion
This study demonstrates the feasibility of using Vs, InSAR-derived surface deformation, and groundwater-level data to map the storage properties of a confined aquifer at high spatial resolution (Figure 9). Therefore, we can estimate the spatiotemporal variation of groundwater level based on geodetic data (i.e., surface displacement derived from InSAR or GNSS). Our approach could offer great potential for improving the groundwater flow modeling by using the estimated Sk. However, the successful application of the approach requires consideration of the investigation depth of Vs estimated from microtremor surveys. As the Vs at a water level depth is used to construct a relationship between Sk and Vs, microtremor surveys need to be designed to cover the depth of the water level in a study area.

Discussion
This study demonstrates the feasibility of using Vs, InSAR-derived surface deformation, and groundwater-level data to map the storage properties of a confined aquifer at high spatial resolution ( Figure 9). Therefore, we can estimate the spatiotemporal variation of groundwater level based on geodetic data (i.e., surface displacement derived from InSAR or GNSS). Our approach could offer great potential for improving the groundwater flow modeling by using the estimated Sk. However, the successful application of the approach requires consideration of the investigation depth of Vs estimated from microtremor surveys. As the Vs at a water level depth is used to construct a relationship between Sk and Vs, microtremor surveys need to be designed to cover the depth of the water level in a study area.
In general, the uncertainty of Sk occurs for the following reasons.
(1) The Sk is biased if the horizontal displacement is not considered in the Sk calculation [81]. In our study area, horizontal fluid movement occurs over the groundwater system [23]. (2) Monitoring wells observe water from one aquifer in the saturated confined aquifer system [81]. Hoffmann et al. [23] mentioned that the estimated Sk would be inaccurate if hydraulic heads at piezometers do not represent the average local condition in the groundwater system. Most of the wells in this study correspond to a single confined aquifer based on the strainer depth information. Although there is a possibility In general, the uncertainty of Sk occurs for the following reasons.
(1) The Sk is biased if the horizontal displacement is not considered in the Sk calculation [81]. In our study area, horizontal fluid movement occurs over the groundwater system [23]. (2) Monitoring wells observe water from one aquifer in the saturated confined aquifer system [81]. Hoffmann et al. [23] mentioned that the estimated Sk would be inaccurate if hydraulic heads at piezometers do not represent the average local condition in the groundwater system. Most of the wells in this study correspond to a single confined aquifer based on the strainer depth information. Although there is a possibility that the temporal variation of the deeper aquifer may generate errors in estimating Sk, we assume the influence of the shallowest confined aquifer is dominant. (3) Error in measurements of InSAR displacement is due to atmospheric phase effects.
In this study, temporal filtering was used to mitigate the error of atmospheric contribution. Because PSI displacements were consistent with the F3 solution of GEONET, GSI, Japan, the error in measurements of InSAR displacement could be minor. (4) Incomplete removal of the long-term subsidence from the land deformation time series for case studies of high subsidence rate. To completely separate long-term trends of subsidence and hydraulic head from their time series, daily or weekly measurements of InSAR and head time series for several years are required [82].
Although these factors generate uncertainty for the estimated Sk, a part of uncertainty was suppressed by removing outliers of InSAR displacements and groundwater-level data using a Z-score.
To validate the Sk-Vs relationship, furthermore, we applied a statistical method called repeated holdout cross-validation. In this method, the original data are split into two data sets. One data set is used for testing (validation) of the model, and another is used as a training data set. Validation and training data sets contain four and nine data points, respectively. Exponential curve fitting was applied to the training data set, and the resulting model was used to evaluate the root mean square error (RMSE) of the testing data set as follows: where Z (x i ) and Z (y i ) are measured and estimated Sk values of the i sampling point, respectively, and N is the total number of observation points. In our case, N corresponds to the number of validation data sets (N = 3). The method is repeated 20 times (trials) to estimate the RMSE of each trial ( Figure 10). The low RMSE (averaged RMSE = 8.1 × 10 −3 ) indicates that the constructed equation (i.e., Equation (5)

Conclusions
This study proposes a new approach to map storage coefficients from surface displacements, groundwater-level data, and a high-resolution Vs model. Its main findings are summarized as follows.

•
The zone of low Vs found by the microtremor survey could have coincided with the Futagawa fault zone; • Sk of the confined aquifer ranges between ~0.03 and 2 × 10 −3 , with an average of 7.23 × 10 −3 , reflecting semi-confined and confined conditions; • An empirical relationship between the Sk and Vs was found, indicating that aquifer compressibility is linked to its stiffness and Vs; • The map of Sk estimated from the empirical relationship correlates with the hydrogeological setting and can be used to estimate the spatiotemporal variation of groundwater-level based on the geodetic data.
In conclusion, our approach can effectively estimate aquifer storage properties from S-wave velocities even where limited groundwater level data are available.  Thus, if we can validate the parameters used in our model, this approach could be used to map the Sk and temporal variation of groundwater levels from the surface deformation (geodetic data) at a lower cost.

Conclusions
This study proposes a new approach to map storage coefficients from surface displacements, groundwater-level data, and a high-resolution Vs model. Its main findings are summarized as follows.

•
The zone of low Vs found by the microtremor survey could have coincided with the Futagawa fault zone; • Sk of the confined aquifer ranges between~0.03 and 2 × 10 −3 , with an average of 7.23 × 10 −3 , reflecting semi-confined and confined conditions; • An empirical relationship between the Sk and Vs was found, indicating that aquifer compressibility is linked to its stiffness and Vs; • The map of Sk estimated from the empirical relationship correlates with the hydrogeological setting and can be used to estimate the spatiotemporal variation of groundwater-level based on the geodetic data.
In conclusion, our approach can effectively estimate aquifer storage properties from S-wave velocities even where limited groundwater level data are available. Funding: The first author was funded by a full scholarship from the Ministry of Higher Education of the Arab Republic of Egypt, through a grant-in-aid by the Egypt-Japan Education Partnership (EJEP) for developing human resources, which is regarded as an initial goal in the long-term development plan in Egypt. This research is partially supported by JSPS KAKENHI (JP20H01997).
Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Data sharing is not applicable.