Estimations of tidal characteristics and aquifer parameters via tide-induced head changes in coastal observation wells

The groundwater fluctuations due to tidal variations at an observation well in a coastal aquifer can be used to determine the tidal characteristics and aquifer parameters without conducting an aquifer test. In this study, a method, comprised of Jeng et al.’s solution (2005) and simulated annealing (SA) algorithm, is developed to determine the coastal aquifer parameters (hydraulic diffusivity, beach slope, and aquifer thickness) as well as the tidal characteristics (bichromatic-tide amplitudes, bichromatic-tide wave frequencies, and tidal phase lag) from the analysis of the tideinduced well-water-level (WWL) data. The synthetic WWL data generated from Jeng et al.’s solution (2005) with assumed parameter values and field data obtained from Barrenjoey beach, Australia, are analyzed. The estimated parameter values obtained from analyzing synthetic WWL data by the present method show good agreements with the previously assumed parameter values. The parameter estimation procedure may however fail in the case of a large shallow water parameter which in fact violates the constraint on the use of Jeng et al.’s solution (2005). In the analysis of field WWL data, the results indicate that the aquifer parameters estimated from the present method with single or multiple well data are significantly different from those given in Nielsen (1990). Inspecting the observed WWL data and the WWL data predicted from Jeng et al.’s solution (2005) reveals that the present method may provide better estimations for the aquifer parameters than those given in Nielsen (1990). Correspondence to: H.-D. Yeh (hdyeh@mail.nctu.edu.tw)


Introduction
In coastal areas, groundwater levels of an aquifer fluctuate with tidal variations.The coastal aquifer parameters can be estimated from analyzing the well-water-level (WWL) data at an observation well without conducting conventional aquifer tests.Numerous investigations have been devoted to the study of hydraulics of tide-induced groundwater variation in coastal aquifers.Nielsen (1990) used a perturbation technique to obtain an analytical solution up to the secondorder of amplitude parameter (α = A/D) for tidal dynamics in sloping sandy beaches.In his model, the shoreline boundary condition at the interface of the beach and the ocean was allowed to vary with the tide height; it was however invalid in the situations of flat beach and/or large tidal range that may cause a seepage point deviated from the shoreline.Li et al. (2000) overcame this conflict by introducing a moving shoreline condition and considered the tides to be bichromatic, namely, tides can be represented by two different wave frequencies.Both models however treated the beach slope as a part of the perturbation parameter, which restricts the applicability of the models to the case of aquifers with large beach slopes.Moreover, both models relied on the Boussinesq equation that was solved only by the zero-order approximation in shallow water expansion.On the other hand, Teo et al. (2003) developed a higher-order analytical solution based on the shallow water expansions for the water table fluctuations induced by the monochromatic tide in a sloping coastal aquifer.Jeng et al. (2005) further considered the effect of bichromatic tide in the development of solution for WWL fluctuations in a sloping coastal aquifer.
Identifying the aquifer parameters and tidal characteristics from the analysis of WWL data can be cast as a minimization problem.Simulated annealing (SA), first proposed by Metropolis et al. (1953), is a technique constructed on the statistical mechanics for solving the optimization problems.The concept of the algorithm is based on simulating the recrystallization of a material in the process of annealing.SA has the capacity of dealing with the complicated problems involving multi-degrees of freedom and several local optima.In hydrological engineering, SA has been widely applied to solve various types of optimization problems (e.g., Marryott et al., 1993;Pardo-Iguzquiza, 1998;Huang and Yeh, 2007;Yeh and Chen, 2007;Chen and Yeh, 2009).
In this study, we propose a method to estimate the coastal aquifer parameters as well as tidal characteristics, including hydraulic diffusivity, beach slope, aquifer thickness, bichromatic-tide amplitudes, bichromatic-tide wave frequencies, and tidal phase lag, from analyzing the tide-induced WWL data.The analytical solution presented by Jeng et al. (2005) along with a set of chosen parameter values is adopted to generate the observed WWL data in the hypothetical case.Then, both synthetic and real field WWL data are analyzed to demonstrate the capability and the limitation of proposed method in the determination of the aquifer parameters and tidal characteristics.

