China Ocean Eng., Vol. 29, No. 4, pp. 489 – 502 © 2015 Chinese Ocean Engineering Society and Springer-Verlag Berlin Heidelberg DOI 10.1007/s13344-015-0034-y, ISSN 0890-5487

Uncertainty and Sensitivity Analyses of the Simulated SeawaterFreshwater Mixing Zones in Steady-State Coastal Aquifers* ZHAO Zhong-wei (赵忠伟)a, 1, ZHAO Jian (赵

坚)b, XIN Pei (辛

沛)a, b,

HUA Guo-fen (华国芬)b and JIN Guang-qiu (金光球)a, b a

State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering, Hohai University, Nanjing 210098, China

b

College of Water Conservancy and Hydropower Engineering, Hohai University, Nanjing 210098, China (Received 20 March 2015; received revised form 18 May 2015; accepted 24 May 2015)

ABSTRACT The uncertainty and sensitivity of predicted positions and thicknesses of seawater-freshwater mixing zones with respect to uncertainties of saturated hydraulic conductivity, porosity, molecular diffusivity, longitudinal and transverse dispersivities were investigated in both head-control and flux-control inland boundary systems. It shows that uncertainties and sensitivities of predicted results vary in different boundary systems. With the same designed matrix of uncertain factors in simulation experiments, the variance of predicted positions and thickness in the flux-control system is much larger than that predicted in the head-control system. In a head-control system, the most sensitive factors for the predicted position of the mixing zone are inland freshwater head and transverse dispersivity. However, the predicted position of the mixing zone is more sensitive to saturated hydraulic conductivity in a flux-control system. In a head-control system, the most sensitive factors for the predicted thickness of the mixing zone include transverse dispersivity, molecular diffusivity, porosity, and longitudinal dispersivity, but the predicted thickness is more sensitive to the saturated hydraulic conductivity in a flux-control system. These findings improve our understandings for the development of seawaterfreshwater mixing zone during seawater intrusion processes, and give technical support for groundwater resource management in coastal aquifers. Key words: seawater intrusion; mixing zone; uncertainty analysis; fractional factorial design; Morris’s OAT design; coastal aquifer

1. Introduction The mixing zone between seawater and freshwater is an important feature of coastal aquifers (Michael et al., 2005). It induces the submarine groundwater discharge (SGD) whose magnitude might be several times larger than that of river system (Wang et al., 2015). The mixing zone transports nutrients, contaminants, and other chemicals to the coastal waters, and has significant impact and implication on coastal environment and ecology (Li, 2013). Its behavior underpins salt exchange between the ocean and coastal aquifer and affects fates of reactive solutes (Simmons, 1992; Moore, 1996). Effective management of coastal groundwater resources requires that the mixing zone is properly characterized and predicted (Werner et al., 2013). Numerical simulation is an effective research tool for predicting a seawaterfreshwater mixing _____________________________________________________

* This project was financially supported by the National Natural Science Foundation of China (Grant Nos. 51309091, 51239003 and 51279045), and the Postdoctoral Science Foundation of China (Grant No. 2012M520989). 1 Corresponding author. E-mail: [email protected] 

490

ZHAO Zhong-wei et al. / China Ocean Eng., 29(4), 2015, 489 – 502  

process (Zheng et al., 2002). However, a numerical model is never a perfect representation of reality (Sanford and Pope, 2010). Uncertainties exist within various components of the model, caused by the simplification in model’s conceptualization (Xia and Li, 2012; Li and Jiao, 2013), the measurement error of input parameters and flow velocities (Ma et al., 2015; Cheng et al., 2012), the spatial variance of initial conditions, and the spatial and temporal variance of boundary conditions (Bear and Cheng, 2010). All these uncertain factors can cause large biases to the model predictions, deviating from the true value. Uncertainty analysis gives technical supports for decision-making in groundwater resources management through quantification of the uncertainties in the predicted mixing zones. Sensitivity analysis investigates the relationships between input and output variables in numerical models, and gives instructions for model simplification and calibration. It improves our understanding on behaviors of the mixing zone during seawater intrusion processes and enables proper assessment of coastal groundwater resources with consideration of risks linked to uncertainties (Zhao et al., 2013). The heterogeneity of aquifers is supposed to be one of the most important factors controlling the development of saltwater wedges and attracts much attention in uncertainty analysis of seawater intrusion simulations (Diersch and Kolditz, 2002; Simmons, 2005; Pool and Carrera, 2011). Based on the sharp interface model, Dagan and Zeitoun (1998) and Naji et al. (1999) gave analytical solutions for variance of seawaterfreshwater interfaces in coastal aquifers with randomly distributed permeability. Using Monte Carlo simulation experiments, Lecca and Cau (2009), Kerrou et al. (2010) and Herckenrath et al. (2011) conducted uncertainty analysis on simulated mixing zones caused by heterogeneity of hydraulic conductivities in coastal aquifers. However, the uncertainties of other parameters have not been well quantified even though they were reported to play important roles on predicted mixing zones, such as longitudinal and transverse dispersivities (Rajabi and Ataie-Ashtiani, 2014). Moreover, the contributions of input factors to uncertainty of predicted results vary with different boundary conditions (Smith, 2004; Werner and Simmons, 2009; Werner et al., 2012). For example, a head-controlled system with inland boundary conditions defined as a special head (a lake or a river in the real system) or a flux-controlled system with inland boundary conditions defined as a specified freshwater influx. Therefore, uncertainty and sensitivity analyses of predicted mixing zones should also investigate the differences caused by different boundary systems. A seawaterfreshwater mixing zone is normally characterized by its position and thickness which are defined with contour lines of 10%, 50% and 90% salt concentrations (Lu et al., 2013). Observed thicknesses of mixing zones vary widely in laboratory experiments, field investigations and numerical simulations (Werner et al., 2013). Based on inverse analysis of Henry’s model, Sanz and Voss (2006) concluded that the thickness of a mixing zone was mainly controlled by porosity, molecular dispersivity and inland freshwater input. However, the hydrodynamic dispersion was not concerned in their model. Abarca et al. (2007) showed that longitudinal and transversal dispersions are equally important in controlling the thickness of a mixing zone. Paster and Dagan (2007) suggested that the steady state thickness of a mixing zone was dominated by local transverse dispersion which explained the narrow mixing zones observed in laboratory experiments. However, field investigations revealed wide mixing zones varying from several meters to kilometers (Xue et al., 1993; Paster et al., 2006). There is a tremendous need to investigate the uncertainty and sensitivity of the predicted thickness of  

