2018 | 71/2 | 83–90 | 3 Figs. | 3 Tabs. | 1 App. | www.geologia-croatica.hr Journal of the Croatian Geological Survey and the Croatian Geological Society 1. INTRODUCTION In many small islands in the oceans and seas composed of hy- draulically permeable geologies, freshwater lenses are the only source of abundant freshwater supply for the residents and oth- ers, as well as for agriculture that is often an important industry on these islands. Such groundwater resources are known or at least considered to be susceptible to natural and anthropogenic threats such as long-term unprecedented sea-level rise, extended drought, and over-abstraction (e.g., TERRY & CHUI, 2012; ISHIDA et al., 2015; BARKEY & BAILEY, 2017). Understand- ing the hydraulic properties of aquifers is essential for appropri- ate management of the groundwater resources that are both vital and vulnerable. One conventional method for studying the hydraulic proper- ties of an aquifer is a pumping test. However, this is not always the best approach for a freshwater-lens aquifer for several rea- sons, including the following: (i) when a pumping test is con- ducted within an area where fresh groundwater is underlain by saltwater, as an environmental concern, the pumping may cause an upflow of the underlying saltwater that can contaminate the freshwater resource, and (ii) when the test is alternatively con- ducted at a site close to the coastline avoiding the fresh-ground- water distribution, technologically, the analysis of the pumping- test data requires some strategies to compensate for large tidal disturbances common in such locations (JHA et al., 2003; CHAT- TOPADHYAY et al., 2015; HOUBEN & POST, 2017). Heterogeneous hydraulic properties of an insular aquifer clarified by a tidal response method with simple decomposition techniques Katsushi Shirahata, Shuhei Yoshimoto, Takeo Tsuchihara and Satoshi Ishida NARO Institute for Rural Engineering, Tsukuba, Ibaraki 305-8609, Japan (shirahatak@affrc.go.jp) doi: 10.4154/gc.2018.06 Abstract Two simple frequency decomposition techniques were used as part of a tidal response method to derive the hydraulic diffusivities of a freshwater-lens aquifer. Digital high-pass filtering can separate the tidal components of diurnal and shorter periods from longer-period components. Discrete Fourier transform can be used to isolate a specific tidal component. These techniques are easy to practice using the built-in functions of spreadsheet software. The applied techniques were each optimized for the frequencies of known major tidal components. Isolation of the spe- cific tidal signals helps to reduce the errors of a basic tidal response method that uses in its cal- culations the amplitude attenuation and phase lag of a simple sinusoidal wave of groundwater fluctuations. Another advantage of the present tidal method is the utilization of two groundwater time series collected from near-shore and relatively inland sites affected by the same ocean tide. The method does not use surface-water observation, thus avoiding errors derived from gene- rally possible surface-water/groundwater boundary effects. The tidal response method with simple decomposition techniques was used to investigate the aquifer properties of an uplifted limestone island located in a subtropical region of Japan. A fresh- water lens is the principal water resource for this island and its sustainable development is de- sired. Significant hydraulic layering has not been reported in the limestone aquifer. Pairs of groundwater-level time-series data collected by simultaneous observations at near-shore and inland sites were analysed by the tidal response method. The results demonstrated heteroge- neous aquifer diffusivity on the island, typically with larger values in the southeastern coastal part than in the northwestern coastal part, which is consistent with the planar distribution of the entire freshwater lens and the position of its maximum thickness that are slightly biased toward the northwestern side. In studies of the hydraulic properties of aquifers bearing wa- ter affected by tidal oscillations, tidal methods have found world- wide application (e.g., CARR & VAN DER KAMP, 1969; KRIVIC, 1982; ALMEIDA & SILVA, 1983; KOIZUMI et al., 1998; TREFRY & JHONSTON, 1998; CORBETT et al., 2000; BANERJEE et al., 2008; ALCOLEA et al., 2009; CAROL et al., 2009; FADILI et al., 2012; MARTIN et al., 2012; AUSTIN et al., 2013; ROTZOLL et al., 2013; JO et al., 2014; PERRIQUET et al., 2014; YANG et al., 2015; NIETO LÓPEZ et al., 2016; GUO et al., 2017). The methods analyse the natural tidal oscillations observed in groundwater height or pressure propagated from adjacent sur- face waters, and do not require anthropogenic pumping. Analysis of the natural tidal response of groundwater fluctuations gives a general idea of the aquifer hydraulic parameters (PERRIQUET et al., 2014) and would be an appropriate and convenient approach for studying a freshwater-lens aquifer. SHIRAHATA et al. (2014) used a tidal response method with a simple decomposition technique, following a fundamental cal- culation procedure for discrete Fourier transform to extract or isolate four known major tidal components. They used a ground- water-level observation time series covering 369 days, a length optimized for the periods of the four major tides. SHIRAHATA et al. (2016) suggested the use of nonrecursive digital filters to separate a semidiurnal to diurnal tidal-period band from a longer- period weather band prior to the quantitative analysis of tidal components. SHIRAHATA et al. (2017) added appropriate time Article history: Manuscript received December 15, 2017 Revised manuscript accepted April 23, 2018 Available online June 21, 2018 Keywords: freshwater-lens aquifer, nonrecursive digital filtering, discrete Fourier transform, tidal response method, hydraulic heterogeneity G eo lo gi a C ro at ic a Geologia Croatica 71/284 series lengths for the isolation of major tidal components, the shortest of which was 29.5 days, based on examination of artifi- cially synthesized various time series. The purpose of the present work was (i) to confirm the va- lidity of nonrecursive digital filtering, simple discrete Fourier transform, and their integration with a tidal response method, suggested by the recent studies mentioned above, through appli- cation to observation time series of real coastal groundwater, and (ii) to estimate the hydraulic properties of the island’s aquifer where sustainable development of a freshwater lens has long been desired. 2. HYDROGEOLOGICAL OUTLINE OF THE STUDY ISLAND Tarama Island is one of the Ryukyu Islands that comprise part of the Japanese Islands. The oval-shaped island measures about 5.8 km × 4.3 km with an area of approximately 19.8 km2 (Fig. 1). Most of the land area is 5–15 m above sea level (a.s.l.). In the east- ern part of the island, an approximately N-S lineament is clearly detected in satellite and aerial photos, with a difference in ground elevation of roughly several metres observed in the field (higher on the east side). The lineament is thought to be a fault line with an approximately vertical fault plane (YAZAKI, 1977). Tarama Island is dominated by Quaternary limestone several tens of metres thick underlain by fine-grained sandstone (OOGA et al., 1974; OHZEKI et al., 2014). The limestone is hydraulically highly permeable with a reported hydraulic conductivity range of 1.05×10−2 to 2.91×10−2 m/s based on a pumping test conducted on the island (YAMADA et al., 2009). The exact site of the pump- ing test is not known. Significant hydraulic layering has not been reported in the limestone aquifer. The underlying fine-grained sandstone is regarded as practically impermeable, judging from the results of in situ permeability tests conducted in investiga- tions by the Okinawa General Bureau, Cabinet Office of Japan. Estimated hydraulic conductivities were in the order of 10–8 to 10–6 m/s. Contour lines of roughly estimated elevations of the top sur- face of the sandstone are, in principle, drawn concentrically with the highest part near the centre of the island and the lower part in the surrounding coastal area, so as to be consistent with test drill- ing records where the top of the sandstone was reached (Fig. 1A). The water table averages several decimetres a.s.l. in the central part of the island. The thickness of the limestone aquifer on this island is roughly estimated as 35 m near the centre of the island and more than 60 m in the coastal area, except within the eastern part of the fault. The above descriptions of the hydrogeological structure of the island are mostly based on unpublished results of investigations conducted by the local branch of the national ad- ministrative office in charge of this area. Owing to the extended permeable limestone, rainwater rea- dily percolates down into the ground and provides no surface water on Tarama Island, except where runoff water on paved sur- faces is artificially collected and stored in agricultural ponds. Naturally, groundwater is the principal freshwater resource of the island. A freshwater lens is the source of domestic-water supply for a population of about 1200 residents. On the other hand, the local branch of the national administrative office is investigating the groundwater resource to make a general plan for agricultural and rural development with an appropriate combination of sur- face-water and groundwater resource development. The freshwa- ter lens covers an area of roughly 10 km2 and is about 7 m thick at its maximum (SHIRAHATA, 2010; ISHIDA et al., 2011) (Fig. 1B). The horizontal positions of the whole freshwater lens and the maximum thickness are slightly offset to the northwest- ern side of the island. 3. TIME-SERIES DATA COLLECTION As described by SHIRAHATA et al. (2017) and referred to later in this paper, the present tidal response method requires a pair of sites for groundwater observations with different distances to the tidal surface water. The present study used time-series data pre- viously reported by SHIRAHATA et al. (2014). The groundwater- level data were collected from five observation sites, one of which is located only 0.02 km inland from the south-southwestern edge of the island (site A) and four others are located further inland at distances of 1.0 to 1.3 km to each nearest coastline of the east- southeastern, southern, western, and north-northwestern sides of the island (sites B, C, D, and E, respectively) (Fig. 1B). The four 1 m 2 m 3 m 4 m5 m 6 m 7 m 6 m 5 m 1 m 1 m 1000 m N site A B C D E B 1000 m N 40 m bsl 33.6 m 38.1 m 39.2 m 49.0 m 50.8 m bsl 50 m 60 m 40 m A Figure 1. Estimated contour map of the top surface of the impermeable sandstone basement beneath the limestone aquifer of Tarama Island (A) and the distribu- tion of the freshwater lens in the aquifer (B). (A) The map of the basement is made by tracing contour lines from a map created in investigations by the Okinawa General Bureau, Cabinet Office of Japan, a map in accordance with test drilling records at five sites (crosses in the map, with the elevations of the top surface of the basement below sea level). (B) The isopachs of the freshwater aquifer (measured EC less than 200 mS/m) are based on measurements taken on 1st October 2008, previously partly drawn by SHIRAHATA & NAGATA (2009). Solid circles are locations of the groundwater-level observation holes, sites A through E, where time-series data were collected. Dotted lines represent the shortest paths from the respective observation sites to the coastline, indicating the expected propagation paths of tidal waves observed in the groundwater fluctuations. G eologia C roatica Shirahata et al.: Heterogeneous hydraulic properties of an insular aquifer clarified by a tidal response method with simple decomposition ... 85 inland sites were each paired with the near-shore site as required for the tidal response method used and four pairs of sites were made (A/B, A/C, A/D, A/E). It was assumed both in the study by SHIRAHATA et al. (2014) and here, as an approximation, that the tidal oscillations of groundwater levels 0.02 km inland of the above four coastlines nearest to the inland observation sites are identical to the oscillations at the near-shore observation site. Hourly sampled time-series data were collected from June 2008 to July 2009 at groundwater-observation holes using sub- mersible water-level data loggers with a built-in sensor. SHIRA- HATA et al. (2014) showed in their figs. 5 and 6 the time-series charts of the observed groundwater levels. The time series in- cluded 1- to 3-h data gaps or missing data on rare occasions. Data gaps with a length of 1 h or one data point in the time series were filled with the average of the previous and next data values. Digi- tal filtering and Fourier transform were applied to consistently hourly time-series data after filling the gaps. 4. FREQUENCY DECOMPOSITION TECHNIQUES AND TIDAL RESPONSE METHOD The decomposition techniques and tidal response method em- ployed in the present study are explained below. The simple dis- crete Fourier-transform technique was used in two different steps for different aims. First it was used to make frequency spectra of the groundwater fluctuations. Then the technique was utilized to provide amplitude and initial phase of one desired tidal compo- nent, which were subsequently used in formulae to calculate aqui- fer diffusivity. The two procedures of the Fourier transform are described separately, with an explanation of digital filtering in between. The basics of the techniques for discrete Fourier trans- form and digital filtering are also mentioned. 4.1. Fourier transform to produce spectra with the amplitudes of major tides Major tidal components, called “tidal constituents,” are identified by the frequency or period of the sinusoidal tide-generating po- tentials. SHIRAHATA et al. (2014) utilized a Fourier technique for isolation of the top four major tidal constituents, M2, K1, S2, and O1, from 8856-h-long hourly time-series data. The isolation means the determination of the amplitude and initial phase of the sinusoidal tidal oscillations. The technique simply followed (part of) the fundamental formulae for discrete Fourier transform and was completed using spreadsheet software with built-in functions for trigonometric calculations, averages, square roots, and other basic operations (APPENDIX). The appropriate choice of length of the transformed time series is essential for accurate isolation of the desired tidal components of specific frequencies. There are several recommendable combinations of the transformed time- series lengths and isolated major tidal constituents for the simple Fourier-transform technique (SHIRAHATA et al., 2017). In the present study, year-long time-series data collected at sites A, B, C, and E were firstly subjected to the Fourier transform to describe the frequency compositions of the groundwater fluc- tuations. The data of site D were not employed because they con- tained unfilled 2- and 3-h-long data gaps. Simple Fourier transform of 8856-h-long data from 2008-7-1 0:00 to 2009-7-4 23:00 was performed to produce a spectrum that exhibits amplitudes at the approximate or exact frequencies of major tides O1, P1, K1, M2, S2, and K2. In the obtained spectrum, tidal frequencies will be indicated on the horizontal axis by inte- gers n used in the transform calculations, the integers for divi ding the time-series length (8856) to give quotients close to the spe- cific periods of the tidal components (SHIRAHATA et al., 2017). For example, the amplitude for the M2 frequency will be shown at n = 713, because the period 12.420601 h is very closely ap- proximated by the quotient 8856 divided by 713. In the same man- ner, the amplitudes for the frequencies of major tidal constituents O1, P1, K1, S2, and K2 will be shown at n = 343, 368, 370, 738, and 740, respectively. To supplement the measurement of the amplitude of another major tide, N2 constituent (period: 12.658348 h), 9924-h-long time-series data from 2008-6-12 0:00 to 2009-7-30 11:00 were also transformed. In the spectrum, the amplitudes for the N2 and M2 frequencies will be shown at n = 784 and 799, respectively. 4.2. Nonrecursive digital filtering to remove long-period components Digital filters are used to separate out fluctuation components with a desired range of frequencies from discrete-time time- series data. The digital filter employed in the present study is a nonrecursive type. In other words, one data point of the filtered output discrete time series is calculated by the linear combination of terms of a number sequence that is a finite-length portion of the input time series. The multiplication coefficients of the linear combination comprise the number sequence that represents the applied digital filter. Nonrecursive digital filtering of electromag- netic time-series data with a constant sampling interval is easily achieved using standard spreadsheet software with a built-in function for calculating the sum of the products of the corre- sponding terms of two number sequences. SHIRAHATA et al. (2016) provided a further introduction to nonrecursive digital fil- tering and concrete examples of conventional and new digital fil- ters for tidally fluctuated time-series data. The high-pass filter used in the present study is the reverse of the “LP241H079122kM3” low-pass filter made by SHIRA- HATA et al. (2016). The employed high-pass filter removes com- ponents with periods longer than 2 days, but accurately preserves in the output hourly time series the major diurnal and semidiurnal tidal constituents (Q1, O1, P1, K1, N2, M2, S2, and K2) contained in the input hourly time series. This filter has a length of 241 terms or hours and outputs time series shorter than the finite in- put time series by 10 days. Groundwater-level time-series data of the five observation sites underwent the filtering. Because the subsequent Fourier transform requires data covering 29.5 days, segments of 39.5-day-long time series were selected to be filtered. Two segments were extracted from the time series of each obser- vation site. From the data acquired at five sites from late 2008 to early 2009, the time series from 2008-9-6 0:00 to 2008-10-15 11:00 and the time series from 2009-1-11 0:00 to 2009-2-19 11:00 were extracted. These segments of time series were selected be- cause of the absence of data gaps in the raw data and the relatively large disturbance of water levels caused by the weather, namely atmospheric pressure changes, in order to clearly demonstrate the effect of digital filtering. 4.3. Fourier transform to isolate a tidal constituent for tidal response method The simple Fourier-transform technique was again used to isolate a tidal component for the tidal response method from short time- series data of the five observation sites. In the present step to iso- late a tidal constituent, part of the calculations for the Fourier transform, to determine one component, was utilized. Transform of time series with a length of 708 h, the shortest of the time-se- ries lengths that SHIRAHATA et al. (2017) recommended, was G eo lo gi a C ro at ic a Geologia Croatica 71/286 adopted and M2 constituent was isolated. Ten time-series seg- ments of the five observation sites, each 708 h long after shorte- ned by 240 h by high-pass filtering, were transformed to give the amplitudes and initial phases of the oscillations of M2 constituent contained in the groundwater fluctuations. Table 1 summarizes the time-series data subjected to Fourier transform in the two steps described above (see 4.1. for the first step). Each transformed time series is specified by a combination of the observation site and the start and end of the transformed period. Distances of the observation sites to each nearest coast- line are also shown. The distances were used in the tidal response method explained in the following section. 4.4. Basic tidal response method to derive aquifer hydraulic diffusivity Two basic formulae rearranged from FERRIS (1951) were used in this study (SHIRAHATA et al., 2017): D = (π/t)·(XB − XA)2/(ln(hA/hB))2 and D = (π/t)·(XB − XA)2/(Δω)2 , where D is the desired hydraulic diffusivity of the aquifer (equal to the ratio of transmissivity T to storage coefficient S, or to the ratio of conductivity K to specific storage Ss) [m2/s], t is the pe- riod of the tidal constituent under consideration [s], XA and XB are the distances of the near-shore site and inland site nearest to each boundary with the same tidal surface water (XA < XB) [m], hA and hB are the constituent amplitudes at the near-shore site and inland site (hA > hB) [m], and Δω is the phase lag of the inland site to the near-shore site (Δω > 0) [rad]. In the present case, the period t in the formula was set to 44714.164, the period of M2 constituent in units of seconds. Tidal- amplitude ratio (hA/hB) and phase shift (Δω) for M2 constituent were calculated from the outputs of the simple Fourier transform of the 708-h-long time series of the paired two sites. For each pair of near-shore and inland sites, four diffusivity values were calcu- lated from the two different segments of time series using the two formulae above. The four values were averaged to derive a de- finitive aquifer diffusivity for the pair of observation sites. 5. RESULTS 5.1. Major tidal components in water fluctuations Figure 2 shows amplitude-frequency spectra obtained from the time series of groundwater levels at site A using the simple Fou- rier-transform calculations. In the spectrum of the 8856-h-long time series (Fig. 2, top panel), it is shown that the fluctuations contain large signals of semidiurnal tidal constituents, especially M2 with an amplitude of 0.466 m (at n = 713). The second-largest component with an amplitude of 0.196 m at n = 738 corresponds to S2 constituent. The amplitude for the frequency of K2 consti- tu ent (n = 740) is 0.068 m. A rise in the spectrum around n = 700 with a peak amplitude of 0.063 m corresponds to N2 tide, but the amplitude for N2 frequency is not directly shown in this spectrum produced from the 8856-h time series. Instead, using the Fourier transform of the 9924-h time series (Fig. 2, bottom panel), the amplitude for the N2 frequency (n = 784) was determined as 0.087 m. Diurnal tidal components are relatively insignificant with amplitudes of 0.174, 0.055, and 0.191 m for the frequencies of O1, P1, and K1 (n = 343, 368, and 370 in the spectrum of the 8856-h time series), respectively. Table 2 lists the amplitudes for the frequencies of the major diurnal and semidiurnal tidal constituents estimated by the Fou- rier transform of time series of specified lengths. As an example of the inland sites, groundwater fluctuations at site C have signals for the frequencies of O1, P1, K1, N2, M2, S2, and K2 constituents with amplitudes of 0.026, 0.006, 0.027, 0.005, 0.029, 0.011, and 0.005 m, respectively. The diurnal tidal oscillations were attenu- ated compared to the near-shore site (A) with amplitude ratios between 11% and 16%, whereas for semidiurnal constituents the amplitude ratios were between 5% and 7%. Among the inland sites (B, C, and E) with similar distances to the coast (see Table 1), the amplitudes of the constituents at site B, the eastern site, are four to eight times those at site E, the north- western site. This indicates heterogeneity in the hydraulic pro- perties of the aquifer of the island, which is quantitatively eluci- dated later. 5.2. Effects of digital filtering and tidal-component isolation on time series Figure 3 exemplifies the effects of digital filtering and isolation of a major tidal component. It shows the sets of time-series line charts of the original groundwater-level observation data, low- pass-filtered and high-pass-filtered data, and the chart of the iso- lated M2 constituent, the latter reproduced from the outputs of the 708-h Fourier transform, for sites A and D. The digital filter- ing effectively separated semidiurnal to diurnal tidal signals (“High-passed”) from longer-period changes (“Low-passed”). Specifically, as an example of the longer-period changes, the low- passed data of site A (Fig. 3, top panel) reveals two rises, each a few days long with peaks on the 12th and 28th September that were unclear from the original data chart. The rises coincided with the approach of typhoons No. 13 and No. 15 of this year, re- Table 1. Observation sites and periods of Fourier-transformed time-series data. Site Closer side / Distance to coast Transformed period to produce spectrum Transformed period to isolate M2 Start / End (yyyy-m-d hh:) Start / End (yyyy-m-d hh:) Start / End (yyyy-m-d hh:) Start / End (yyyy-m-d hh:) A SSW 2008-7-1 00: 2008-6-12 00: 2008-9-11 00: 2009-1-16 00: 0.02 km 2009-7-4 23: 2009-7-30 11: 2008-10-10 11: 2009-2-14 11: B ESE 2008-7-1 00: 2008-6-12 00: 2008-9-11 00: 2009-1-16 00: 1.07 km 2009-7-4 23: 2009-7-30 11: 2008-10-10 11: 2009-2-14 11: C S 2008-7-1 00: 2008-6-12 00: 2008-9-11 00: 2009-1-16 00: 1.05 km 2009-7-4 23: 2009-7-30 11: 2008-10-10 11: 2009-2-14 11: D W 2008-9-11 00: 2009-1-16 00: 1.24 km 2008-10-10 11: 2009-2-14 11: E NNW 2008-7-1 00: 2008-6-12 00: 2008-9-11 00: 2009-1-16 00: 1.22 km 2009-7-4 23: 2009-7-30 11: 2008-10-10 11: 2009-2-14 11: G eologia C roatica Shirahata et al.: Heterogeneous hydraulic properties of an insular aquifer clarified by a tidal response method with simple decomposition ... 87 spectively. The intensive low-pressure system would have af- fected the water table of the permeable island through the inter- mediary of changes in height of the surrounding surface water (VACHER, 1978a). For all five sites A through E, including not shown, line charts of the original data, high-passed data, and the isolated con- stituent generally match the timing of periodic peaks and troughs, demonstrating that the digital filtering and Fourier transform of 708-h-long time series successfully isolated the desired M2 con- stituent. 5.3. Aquifer hydraulic parameters Table 3 lists aquifer diffusivities calculated by the tidal response method using amplitudes and initial phases of the M2 constituent isolated from the 708-h-long time series. The aquifer diffusivities are close to the values previously derived from the 8856-h time 0.0 0.1 0.2 0.3 0.4 0.5 1 101 201 301 401 501 601 701 801 Am pl itu de (m ) n (frequency = n/8856 cph) Site A, 8856-h time series n = 738: S2 0.196 m n = 713: M2 0.466 m n = 700 0.063 m n = 370: K1 0.191 m n = 368: P1 0.055 m n = 343: O1 0.174 m n = 740: K2 0.068 m 0.0 0.1 0.2 0.3 0.4 0.5 1 101 201 301 401 501 601 701 801 901 Am pl itu de (m ) n (frequency = n/9924 cph) Site A, 9924-h time series n = 799: M2 0.466 m n = 784: N2 0.087 m Figure 2. Frequency spectra of groundwater-level fluctuations at site A produced using the simple Fourier-transform technique. The two panels show spectra ob- tained from 8856- and 9924-h-long time series. Frequencies in cycles per hour (cph) are given by the quotients of the integers n on the horizontal axis divided by the transformed time-series length (8856 or 9924). Table 2. Amplitudes of major tidal constituents determined using Fourier transform of 8856- and 9924-h-long time series. Site Time-series length Amplitude of tidal constituent (m) O1 P1 K1 N2 M2 S2 K2 A 8856 h 9924 h 0.174 0.055 0.191 0.087 0.466 0.466 0.196 0.068 B 8856 h 9924 h 0.036 0.009 0.038 0.008 0.045 0.044 0.017 0.007 C 8856 h 9924 h 0.026 0.006 0.027 0.005 0.029 0.030 0.011 0.005 E 8856 h 9924 h 0.009 0.002 0.008 0.001 0.006 0.006 0.002 0.001 Figure 3. Examples of sets of time-series line charts of the original groundwater-level observation data, low-pass-filtered data, high-pass-filtered data, and the chart of isolated M2 constituent. Top and bottom panels display sets of data of sites A and D, respectively. -2 -1 0 1 2 3 -3 -2 -1 0 1 2 9/6 00:00 9/11 00:00 9/16 00:00 9/21 00:00 9/26 00:00 10/1 00:00 10/6 00:00 10/11 00:00 W at er le ve l ( m ) W at er le ve l ( m a sl ) Site A year 2008 Original data (left axis) Low-passed (left axis) High-passed (right axis) Isolated M2 (right axis) -0.4 -0.3 -0.2 -0.1 0.0 0.1 0.2 0.3 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 9/6 00:00 9/11 00:00 9/16 00:00 9/21 00:00 9/26 00:00 10/1 00:00 10/6 00:00 10/11 00:00 W at er le ve l ( m ) W at er le ve l ( m a sl ) Site D year 2008 Original data (left axis) Low-passed (left axis) High-passed (right axis) Isolated M2 (right axis) G eo lo gi a C ro at ic a Geologia Croatica 71/288 ssary to at least increase the number of case studies using the present method, including application to confined aquifers. The definitive diffusivity of the east-southeastern part of the aquifer (represented by the result for the pair of sites A/B) is two to three times that of the northwestern part of the aquifer (A/E pair). This difference is beyond the variations of the tentative dif- fusivities. The same magnitude relationship naturally holds for the estimated conductivities. This relationship is clearly consis- tent with the areal distribution of the freshwater lens offset to the northwestern side of the island (Fig. 1B), in a similar way to the case of Bermuda Island presented by VACHER (1978b). 7. CONCLUDING REMARKS This paper reported examples of applying simple techniques for processing tidally fluctuated groundwater time-series data, one to separate high- and low-frequency signals and the other to iso- late known major tidal components. The techniques were further used as part of a simple tidal response method that has recently been improved for application to relatively short observation data. The improved method combined with the two decomposition techniques was first applied to real groundwater observation data in the present work. The method provided convincing hydraulic parameters that compare well with the result of a pumping test and are consistent with the distribution of the freshwater lens on the study island. The present results showed the heterogeneity of the hydraulic properties of the aquifer on the island, which is now being used as the basis for the investigation of sustainable ground- water development on the island. The techniques and method are easy to practice using the built-in functions of standard spreadsheet software and are highly applicable to digitally recorded observation time-series data. They should be useful for investigating aquifer hydraulic proper- ties in insular and coastal areas where groundwater is often the principal, or sometimes the only, water resource. ACKNOWLEDGEMENT The authors are grateful to the Land Improvement General Office and the successive Officers for Policy Planning in the Land Im- provement Division of the Okinawa General Bureau, Cabinet Of- fice of Japan for their support and collaboration in field and desk research. Thanks are also due to local government officials and residents of Tarama Island for their understanding and coopera- tion in the field study. This work was partly supported by a re- search project funded by the Ministry of Agriculture, Forestry and Fisheries of Japan (Development of mitigation and adaptation technologies to climate change in the sectors of agriculture, fore- stry, and fisheries, 91150), and also by JSPS KAKENHI Grant No. 17K08011. series (SHIRAHATA et al., 2014). This confirms the validity of the present method based on the analysis of 708-h time series. The definitive aquifer diffusivities for the four pairs of sites A/B, A/C, A/D, and A/E, determined as the average of the four tentative values, range from 6.1 to 15.5 m2/s. Assuming the sto- rage coefficient of the limestone aquifer of this island as roughly 0.1, referring to the value used in public works on a neighbouring island to construct underground dams in the contemporaneous limestone formation, and assuming the aquifer thickness as roughly 55 m, based on the contour map of the top surface of the underlying impermeable stratum drawn in investigations by the national administrative office, hydraulic conductivities are calcu- lated as 1.1×10−2 to 2.8×10−2 m/s. 6. DISCUSSION The estimated hydraulic conductivities using the present method fall in the same order of magnitude as the range of values based on a pumping test referred to earlier (1.05×10−2 to 2.91×10−2 m/s; YAMADA et al., 2009). Unfortunately, the result of the pumping test was only briefly mentioned, and the exact test location is un- published. It is impossible to make a detailed comparison be- tween the results from the present tidal response method and the previous pumping test. Although there seems to be nothing amiss with observation data or the isolation of a tidal constituent (5.2.), the tentative hy- draulic diffusivities for each site pair show a consistent bias; dif- fusivities derived from the amplitude ratio are smaller than those from the phase lag. The same tendency is found in the results of SHIRAHATA et al. (2014), who used a tidal response method based on 8856-h time-series Fourier decomposition. Thus, the detected bias does not imply a defect specific to the present method based on 708-h time-series data. Among previous stu- dies using tidal response methods, the same bias between attenua- tion-based and lag-based aquifer parameters, with much more significance than that of the present case, is reported by ER- SKINE (1991), who attributed the variations in the calculated aq- uifer parameters to the damping effect of a phreatic surface on fluctuation amplitudes. The same bias, to an extent slightly lesser than that of ERSKINE (1991), is also detected by TREFRY & BEKELE (2004) between the amplitude attenuation and time de- lay of tidal signals in their study aquifer. They deduced that the less than ideal relationship between the attenuation and delay was due to horizontal layering in the aquifer properties. In the present study, there is no definitive evidence as to the reason for the dif- ference in the tentative hydraulic diffusivities for each site pair, but the effect of a phreatic surface or unreported hydraulic layer- ing is a possible cause. To elucidate the reason, it would be nece- Table 3. Aquifer hydraulic diffusivities estimated by the simple tidal response method using M2 constituent isolated by Fourier transform of 708-h-long time series. Site pair Length / Location Start of time series (Near-shore / Inland) of target aquifer (yyyy-m-d hh:) M2 Amp M2 IP M2 Amp M2 IP derived from Amp derived from IP 1.05 km 2008-9-11 00: 0.440 6.048 0.043 3.885 14.4 16.6 15.5 / ESE 2009-1-16 00: 0.461 2.281 0.044 0.131 14.1 16.8 1.4 ( 9.0%) 1.03 km 2008-9-11 00: 0.440 6.048 0.028 3.564 9.9 12.1 10.8 / S 2009-1-16 00: 0.461 2.281 0.028 6.055 9.4 11.8 1.3 (12.4%) 1.22 km 2008-9-11 00: 0.440 6.048 0.023 3.497 12.0 16.1 14.1 / W 2009-1-16 00: 0.461 2.281 0.024 6.031 11.9 16.3 2.4 (17.4%) 1.20 km 2008-9-11 00: 0.440 6.048 0.005 2.208 5.1 6.9 6.1 / NW 2009-1-16 00: 0.461 2.281 0.006 4.764 5.3 7.0 1.0 (16.8%) Amp: amplitude (m), IP: initial phase (rad), D: hydraulic diffusivity (m2/s), Avg: average, SD: standard deviation (Avg / 1 SD of left four) Definitive DNear-shore site Inland site Tentative D A / B A / C A / D A / E G eologia C roatica Shirahata et al.: Heterogeneous hydraulic properties of an insular aquifer clarified by a tidal response method with simple decomposition ... 89 KRIVIC, P. (1982): Transmission des ondes de marée à travers l’aquifère côtier de Kras [Transmission of tidal waves across the coastal aquifer of Kras – in French].– Ge- ologija, 25/2, 309–325. MARTIN, J.B., GULLEY, J. & SPELLMAN, P. (2012): Tidal pumping of water between Bahamian blue holes, aquifers, and the ocean.– Journal of Hydrology, 416–417, 28–38. doi: 10.1016/j.jhydrol.2011.11.033 NIETO LÓPEZ, J.M., ANDREO NAVARRO, B. & MUDARRA MARTÍNEZ, M. (2016): Hydrogeological parameters assessment by tidal influence analysis in the coastal aquifers of Bajo Guadalhorce (Malaga province, southern Spain) [in Spa- nish, with an English abstract].– Geogaceta, 59, 39–42. OHZEKI, M., IMAI, R., TAKAYANAGI, H. & IRYU, Y. (2014): Stratigraphy and geo- logic age of the Ryukyu Group on Tarama-jima, Ryukyu Islands, Japan [in Japa- nese].– Abstracts of the 121st Annual Meeting of the Geological Society of Japan, R5-P-19. OOGA, H., FURUKAWA, H. OGURA, I. & NISHIDA, T. (1974): Groundwater of Tarama Island, Okinawa Prefecture [in Japanese].– Abstracts of the 81st Annual Meeting of the Geological Society of Japan, 368. PERRIQUET, M., LEONARDI, V., HENRY, T. & JOURDE, H. (2014): Saltwater wedge variation in a non-anthropogenic coastal karst aquifer influenced by a strong tidal range (Burren, Ireland).– Journal of Hydrology, 519B, 2350–2365. doi: 10.1016/j. jhydrol.2014.10.006 ROTZOLL, K., GINGERICH, S.B., JENSON, J.W. & EL-KADI, A.I. (2013): Estima- ting hydraulic properties from tidal attenuation in the Northern Guam Lens Aquifer, territory of Guam, USA.– Hydrogeology Journal, 21/3, 643–654. doi: 10.1007/ s10040-012-0949-9 SHIRAHATA, K. (2010): Analysis of water balance in formation of freshwater lens through electric conductivity measurement [in Japanese].– Journal of the Japanese Society of Irrigation, Drainage and Rural Engineering, 78/6, 514–515. SHIRAHATA, K. & NAGATA, J. (2009): A study of freshwater lens in Taramajima Is- land for a large-scale water resource development in the future [in Japanese].– Geotechnical Engineering Magazine, 57/9, 42. SHIRAHATA, K., ISHIDA, S., YOSHIMOTO, S. & TSUCHIHARA, T. (2014): New simple method for estimating hydraulic properties of a freshwater-lens aquifer by analysis of tidal groundwater fluctuations [in Japanese, with an English summa- ry].– Technical Report of the National Institute for Rural Engineering, 215, 141–154. SHIRAHATA, K., YOSHIMOTO, S., TSUCHIHARA, T. & ISHIDA, S. (2016): Digital filters to eliminate or separate tidal components in groundwater observation time- series data.– Japan Agricultural Research Quarterly: JARQ, 50/3, 241–252. doi: 10.6090/jarq.50.241 SHIRAHATA, K., YOSHIMOTO, S., TSUCHIHARA, T. & ISHIDA, S. (2017): Im- provements in a simple harmonic analysis of groundwater time series based on error analysis on simulated data of specified lengths.– Paddy and Water Environ- ment, 15/1, 19–36. doi: 10.1007/s10333-016-0525-3 TERRY, J.P. & CHUI, T.F.M. (2012): Evaluating the fate of freshwater lenses on atoll islands after eustatic sea-level rise and cyclone-driven inundation: A modelling approach.– Global and Planetary Change, 88–89, 76–84. doi: 10.1016/j.glopla- cha.2012.03.008 TREFRY, M.G. & BEKELE, E. (2004): Structural characterization of an island aquifer via tidal methods.– Water Resources Research, 40/1, W01505. doi: 10.1029/2003WR002003 TREFRY, M.G. & JHONSTON, C.D. (1998): Pumping test analysis for a tidally forced aquifer.– Ground Water, 36/3, 427–433. doi: 10.1111/j.1745-6584.1998.tb02813.x VACHER, H.L. (1978a): Hydrology of small oceanic islands-influence of atmospheric pressure on the water table.– Ground Water, 16/6, 417–423. doi: 10.1111/j.1745-6584.1978.tb03256.x VACHER, H.L. (1978b): Hydrogeology of Bermuda–Significance of an across-the-island variation in permeability.– Journal of Hydrology, 39/3–4, 207–226. doi: 10.1016/0022-1694(78)90001-X YAMADA, S., YONAHARA, N. & SOBUE, H. (2009): The Quaternary coral reef com- plex deposits (Ryukyu Group) and hydrogeologic features on Tarama-jima, Ok- inawa Prefecture, Japan [in Japanese].– Abstracts of the 116th Annual Meeting of the Geological Society of Japan, 83. YANG, H., SHIMADA, J., MATSUDA, H., KAGABU, M. & DONG, L. (2015): Evalu- ation of a freshwater lens configuration using a time series analysis of a ground- water level and an electric conductivity in Minami-daito Island, Okinawa Prefec- ture, Japan [in Japanese, with an English abstract].– Journal of Groundwater Hydrology, 57/2, 187–205. doi: 10.5917/jagh.57.187 YAZAKI, K. (1977): Geology of the Taramashima District [in Japanese, with an English abstract].– Geological Survey of Japan, Kawasaki, 28p. REFERENCES ALCOLEA, A., RENARD, P., MARIETHOZ, G. & BERTONE, F. (2009): Reducing the impact of a desalination plant using stochastic modelling and optimization techniques.– Journal of Hydrology, 365/3–4, 275–288. doi: 10.1016/j.jhy- drol.2008.11.034 ALMEIDA, C. & SILVA, M.L (1983): Novas observações sobre o efeito de maré em aquíferos costeiros do Algarve [New observations of the tidal effect in coastal aq- uifers of Algarve – in Portuguese, with an English abstract].– Boletim da Sociedade Geológica Portugal, 24, 289–293. AUSTIN, M.J., MASSELINK, G., MCCALL, R.T. & POATE, T.G. (2013): Groundwa- ter dynamics in coastal gravel barriers backed by freshwater lagoons and the po- tential for saline intrusion: Two cases from the UK.– Journal of Marine Systems, 123–124, 19–32. doi: 10.1016/j.jmarsys.2013.04.004 BANERJEE, P., SARWADE, D. & SINGH, V.S. (2008): Characterization of an island aquifer from tidal response.– Environmental Geology, 55/4, 901–906. doi: 10.1007/s00254-007-1041-y BARKEY, B.L. & BAILEY, R.T. (2017): Estimating the impact of drought on ground- water resources of the Marshall Islands.– Water, 9/1, 41. doi: 10.3390/w9010041 CAROL, E.S., KRUSE, E.E., POUSA, J.L. & ROIG, A.R. (2009): Determination of he- terogeneities in the hydraulic properties of a phreatic aquifer from tidal level fluc- tuations: a case in Argentina.– Hydrogeology Journal, 17/7, 1727–1732. doi: 10.1007/s10040-009-0478-3 CARR, P.A. & VAN DER KAMP, G.S. (1969): Determining aquifer characteristics by the tidal method.– Water Resources Research, 5/5, 1023–1031. doi: 10.1029/ WR005i005p01023 CHATTOPADHYAY, P.B., VEDANTI, N. & SINGH, V.S. (2015): A conceptual nume- rical model to simulate aquifer parameters.– Water Resources Management, 29/3, 771–784. doi: 10.1007/s11269-014-0841-6 CORBETT, D.R., DILLON, K. & BURNETT, W. (2000): Tracing groundwater flow on a barrier island in the north-east Gulf of Mexico.– Estuarine, Coastal and Shelf Science, 51/2, 227–242. doi: 10.1006/ecss.2000.0606 ERSKINE, A.D. (1991): The effect of tidal fluctuation on a coastal aquifer in the UK.– Ground Water, 29/4, 556–562. doi: 10.1111/j.1745-6584.1991.tb00547.x FADILI, A., MEHDI, K., RISS, J., MALAURENT, P.H., BOUTAYEB, K. & GUESSIR, H. (2012): Oceanic tidal influence on the piezometric level variation of the coast- al karst aquifer of Oualidia (Morocco) [in French, with an English abstract].– Af- rica Geoscience Review, 19/3, 135–150. FERRIS, J.G. (1951): Cyclic fluctuations of water level as a basis for determining aqui- fer transmissibility.– International Association of Scientific Hydrology, Publication 33, 148–155. GUO, M., WAN, J., JIANG, F. & HUANG, K. (2017): Estimating unconfined aquifer parameters based on groundwater tidal effect [in Chinese, with an English abstract] .– Earth Science, 42/1, 155–160. doi: 10.3799/dqkx.2017.012 HOUBEN, G. & POST, V.E.A. (2017): The first field-based descriptions of pumping- induced saltwater intrusion and upconing.– Hydrogeology Journal, 25/1, 243–247. doi: 10.1007/s10040-016-1476-x ISHIDA, S., TSUCHIHARA, T., YOSHIMOTO, S., MINAKAWA, H., MASUMOTO, T. & IMAIZUMI, M. (2011): Estimate of the amount of freshwater lens of Tarama Island, Japan [in Japanese, with an English abstract].– Irrigation, Drainage and Rural Engineering Journal, 273, 7–18. ISHIDA, S., YOSHIMOTO, S., KODA, K., KOBAYASHI, T., SHIRAHATA, K. & TSUCHIHARA, T. (2015): Salt water intrusion into groundwater and problem on Vava’u Island and Lifuka Island, Kingdom of Tonga [in Japanese, with an English abstract].– Technical Report of the National Institute for Rural Engineering, 217, 1–12. JHA, M.K., KAMII, Y. & CHIKAMORI, K. (2003): On the estimation of phreatic aq- uifer parameters by the tidal response technique.– Water Resources Management, 17/1, 69–88. doi: 10.1023/A:1023018107685 JO, S.-B., JEON, B.-C., PARK, E.-G., CHOI, K.-J., SONG, S.-H. & KIM, G.-P. (2014): Estimation of hydraulic characteristics and prediction of groundwater level in the eastern coastal aquifer of Jeju Island [in Korean, with an English abstract].– Jour- nal of Environmental Science International, 23/4, 661–672. doi: 10.5322/JE- SI.2014.4.661 KOIZUMI, N., KITAGAWA, Y., KAZAHAYA, K. & TAKAHASHI, M. (1998): Vol- canic gas concentration and aquifer permeability estimated from tidal fluctuations in groundwater level: Case of Koshimizu Well in Izu-Oshima, Japan.– Geophysi- cal Research Letters, 25/12, 2237–2240. doi: 10.1029/98GL01409 G eo lo gi a C ro at ic a Geologia Croatica 71/290 APPENDIX Example of spreadsheet for isolation of four major tidal compo- nents from 8856-h-long time series, using part of the fundamen- tal formulae for discrete Fourier transform (SHIRAHATA et al., 2014). If comparable calculations are made for consecutive natu- ral numbers, instead of only four numbers (343, 370, 713, 738) used in this spreadsheet, a Fourier spectrum as shown in Fig. 2 (top panel) can be made.