Methodology
The method of least squares minimizes the sum of squared residuals between the observed WWL data and predicted WWL data from Jeng et al. 's solution (2005).If the observed WWL data in a specific well is analyzed, the objective function f being minimized can be expressed as where n is the number of WWL data and h p (t m ) and h o (t m ) are the predicted and observed WWL data at time t m , respectively.The SA algorithm is then applied to find the best estimates of the tidal characteristics and aquifer parameters that can minimize the objective function value.

Field data simulator
Figure 1 illustrates the groundwater fluctuations in response to tidal variations in a coastal aquifer.Jeng et al. (2005) assumed that the flow in a rigid porous medium is homogeneous and incompressible.The tides are assumed to be bichromatic, which can be synthesized by the superposition of two different wave frequencies.The water table height at the boundary of ocean and coast equals tidal oscillation (i.e., no seepage face); that is, where h is the tide-induced water table height; D is the average height of the water table, which can be regarded as the aquifer thickness; A 1 and A 2 are the amplitudes of  bichromatic-tide variations; ω 1 and ω 2 are the bichromatictide wave frequencies; δ 1 and δ 2 are the tidal phases.In addition, the boundary located at x 0 is related to the tidal oscillation as: where β denotes the beach slope.Jeng et al. (2005) gave the solution for the tide-induced water table height as where α = A 1 D is an amplitude parameter and ε = n e ω 1 D 2K is a shallow water parameter with n e and K denoting the effective porosity and hydraulic conductivity, respectively.The coefficients, H 01 , H 02 , H 11 , H 12 , and H 21 , are defined in Appendix A. Equation ( 4) is used either to calculate the values of the predicted WWL data h p (t m ) in Eq. (1) or to generate the observed WWL data h o (t m ) for the hypothetical cases.

Simulated annealing
The problem of parameter estimation involves multi-degrees nonlinear optimization and may contain several local optima in the problem domain.The SA algorithm is applied to find a set of trial solution for the unknown parameters that minimize the objective function (i.e., Eq. 1).The search algorithm starts from an initial value, referred to as the current solution TS p , and moves to a new trial solution TS p for parameter p generated by the following equation: where RD is a random number generated from a uniform (0, 1) distribution and V M p is a step length vector of the parameter p.If the trial solution is out of specified upper and lower bounds, an alternative approach for creating a new trial solution within the bounds is  Teo et al. (2003) also indicated that the shallow water parameter is usually small in real environments and suggested its value ranging from 0.1 to 0.6 in their simulations.Therefore, an additional constraint, ε < 0.6, was imposed during the search of a set of trial solutions for k/n e , ω 1 , and D to ensure that the constraint on the shallow water parameter ε is not violated.
In SA, the Metropolis criterion (1953) describes the acceptation probability of the change from the current solution i to the trial solution j .The criterion making SA have the ability to escape from local optimum is expressed as where T e represents the system temperature in SA.The acceptation probability of an inferior trial solution becomes smaller when the system cools down.
In the following, SA is used to analyze the aquifer parameters such as hydraulic diffusivity (K/n e ), beach slope and aquifer thickness as well as the tidal characteristics including tidal amplitudes, tidal wave frequencies, and tidal phase lag simultaneously.The tidal phase lag of the main harmonic constituent (δ 1 ) is excluded from the search since it is set as zero in both hypothetical and field cases.Each parameter has its own lower and upper bounds; that is, 1 m day −1 and 10 4 m day −1 for K/n e , 10 −4 and π/2 for β, 1 m and 100 m for D, 10 −4 m and 10 m for A 1 , 0 m and 10 m for A 2 , 10 −4 day −1 and 4π day −1 for ω 1 , 0 day −1 and 4π day −1 for ω 2 , and 0 and π for δ 2 .The lower bounds of A 1 and ω 1 are 10 −4 m and 10 −4 day −1 , respectively, rather than zero to prevent the denominators of λ and ω, which appear in the coefficients of Eq. ( 4) and defined in Appendix A as λ = A 2 A 1 and ω = ω 2 /ω 1 , being zero.Each parameter begins with the averaged value of the upper and lower bounds.Additionally, the SA algorithm starts at an initial temperature of 5.The system temperature reduces with a cooling rate of 0.85 after the searching of 1600 trial solutions.The algorithm terminates when the difference of the best-so-far objective functions between two consecutive temperatures is less than 10 −6 for four consecutive times or the iteration number exceeds 2 × 10 7 .