ZHAO Zhong-wei et al. / China Ocean Eng., 29(4), 2015, 489 – 502 

491

mixing zones in different boundary systems. Based on a hypothetical coastal sandy aquifer, an uncertainty and sensitivity analysis on the predicted seawaterfreshwater mixing zones was made in a head-control and a flux-control system respectively, in consideration of the uncertainties of saturated hydraulic conductivity, porosity, molecular diffusivity, and longitudinal and transverse dispersivities. The simulated mixing zone was characterized by its position (S) and thickness (D) at a depth (as defined in the Appendix). A group of simulation experiments were designed based on the Fractional Factorial Design (FFD) method in order to quantify the uncertainty of the predicted positions and thicknesses of seawaterfreshwater mixing zones. Sensitivity analysis was carried out based on simulation experiments designed by the Morris’s One-at-a-time (Morris’s OAT) method. The numerical model, FFD and Morris’s OAT methods are introduced in Section 2. Uncertainty and sensitivity analyses results are described in Section 3.

2. Methodology 2.1 Numerical Model The variably saturated variable-density groundwater flow and transport modeling package SUTRA (Voss and Provost, 2010) was used in the simulation experiments. The compressibility of soil and fluid was neglected (Xin et al., 2009). The salt adsorption of the sand was assumed to be negligible (Kuan et al., 2012). As a coarse sandy aquifer was investigated, the capillary effect is weak and the unsaturated groundwater flow is negligible compared with the saturated groundwater flow. Therefore, a two-dimensional saturated variable-density flow is employed and described as follows (Voss and Provost, 2010):  ( )     v  , t v

Ks



 p    z , g   

(1a) (1b)

where ε is the porosity of solid matrix; ρ is the density of fluid that varies with salt concentration  C with ρ0 being the freshwater density, C being the salt concentration and according to    0  C  (  714.3 kg/m3) describing the linear relationship between fluid density and salt concentration; t C is time; v is the average velocity; Ks is the saturated hydraulic conductivity; p is the pressure; g is the magnitude of gravitational acceleration; and z is the elevation. The governing equation of solute transport based on mass balance is shown below:    C  t

     vC       Dm I + D   C  ,

(2)

where Dm is the molecular diffusivity; I is the identity tensor; D is the dispersion tensor whose element is given by Di , j  f  L ,  T , v, vi , j  , (i, j = x, y, z) with αL being the longitudinal dispersivity, αT being

the transverse dispersivity, v being the magnitude of velocity and vi,j being the magnitude of velocity  

492

ZHAO Zhong-wei et al. / China Ocean Eng., 29(4), 2015, 489 – 502  

component. A hypothetical coastal sandy aquifer was investigated in simulation experiments. The scale of the aquifer was only hundreds of meters. The aquifer was assumed to be homogeneous with heterogeneity characterized by uncertainty of hydraulic conductivity and porosity. A two-dimensional model was built in cross-shore direction, assuming that the flow and solute transport in alongshore direction were negligible (Fig. 1).

A (m)

B (m)

C (m)

D (m)

(100, 3)

(50, 3)

(50, 3)

(500, 3)

E (m)

F (m)

(500, 30) (100, 30)

Fig. 1. Model geometry: the coordinates of domain reference points are shown in the table (x, z).