Generation of hypothetical data
Six hypothetical scenarios are designed herein to demonstrate the capability and limitation of proposed method in the determination of the aquifer parameters and tidal characteristics.In each scenario, five cases named from cases a to e are analyzed.Case a contains the noise-free WWL data generated by Eq. ( 4), while other cases have the data with measurement errors produced by where h ub (x,t m ) denotes the noise-free WWL data; ϕ is chosen to be 1 % and represents the accuracy of field measurements on the order of centimeter; RN (m) represents the random number on the order m.Four sets of random numbers being adopted in cases b to e are normally distributed and generated by the routine RNNOF of IMSL (2003).
3 Results and discussion

Analysis of hypothetical data
Table 1 shows the estimated results form the analyses of three scenarios for synthetic WWL data.Scenarios 1-3 respectively represent the problem that the WWL data being observed at the well located at 5 m, 10 m, and 20 m away from the intersection point of the beach surface and mean sea level.
Table 1 provides the root mean squared error (RMSE) between the predicted and observed WWL data as an index to examine the accuracy of the prediction, which is defined as Moreover, the standard deviation (SD), relative error (RE), 95 % lower-limit and upper-limit of confidence interval (95 % LLCI and 95 % ULCI) are also given in Table 1.
The LLCI and ULCI are calculated using the formula ȳ ± s ȳ t 4,0.025where ȳ is the mean value of estimated parameters from cases a to e; s ȳ is the estimated standard error of the mean; t 4,0.025 is t statistic with degrees of freedom equaling 4 and 95 % confidence interval obtained from a t-distribution table as 2.776.In Table 1, the WWL data were generated based on the parameters K/n e , β, D, A 1 , A 2 , ω 1 , ω 2 and δ 2 being given as 500 m day −1 , π/3 (1.047), 25 m, 2 m, 1 m, 4πday −1 , 2πday −1 , and π/4, respectively.The estimated results are fairly close to the target values.In addition, the target values of estimated parameters are all within their 95 % confidence interval and have small RMSE values on the order of 10 −4 to 10 −3 .These results indicate good matches between the predicted and synthetic WWL data.The CPU time listed in the right column is the computing time for parameter estimation in each case using a personal computer with Genuine Intel CPU 2140 @ 1.60 GHz and 1 GB RAM.
The searches for the optimal results for all the cases are done within two minutes.Table 2 provides the correlation matrix of the eight estimated parameters in scenarios 1-3.Only the