The model was 600 m long from the vertical section of seaward boundary (AF) to inland boundary (DE). The top surface of the model (CD) was at a height of 3 m above the mean sea level. An impermeable base (EF) was at depth of 30 m below the mean sea level. A sloping beach (BC) was adopted for the foreshore section of seabed. The slope extended offshore with a horizontal seabed platform (AB) at the depth of 3 m below the mean sea level. The inland boundary (DE) condition was defined as a specified freshwater influx (Q) in a flux-controlled system (Qs) or a specified freshwater head (H) in a head-controlled system (Hs) with zero solute concentration. The boundary conditions for beach (BC) and seabed platform (AB) under the sea level were defined as a specified head (H) with salt concentration 35 ppt (mass fraction, parts per thousand). No area recharge or evapotranspiration was considered and hence the top surface (CD) was specified as no water and no solute flux boundary. The aquifer base (EF) and the vertial part of seaward boundary (AF) were set to be impermeable with no water and solute flux. Table 1 Parameter Ks (m/s) αL (m) αT (m) ε Dm (m) Q (m2/day) H (m)

Input parameters for models Based model

Uncertainty model

1.157104 1.00 0.1 0.30 1.0109 1.65 2.5

1.157104 5.787104 0.251.00 0.0250.1 0.300.45 1.0108 1.01010 1.651.98 2.53.0

Two base models were set up to represent two “real conditions” with the same deterministic input parameters, one for a head-control system and the other for a flux-control system (Table 1). The two models ran for 900 days until reaching steady states, both starting with initial zero hydraulic head and solute concentration based on Glover’s (1959) sharp interface solutions. The steady state of the  

ZHAO Zhong-wei et al. / China Ocean Eng., 29(4), 2015, 489 – 502 

493

pressure and solute distributions of the two models were proven to be closely consistent and subsequently taken as initial conditions for simulation experiments with the same inland boundary conditions, respectively. The Péclet numbers Pe ≤ 4 (Voss and Provost, 2010) and Courant numbers Cr< 1 (Robinson et al., 2007) were satisfied to avoid numerical dispersion and oscillations. Tests were conducted to make sure that the simulations met the convergence criteria and were independent of time steps and mesh sizes. All simulations were run until both hydraulic head and salt concentration solutions reached a quasi-steady state (Xin et al., 2010). 2.2 Fractional Factorial Design Factorial Design (FD) is an experiment-design method aimed at estimating the contribution of each factor to behaviors of a system, including the effects of each factor and their interplay (Box et al., 1978). It is recognized as an effective screening tool for selection of important factors in simulation experiments (Saltelli et al., 2000). In FD, the experiments are designed with full combinations of possible input values at different levels. The number of simulation experiments increases geometrically with the number of input factors, which creates a large burden for uncertainty analysis of seawater intrusion simulations. FFD assumes some higher order interplay to be unimportant and negligible. Thus, the number of evaluated combined effects is effectively reduced compared with that in a full FD. The number of simulation runs is reduced significantly as well, especially when the number of uncertain factors is very large. Two-level FFD is used for uncertainty analysis. The main effects of each variable and interplay with other variables are estimated by the difference between averaged responses at the two levels. The total effect of an input factor on the response includes the main effect and the interplay involving it. The main effect, interplay and total effect of a factor give quantitative evaluations of relative importance of the factor contributing to the total variance of the predicted results, with respect to single effect, interplay with other factors and total effect respectively. Thirty-two simulation experiments were designed with FFD in each boundary system to investigate the uncertainty of the predicted positions and thicknesses of the mixing zone and the contributions of the uncertain factors. 2.3 Morris’s OAT Design Morris’s OAT design is a global sensitivity analysis method which covers the entire value space of a factor (Morris, 1991). It estimates the sensitivity of a factor by computing a number of local measures at different evaluated points of the input space and then takes their average. It reduces dependence of the sensitivities on specific points which local experiments suffer from (Saltelli et al., 2000). Moreover, the Morris’s OAT design method is more computationally efficient, requiring few simulation runs to obtain a similar result. In Morris’s OAT design, the elementary effect of a factor (xi, i=1, 2, …, k) is used to characterize the variation of the simulated result (y) caused by a unit change of the factor (Morris, 1991). In order to compare the effects of different results with respect to a single input factor or effects of a single result with respect to different input factors, a scaled elementary effect is introduced as follows:

di ( x )   

y  x1 , , xi 1 , xi  , xi 1 , , xk   y  x  

  xi ,max  xi ,min  ,

(3)

494

ZHAO Zhong-wei et al. / China Ocean Eng., 29(4), 2015, 489 – 502  

where Δ is the perturbation at the point x=(x1,…, xi, …, xk), (i=1, 2, …, k), xi,max and xi,min are the maximum and minimum values of xi, respectively. The scaled elementary effect provides a possible way for comparing the importance of a single result with respect to different input factors whose values might be different in orders of magnitude, such as the hydraulic conductivity and longitudinal dispersivity. Morris’s OAT design method is based on the construction of a randomly selected orientation matrix B* which provides one elementary effect per factor. The orientation matrix B* is given by: 1   B*   J k 1,1 x *   2 B  J k 1, k  D*  J k 1, k  Δ  P * , (4) 2   where Jk+1,1 is a (k+1)×k dimensional matrix with all elements equal to “1”; x*=[x1, x2, …, xi, …, xk] is a randomly selected base value; B is a specified (k+1)×k dimensional lower triangular matrix with element “1” when the index of the row is bigger than the index of the column; D* is a k-dimensional

diagonal matrix whose diagonal element is “1” or “1” with equal probability; P* is the matrix modified from a k-dimensional unit matrix with rows randomly exchanged. The variance of the factor xi is divided into p levels (p is even) and the perturbation of xi is defined as: p  ( xi ,max  xi ,min ) . (5) 2( p  1) The orientation matrix B* is sampled r times independently with different starting points to obtain r (number of) evaluations of the elementary effects for each factor. Thus, rk elementary effects are evaluated with r×(k+1) simulation runs, which normally require 2rk simulation runs. The computational efficiency (defined by the number of elementary effect divided by the number of the simulation runs) is increased from 1/2 to k/(k+1). Morris’s OAT method was employed to design simulation experiments for sensitivity analysis of the predicted mixing zone. The variances of input factors are divided into 2 levels (p = 2), and the orientation matrix is sampled 4 times (r=4) with perturbation Δ=1/3. Twenty-eight experiments were designed in each boundary system to investigate the sensitivities of the simulated mixing zone with respect to the uncertain factors.

3. Results and Discussions 3.1 Uncertainty of the Predicted Mixing Zones and Contributions of Uncertain Factors The value ranges of the predicted S and D and the contributions of the uncertain factors are closely related to inland boundary conditions (Figs. 2 and 3). The value range of the predicted S at the

bottom (z=30 m) is [52.1, 78.1] (unit: m) in the head-control system. More uncertainty exists for the predicted S in the flux-control system compared with the results in the head-control system, with the value range of [44.7, 302.6] (unit: m) at the bottom, nearly nine times larger than that in the head-control system. The value range of predicted D at the bottom (z=30 m) is [0.8, 6.2] (unit: m) in the head-control system. While this range in the flux-control system is [0.8, 16.2] (unit: m), about 10 meters larger than that in the head-control system. In the head-control system (Hs), the most important factors underpinning the uncertainty of the predicted S at the bottom include freshwater hydraulic head, transverse dispersivity, and saturated  

ZHAO Zhong-wei et al. / China Ocean Eng., 29(4), 2015, 489 – 502 

495

hydraulic conductivity (ordered by the magnitudes of total effects as shown in Table 2). The saturated hydraulic conductivity has a positive effect on the predicted S (the predicted result increases with the increasing input parameter), while the freshwater hydraulic head and transverse dispersivity have negative effects on the predicted S (the predicted result decreases with the increase of input parameter). Most contributions of the total effects for the freshwater hydraulic head and the transverse dispersivity are given by the main effects of themselves (Fig. 4a). The contributions of the total effect for the saturated hydraulic conductivity are mostly given by its main effect and the interactions with freshwater hydraulic head, porosity, and longitudinal and transverse disperivities (Fig. 4a).

Fig. 2. Contour lines of 50% salt concentrations in head-control and flux-control systems. Note: The base cloud is the result of the base model; the black and white lines in the head-control system are the contour lines of 50% salt concentrations for cases whose hydraulic head of the inland boundaries were designed with the upper and lower level, respectively; the yellow and white lines in the flux-control system are the contour lines of 50% salt concentrations for cases whose saturated hydraulic conductivity were designed with the upper and lower level, respectively; the contour lines near the depth z = 30 m (in rectangular region) are zoomed in.

Fig. 3. Thickness of seawaterfreshwater mixing zone at different depths. Note: The solid line stands for uncertainty range for the thickness of the mixing zone; the dot stands for the thickness of the mixing zone in base models; the asterisk stands for the predicted results in simulation experiments.

Table 2 Total effects on the predicted S and D at the bottom (z=30 m) Uncertain factor Ks αL αT ε Dm H or Q

S Hs 16.6 5.0 55.2 2.0 1.0 155.0

D Qs 1962.2 11.4 345.0 13.8 36.2 452.6

Hs 8.2 7.3 26.2 0.8 1.1 5.2

Qs 52.4 5.5 42.3 4.2 14.1 11.3

In the flux-control system (Qs), the most important factors contributing to total uncertainty of predicted S at the bottom include saturated hydraulic conductivity, freshwater influx, and transverse dispersivity (ordered by magnitudes of total effects as shown in Table 2). Compared with results in the head-control system, the saturated hydraulic conductivity plays a more important role on uncertainty of  

ZHAO Zhong-wei et al. / China Ocean Eng., 29(4), 2015, 489 – 502  

496

predicted S. Similarly to results in the head-control system, the predicted S increases with the increasing saturated hydraulic conductivity, and decreases with the increase of the freshwater influx and transverse dispersivity (Table 2). Most contributions of total effects for the freshwater hydraulic conductivity, the freshwater influx, and the transverse dispersivity come from the main effects of themselves and the interplay with each other (Fig. 4b). Ks αL αT ε Dm H

13.5 5.5 3.9 4.50 1.9 4.9 Kf

12.1 2.9 3.10 0.1 1.9 αL

65.9 0.70 2.1 8.9 αT

5.50 0.70 2.10 ε

5.7 1.5 Dm

164.5 H

Ks αL αT ε Dm Q

(a) Effect on S in head-control system Ks αL αT ε Dm H

4.9 1.8 2.2 1.14 3.4 0.8 Kf