Estimated results
Aquifer parameters Tidal characteristics Table 3 examines the effect of various ε on the estimated results.Scenarios 4 and 5 have the same target parameter values and well location as scenario 2 except that K/n e become 50 m day −1 and 5000 m day −1 representing the cases of ε being 1.772 and 0.177, respectively.For the cases with an extremely large ε like scenario 4, the estimated results show large deviations on most of the estimated parameter values and the RMSEs are on the order of 10 −1 m.The target values of K/n e , ω 1 and δ 2 , especially K/n e , are out of their 95 % confidence interval.In scenario 4, the algorithm was terminated when the difference of the best-so-far objective function value between two consecutive temperatures was less than 10 −6 for four consecutive times.The total iteration numbers didn't exceed the maximum iteration number allowed in the algorithm.Such a poor result is likely due to the constraint imposed on the shallow water parameter  during the parameter search in SA, which makes it impossible to find out the target parameters.On the other hand, in scenario 5 with a small ε, the proposed method gives the estimations of the aquifer parameters and tidal characteristics close to their target ones.Only the target value of β is out of its 95 % confidence interval.
Table 4 shows the estimated parameter values based on SA; yet, Nielsen's solution (1990), in lieu of Jeng et al. 's solution (2005), was adopted to fit the WWL data.Scenario 6 in Table 4 has the same synthetic WWL data, which were generated by Jeng et al. 's solution (2005), as scenario 2. Nielsen's solution (1990) considers the monochromatic-tide effect and is expressed as where k is the wave number defined as k = ε/D and ξ is the perturbation parameter given by ξ = kA 1 cot(β).Note that the constraint on the shallow water parameter ε to be less than 0.6 is released in scenario 6.The estimated values of A 1 and ω 1 are fairly close to the target values of main harmonic constituent of bichromatic tide.However, the estimated values of K/n e and β are larger than their target values.Figure 2 shows the synthetic heads and predicted heads in scenarios 2 and 6.The pattern of the predicted WWL data in scenario 6 can be depicted by a single cosine function.In contrast, the predicted WWL data in scenario 2 has the period characteristic of each monochromatic tide as well as a new period of resultant tide.This figure indicates that the predicted heads in scenario 2 match well with the synthetic heads.However, the predicted heads in scenario 6 significantly differ from the synthetic heads because Nielsen's solution (1990) only considers the monochromatic-tide effect.Naturally, Jeng et al. 's solution (2005), which considers additional harmonic constituents, is more flexible to fit complicated field data than Nielsen's solution (1990).

Analysis of field data
Nielsen (1990) presented a hydraulic head solution in a coastal aquifer affected by the tide fluctuations and analyzed the WWL measurements observed at Barrenjoey beach in Australia.He provided the WWL at wells 1 to 11 over 26 h in his Table 1.Wells 7-11 were located at 6.6 m, 9.1 m, 11.6 m, 14.1 m, and 16.6 m, respectively, away from the intersection of the beach surface and mean sea level.The aquifer thickness D was 0.51 m.The tidal amplitudes A 1 and A 2 were 0.516 m and 0.014 m, respectively, while the tidal wave frequencies ω 1 and ω 2 were 0.513 h −1 and 1.026 h −1 , respectively.Moreover, the tidal phases δ 1 and δ 2 were 0 and 1.303, respectively.The wave number k was determined by taking the average of the regression results of k in damping term and phase lags term in those wells to the hydraulic head solution.Thereafter, the value of K/n e was obtained as 2076 m day −1 after inserting the known values of k, ω 1 and D to the formula of wave number.The estimated value of β was 0.09967, which was calculated from the difference in sand-level heights at wells 1 and 11 and divided by the distance between these two wells.
Table 5 shows the estimated results of aquifer parameters from the analysis of WWL data measured at Barrenjoey beach.The tidal characteristics are assumed to be the same as those given in Nielsen (1990).The lower and upper bounds of each parameters used in SA are the same as those mentioned in Sect.2.2, except that the lower bound of D is reset as 0.1 m.The estimated values of K/n e based on the proposed method at all wells is significantly smaller than that given in Nielsen (1990).A composite analysis for the WWL data collected from the five wells is also performed.The estimated values of k/n e , β, and D obviously differ from those obtained from the single-well WWL data analysis or those given in Nielsen (1990).The estimated results based on the proposed method with smaller RMSE values for the three aquifer parameters in either single-well or multi-well analysis might be better than those given in Nielsen (1990).However, it is rather difficult to identify that at which well the WWL data analysis gives the most reasonable estimates of K/n e , β, and D.
Figure 3 shows the observed WWL given in Nielsen (1990) represented by open circles.The figure also shows the predicted WWL produced based on Nielsen's parameters and Table 2. Correlation matrix of estimated parameters in scenarios 1-3.The RMSE values calculated with the parameters from our proposed method.b The RMSE values calculated with the parameters given in Nielsen (1990).solution (1990) are plotted as dash-dot-dash lines, while the predicted WWLs produced by Jeng et al. 's solution (2005) with the parameters determined by the present method with single-well and multi-well data are drawn as solid and dashed lines, respectively.The figure indicates that the magnitudes of WWL fluctuations decrease landward.At each well, it seems that both solid and dashed lines are closer to the observed WWL data (in circle) than the dash-dot-dash line, demonstrating the capacity of the present method in the estimation of tidal characteristics and aquifer parameters.On the other hand, the greatest discrepancies between the real WWL data and predicted WWL data occur at low tide, i.e., around t = 8 h to t = 10 h.Nielsen (1990, Fig. 7) also addressed the same problem when comparing the real data (also used in this study) with the data generated by his solution.He explained that the discrepancy might be due to the boundary condition at x 0 , where the assumption of no seepage face contradicts the reality.