8.3 0.4 0.54 0.2 0.3 αL

26.2 0.50 1.9 0.2 αT

1.64 1.05 0.74 ε

5.4 0.3 Dm

(c) Effect on D in head-control system

2345.3 5.1 169.9 6.30 23.5 178.3 Kf

9.9 3.7 221.9 1.50 5.30 8.30 16.5 0.7 3.50 2.1 21.3 0.50 αL αT ε

29.1 4.1 302.3 Dm Q

(b) Effect on S in flux-control system

4.5 H

Ks αL αT ε Dm Q

41.2 1.3 10.1 0.66 4.3 2.6 Kf

5.4 1.2 0.36 0.2 0.4 αL

35.1 2.03 0.37 1.98 2.6 1.1 0.51 αT ε

11.2 1.0 Dm

5.7 Q

(d) Effect on D in flux-control system

Fig. 4. Main effects of the uncertain factors on predicted S and D at the bottom (z = 30 m) and their interplays.

The result of uncertainty analysis gives an assessment on the value range of predicted S and D, as well as contributions of uncertain factors to uncertainty of predicted results. The contributions are measurements by relative importance of the uncertain factors to uncertainty of the predicted results. As the predicted results between the head-control and flux-control systems varied significantly, the total effect, main effect and interplay of an input factor may differ in orders of magnitude in the headcontrol and the flux-control systems. Therefore, it is not reasonable to compare values of the total effect, main effect and interplay of the same uncertain factor in different inland boundary systems. 3.2 Sensitivity of the Predicted Mixing Zone with Respect to Uncertain Factors The sensitivities of the predicted results with respect to uncertain factors were evaluated by scaled element effects of input factors on the predicted results. The results show that only the saturated hydraulic conductivity has positive element effects on the predicted S in both head-control and flux-control systems (Figs. 5a and 5b), which means that the predicted S increases with the increase of saturated hydraulic conductivity but decreases with the increase of other factors. In a head-control system, the magnitudes of scaled elementary effects for saturated hydraulic conductivity are quite small (Fig. 5a). It indicates that the predicted S is almost independent of saturated hydraulic conductivity. The scaled elementary effects on the predicted S in a flux-control system are much larger than that in a head-control system (Fig. 5b). It means that the saturated hydraulic conductivity plays a more important role to the predicted S in a flux-control system than that in a head-control system. The predicted D decreases with the increase of the saturated hydraulic conductivity, but increases with the  

ZHAO Zhong-wei et al. / China Ocean Eng., 29(4), 2015, 489 – 502 

497

increase of other factors in a head-control system (Fig. 5c). However, the predicted D increases with all input factors in a flux-control system (Fig. 5d).

Fig. 5. Scaled elementary effects (x) on predicted S and D at different depths (z) in head-control and flux-control systems.

The scaled elementary effects of input factors on the predicted S and D vary with depth. It shows that the scaled elementary effects of the saturated hydraulic conductivity on the predicted S increase with the increasing depth in both head-control and flux-control systems. However, the scaled elementary effects of other factors decrease with the increasing depth in both systems (Figs. 5a5d). It indicates that the predicted S becomes more sensitive to saturated hydraulic conductivity with the increase of depth. While, the predicted D becomes more sensivity to other factors in a head-control system. It is reasonable to compare the scaled elementary effects of different uncertain factors on the predicted results at the same depth. The scaled elementary effects of uncertain factors on the predicted S and D at depth z = 25 m are selected for sensitivities comparison of the predicted results with respect to uncertain factors (Table 3). Table 3

Scaled elementary effects on the predicted S and D (z = 25 m)

Uncertain factor Ks αL αT ε Dm H or Q

 

S Hs 0.2 0.3 3.7 0.8 0.8 10.8

D Qs 136.3 0.6 15.8 4.3 5.7 27.8

Hs 0.2 0.5 1.6 0.5 0.7 0.2

Qs 2.9 0.4 2.2 0.7 1.2 0.1

498

ZHAO Zhong-wei et al. / China Ocean Eng., 29(4), 2015, 489 – 502  