Conclusions
The  Nielsen (1990).The estimated results of aquifer parameters and tidal characteristics from the present method are fairly close to those of target ones in the analysis of synthetic WWL data.In addition, Jeng et al. 's solution (2005) gives good predictions of the WWL fluctuations induced by the bichromatic tide in sloping coastal aquifer in the shallow-water cases.
When analyzing the field WWL data, the aquifer parameters are estimated alone by setting the tidal characteristics as known.The estimated aquifer diffusivity and aquifer thickness based on single-well and multi-well analyses are obviously different from those given in Nielsen (1990).The comparisons of the observed and predicted WWL data show that the present method gives better predictions to the observed WWL data than the predicted results based on the solution and parameter values given in Nielsen (1990).

21Figure 1 .
Figure 1.Illustration of water table fluctuations in response to tidal variation in a coastal aquifer.

Fig. 1 .
Fig. 1.Illustration of water table fluctuations in response to tidal variation in a coastal aquifer.

Fig. 2 .
Fig. 2. Comparisons of synthetic heads and predicted heads in scenarios 2 and 6.The synthetic heads in scenario 2 are analyzed based on Jeng et al.'s solution (2005) and those in scenario 6 are analyzed based on Nielsen's solution (1990).

Fig. 3 .
Fig. 3. Plots of observed WWL given inNielsen (1990), predicted WWL produced by the solution and parameters inNielsen (1990), and predicted WWLs produced by Jeng et al.'s solution (2005)  with the parameters determined by the present method with single-well and multi-well data.

Table 1 .
Estimated results form synthetic WWL data.Scenarios 1-3 denote the wells located at x = 5, 10, and 20 m, respectively.The target values of the parameters are K/n e = 500 m day −1

Table 3 .
Estimated results form synthetic WWL data.Scenarios 4 and 5 have the same target parameter values and well location as scenario 2 except that K/n e becomes 50 m day −1 when ε = 1.772 and 5000 m day −1 when ε = 0.177, respectively.

Table 4 .
The results estimated based on Nielsens solution (1990) with the synthetic WWL data generated from Jeng et al.'s solution (2005).

Table 5 .
Nielsen (1990)lts using the aquifer parameters fromNielsen (1990)and the proposed method based on the WWL data at Barrenjoey beach in Australia.
method, based on coupling Jeng et al.'s solution (2005) with SA algorithm, is developed to simultaneously estimate the hydraulic diffusivity, beach slope, aquifer thickness, tidal amplitudes, tidal wave frequencies, and tidal phase lag.The method is used to analyze the synthetic WWL data generated Hydrol.Earth Syst.Sci., 15, 1473-1482, 2011 www.hydrol-earth-syst-sci.net/15/1473/2011/ by Jeng et al.'s solution (2005) and the field WWL data, collected from Barrenjoey beach in Australia, presented in