The sensitivity of the predicted S at depth z = 25 m with respect to uncertain factors varies in different inland boundary systems. The predicted S at this depth is more sensitive to freshwater hydraulic head and transverse dispersivity in the head-control system, but not sensitive to the saturated hydraulic conductivity. While the most sensitive factor for the predicted S at this depth is the saturated hydraulic conductivity in the flux-control system, and then followed by freshwater hydraulic head and transverse dispersivity. Similarly as uncertainty analysis results, the predicted S increases with the increasing saturated hydraulic conductivity and decreases with the increase of other uncertain factors both in the head-control and the flux-control systems. The sensitive factors for the predicted D at depth z = 25 m in the head-control system include transverse dispersivity, molecular diffusivity, porosity, and longitudinal dispersivity (with the decreasing magnitude of scaled elementary effects), without the saturated hydraulic conductivity. However, saturated hydraulic conductivity is the most sensitive factor for the predicted D at this depth in the flux-control system, and then followed by transverse dispersivity, molecular diffusivity, porosity, and longitudinal dispersivity. The four sensitive factors have positive effect on the predicted D both in the head-control and the flux-control systems, which means that the predicted D increases with all the four sensitive factors. As known from the GhybenHerzberg relationship, the steady state position of sharp interface between seawater and freshwater at a depth is mainly controlled by the hydraulic gradient between the inland and seawater boundaries. The steady state hydraulic gradient is mainly controlled by the hydraulic head between the inland and seaward boundaries in a head-control system, not too much influenced by the geometric characters of the aquifer, such as hydraulic conductivity. Therefore, the steady state position of sharp interface in a head-control system is not influenced by the hydraulic conductivity. However, the hydraulic gradient is strongly affected by the saturated hydraulic conductivity and influx in a flux-controlled system. Therefore, the steady state position of the mixing zone is more sensitive to the saturated hydraulic conductivity in the flux-controlled system compared with that in a head-control system. In steady state homogeneous aquifers, the thickness of a mixing zone is affected by longitudinal and transverse dispersions, molecular diffusion, freshwater influx, and the density gradient between seawater and freshwater (Volker and Rushton, 1982). As revealed in sensitivity analysis, transverse dispersivity is the most sensitive factor affecting the predicted thickness of the mixing zone in a head-control system. Other sensitive factors include the molecular diffusivity, longitudinal dispersivity, and porosity. The predicted thickness of the mixing zone is more sensitive to the saturated hydraulic conductivity than the transverse dispersivity in a flux-control system. Therefore, sensitivity analysis of the predicted thickness of a mixing zone should also consider the interactions between the aquifer properties and hydraulic boundary conditions.

4. Conclusions Uncertainty and sensitivity analyses for the predicted seawaterfreshwater mixing zones are important processes for decision makers in groundwater resources management. This study made an  

ZHAO Zhong-wei et al. / China Ocean Eng., 29(4), 2015, 489 – 502 

499

assessment on uncertainties of the predicted positions and thicknesses of seawaterfreshwater mixing zones caused by the uncertainty of the saturated hydraulic conductivity, porosity, molecular diffusivity, and longitudinal and transverse dispersivities, and evaluated the sensitivities of predicted results with respect to these uncertain factors. Two different inland boundary systems, a head-control system and a flux-control system, were investigated. The contributions of the uncertain factors to uncertainty of predicted position and thickness of the mixing zone were determined by the total effects of the uncertain factors. The sensitivities of the predicted results to the uncertain factors were measured by the scaled elementary effects of the uncertain factors on the uncertainty of the predicted results. The uncertainty of the predicted positions and thicknesses of the mixing zone varies with different boundary conditions. With the same designed matrix of uncertain factors, the value range of the predicted positions in the flux-controlled system is much larger than that predicted in the headcontrolled system. The contributions of the uncertain factors to the uncertainty of the predicted position and thickness of the mixing zone are closely related to boundary conditions as well. In the head-control system, the most important factor contributing to the uncertainty of the predicted position of seawaterfreshwater mixing zone is the hydraulic gradient between the inland and seaward boundary. While the most important factor contributing to the uncertainty of the predicted position of seawaterfreshwater mixing zone changed to the saturated hydraulic conductivity in the flux-control system. It is necessary to investigate whether the boundary conditions are properly defined and have no influence on the uncertainty analysis results. Otherwise the uncertainty analysis should consider uncertainties caused by different boundary systems. The inland boundary condition also influences the sensitivities of the predicted positions and thicknesses of the mixing zones with respect to input factors. According to the evaluated scaled elementary effects at depth z = 25 m, the most sensitive factors for the predicted position of the mixing zone at this depth include inland freshwater hydraulic head and transverse dispersivity in the head-control system. However, the predicted position at this depth is more sensitive to the saturated hydraulic conductivity and then followed by freshwater influx and transverse dispersivity in the flux-control system. The most sensitive factors for the predicted thickness of the mixing zone at this depth in the head-control system include transverse dispersivity, molecular diffusivity, porosity, and longitudinal dispersivity. However, the saturated hydraulic conductivity plays a more important role in predicting the thickness of the mixing zone at this depth in the flux-control system. The sensitivity analysis improves our understandings for the development of seawaterfreshwater mixing zone during seawater intrusion both in a head-control system and a flux-control system. The FFD simulation experiment is a useful method for uncertainty analysis of the predicted results, particularly when only limited information is available about the values of input factors. If statistical distributions of input factors are provided, formal stochastic methods such as the Monte Carlo analysis should be applied to compare with and test the results from the FFD simulation experiments. A hydraulically homogeneous and isotropic coastal aquifer was assumed in the present simulations and uncertainty analysis. This assumption neglects the spatial variability and anisotropy of the aquifer properties on the predicted mixing zones (Smith, 2004), and was made to allow the study to  

500

ZHAO Zhong-wei et al. / China Ocean Eng., 29(4), 2015, 489 – 502  

focus on other uncertainty factors including particularly the boundary conditions. The hydrodynamic boundary conditions also influence the formation of the mixing zone, such as seasonal variance of groundwater discharge, tidal and wave forces (Michael et al., 2005; Robinson et al., 2007; Li et al., 2008; Xin et al., 2010). Further studies should explore the combined effects of aquifer heterogeneity, seasonal variations of freshwater discharge, and tidal and wave oscillations. These studies are likely to be computationally intensive and should consider the application of parallel computing technologies. Appendix: Definition for Position and Thickness of the Mixing Zone The position (S) of a mixing zone at a depth is given by the horizontal distance from the coastal line to the position of the contour line with 50% salt concentration at the depth. A positive value of S indicates that the mixing zone intrudes landward from the coastal line, while a negative value indicates that the mixing zone is pushed seaward from the coastal line by the freshwater tube. The definition of the thickness (D) for a mixing zone follows that introduced by Lu et al. (2013). The details for the definition are shown in Fig. A: ABC and DEF are the contour lines with 90% and 10% salt concentrations, respectively; A and D, B and E, and C and F are at depth z1, z2, and z3, respectively. Following the definition of the interface-averaged thickness (Lu et al., 2013), the thickness for the mixing zone ABED is calculated as:

W1 

A1

, (A) DE where A1 is the area of ABED and A'B'ED is the parallelogram with the area equal to A1. Similarly, the thickness for the mixing zone BCFE is calculated as: W2 

A2

, (B) BC where A2 is the area of BCFE, and BCF'E' is the parallelogram with area equal to A2. The thickness of the mixing zone at depth z = z2 is defined as the average thickness of W1 and W2 with deviator Δz = z1 – z2 = z3 – z2, that is: W

1 (W1  W2 ) . 2

(C)

Fig. A. Characteristics for the thickness of a mixing zone.

References Abarca, E., Carrera, J., Sánchez-Vila, X. and Dentz, M., 2007. Anisotropic dispersive Henry problem, Adv. Water Resour., 30(4): 913–926.  

ZHAO Zhong-wei et al. / China Ocean Eng., 29(4), 2015, 489 – 502 

501

Bear, J. and Cheng, A. H.-D., 2010. Modeling Groundwater Flow and Contaminant Transport, Springer Dordrecht, Heidelberg, London, New York. Box, G. E., Hunter, W. G. and Hunter, J. S., 1978. Statistics for Experimenters: An Introduction to Design, Data Analysis, and Model Building, Wiley, New York. Cheng, X. J., Zhan, W., Guo, Z. R. and Yuan, L. R., 2012. A modeling study on saltwater intrusion to western four watercourses in the Pearl River Estuary, China Ocean Eng., 26(4): 575590. Dagan, G. and Zeitoun, D. G., 1998. Seawaterfreshwater interface in a stratified aquifer of random permeability distribution, J. Contam. Hydrol., 29(3): 185–203. Diersch, H. J. G. and Kolditz, O., 2002. Variable-density flow and transport in porous media: Approaches and challenges, Adv. Water Resour., 25(8): 899–944. Glover, R., 1959. The pattern of fresh-water flow in a coastal aquifer, J. Geophys. Res., 64(4): 457–459. Herckenrath, D., Langevin, C. D. and Doherty, J., 2011. Predictive uncertainty analysis of a saltwater intrusion model using null-space Monte Carlo, Water Resour. Res., 47(5): W05504. Kerrou, J., Renard, P., Lecca, G. and Tarhouni, J., 2010. Grid-enabled Monte Carlo analysis of the impacts of uncertain discharge rates on seawater intrusion in the Korba aquifer (Tunisia), Hydrological Sciences Journal-Journal Des Sciences Hydrologiques, 55(8): 1325–1336. Kuan, W. K., Jin, G., Xin, P., Robonson, C., Gibbes, B. and Li, L., 2012. Tidal influence on seawater intrusion in unconfined coastal aquifers, Water Resour. Res., 48(2): W02502. Lecca, G. and Cau, P., 2009. Using a Monte Carlo approach to evaluate seawater intrusion in the Oristano coastal aquifer: A case study from the AQUAGRID collaborative computing platform, Phys. Chem. Earth, 34(10-12): 654–661. Li, H. and Jiao, J., 2013. Quantifying tidal contribution to submarine groundwater discharges: A review, Chinese Sci. Bull., 58(25): 3053–3059. Li, H., Boufadel, M. C. and James, W. W., 2008. Tide-induced seawater-groundwater circulation in shallow beach aquifers, J. Hydrol., 352(1-2): 211–224. Lu, C., Chen, Y., Zhang, C. and Luo, J., 2013. Steady-state freshwater–seawater mixing zone in stratified coastal aquifers, J. Hydrol., 505, 24–34. Ma, Q., Li, H., Wang, X., Wang, C., Wan, L., Wang, X. and Jiang, X., 2015. Estimation of seawater–groundwater exchange rate: Case study in a tidal flat with a large-scale seepage face (Laizhou Bay, China), Hydrogeol. J., 23(2): 265–275. Michael, H. A., Mulligan, A. E. and Harvey, C. F., 2005. Seasonal oscillations in water exchange between aquifers and the coastal ocean, Nature, 436(7054): 1145–1148. Moore, W. S., 1996. Large groundwater inputs to coastal waters revealed by 226Ra enrichments, Nature, 380(6575): 612–614. Morris, M. D., 1991. Factorial sampling plans for preliminary computational experiments, Technometrics, 33(2): 161174. Naji, A., Cheng, A. B. D. and Ouazar, D., 1999. BEM solution of stochastic seawater intrusion problems, Eng. Anal. Bound. Elem., 23(7): 529–537. Paster, A. and Dagan, G., 2007. Mixing at the interface between two fluids in porous media: A boundary-layer solution, J. Fluid Mech., 584(1): 455–472. Paster, A., Dagan, G. and Guttman, J., 2006. The salt-water body in the northern part of Yarkon-Taninim aquifer: Field data analysis, conceptual model and prediction, J. Hydrol., 323(1-4): 154–167. Pool, M. and Carrera, J., 2011. A correction factor to account for mixing in GhybenHerzberg and critical pumping rate approximations of seawater intrusion in coastal aquifers, Water Resour. Res., 47(5): W05506.  

502

ZHAO Zhong-wei et al. / China Ocean Eng., 29(4), 2015, 489 – 502  

Rajabi, M. M. and Ataie-Ashtiani, B., 2014. Sampling efficiency in Monte Carlo based uncertainty propagation strategies: Application in seawater intrusion simulations, Adv. Water Resour., 67, 46–64. Robinson, C., Li, L. and Barry, D., 2007. Effect of tidal forcing on a subterranean estuary, Adv. Water Resour., 30(4): 851–865. Saltelli, A., Chan, K. and Scott, E. M., 2000. Sensitivity Analysis, Wiley, New York. Sanford, W. E. and Pope, J. P., 2010. Current challenges using models to forecast seawater intrusion: Lessons from the Eastern Shore of Virginia, USA, Hydrogeol. J., 18(1): 73–93. Sanz, E. and Voss, C. I., 2006. Inverse modeling for seawater intrusion in coastal aquifers: Insights about parameter sensitivities, variances, correlations and estimation procedures derived from the Henry problem, Adv. Water Res., 29(3): 439–457. Simmons, C. T., 2005. Variable density groundwater flow: From current challenges to future possibilities, Hydrogeol. J., 13(1): 116–119. Simmons, G. M., Jr, 1992. Importance of submarine groundwater discharge (SGWD) and seawater cycling to material flux across sediment/water interfaces in marine environments, Mar. Ecol.-Prog. Ser., 84, 173– 184. Volker, R. E. and Rushton, K. R., 1982. An assessment of the importance of some parameters for seawater intrusion in aquifers and a comparison of dispersive and sharp-interface modelling approaches, J. Hydrol., 56 (3-4): 239–250. Voss, C. I. and Provost, A. M., 2010. SUTRA, A Model for Saturated-Unsaturated, Variable-Density Ground-Water Flow with Solute or Energy Transport, Water-Resources Investigations Report, Reston, Virginia, U.S. Geological Survey. Wang, X., Li, H., Jiao, J. J., Barry, D. A., Li, L., Luo, X., Wang, C., Wan, L., Wang, X., Jiang, X., Ma, Q. and Qu, W., 2015. Submarine fresh groundwater discharge into Laizhou Bay comparable to the Yellow River flux, Sci. Rep., 5, Article No: 8814, doi:10.1038/srep08814. Werner, A. D. and Simmons, C. T., 2009. Impact of sea-level rise on sea water intrusion in coastal aquifers, Ground Water, 47(2): 197–204. Werner, A. D., Bakker, M., Post, V. E. A., Vandenbohede, A., Lu, C. H., Ataie-Ashtiani, B., Simmons, C. T. and Barry, D. A., 2013. Seawater intrusion processes, investigation and management: Recent advances and future challenges, Adv. Water Resour., 51, 3–26. Werner, A. D., Ward, J. D., Morgan, L. K., Simmons, C. T., Robinson, N. I. and Teubner, M. D., 2012. Vulnerability indicators of sea water intrusion, Ground Water, 50(1): 48–58. Xia, Y. Q. and Li, H. L., 2012. A combined field and modeling study of groundwater flow in a tidal marsh, Hydrol. Earth Syst. Sci., 16(3): 741–759. Xin, P., Jin, G., Li, L. and Barry, D. A., 2009. Effects of crab burrows on pore water flows in salt marshes, Adv. Water Resour., 32(3): 439–449. Xin, P., Robinson, C., Li, L., Barry, D. A. and Bakhtyar, R., 2010. Effects of wave forcing on a subterranean estuary, Water Resour. Res., 46(12): W12505. Xue, Y., Wu, J., Liu, P., Wang, J., Jiang, Q. and Shi, H., 1993. Sea-water intrusion in the coastal area of Laizhou Bay, China: 1. Distribution of sea-water intrusion and its hydrochemical characteristics, Ground Water, 31(4): 532–537. Zhao, Z. W., Zhao, J. and Fu, C. S., 2013. Uncertainty analysis of seawater intrusion forecasting, Water Sci. Eng., 6(4): 380–391. Zheng, J. H.,Yan, Y. X. and Zhu Y. L., 2002. Three dimensional baroclinic numerical model for simulating fresh and salt water mixing in the Yangtze Estuary, China Ocean Eng., 16(2): 227238.