LIU Jin-bing,JIANG Si-min,*,ZHOU Nian-qing,CAI Yi,CHENG Lu,WANG Zhi-yuan

1 Department of Hydraulic Engineering,College of Civil Engineering,Tongji University,Shanghai 200092,China.

2 State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering,Nanjing Hydraulic Research Institute,Nanjing 210029,China.

Abstract:When groundwater pollution occurs,to come up with an efficient remediation plan,it is particularly important to collect information of contaminant source (location and source strength) and hydraulic conductivity field of the site accurately and quickly.However,the information can not be obtained by direct observation,and can only be derived from limited measurement data.Data assimilation of observations such as head and concentration is often used to estimate parameters of contaminant source.As for hydraulic conductivity field,especially for complex non-Gaussian field,it can be directly estimated by geostatistics method based on limited hard data,while the accuracy is often not high.Better estimation of hydraulic conductivity can be achieved by solving inverse groundwater problem.Therefore,in this study,the multi-point geostatistics method Quick Sampling (QS) is proposed and introduced for the first time and combined with the iterative local updating ensemble smoother(ILUES) to develop a new data assimilation framework QS-ILUES.It helps to solve the contaminant source parameters and non-Gaussian hydraulic conductivity field simultaneously by assimilating hydraulic head and pollutant concentration data.While the pilot points are utilized to reduce the dimension of hydraulic conductivity field,the influence of pilot points’layout and the ensemble size of ILUES algorithm on the inverse simulation results are further explored.

Keywords: Inverse groundwater problem;Data assimilation;Multi-point Geostatistics;Quick Sampling;Non-Gaussian hydraulic conductivity field

Introduction

Accurate and timely investigation of pollution source parameters is very important for formulating a reasonable remediation plan to deal with the pollution occurring in aquifers.The data assimilation method is often used to estimate source parameters through inverse modelling.

The purpose of data assimilation is to reduce the deviation between the predicted and observed values of the simulation model,thus to improve its estimation accuracy and prediction ability.Many researches have achieved good results in this field(Guneshworet al.2018;Jha and Datta,2013;JIANG Si-minet al.2013;XIA Xue-minet al.2019).

Another key factor affecting the migration of contaminants is the permeability of the aquifer,thus accurate characterization of the hydraulic conductivity field also helps to develop a reasonable and efficient pollution remediation plan.However,the actual hydraulic conductivity field cannot be obtained by direct measurements,and can only be estimated from several sampling points in the aquifer.Geostatistical methods are widely used in the estimation of formation thickness,buried depth,oil saturation,and other parameters in the petroleum exploration field (LI Li and WANG Yong-gang,2006).In recent years,geostatistical methods have also been used to quantitatively describe hydrogeological parameters such as hydraulic conductivity (LIU Wen-tinget al.2010;LIU Ling-linget al.2009).Traditional Two-point Geostatistics (TPG) can only be used to estimate geological variables with simple structures,such as the hydraulic conductivity field following Gaussian distribution.However,it is difficult to effectively describe complex geological structures,such as non-Gaussian hydraulic conductivity fields,river channels,etc.(YANG Pei-jie,2014;LUO Hong-meiet al.2015).Therefore,Multi-point Geostatistics (MPG) method was developed to characterize those complex geological structures.Strebelle (2002) proposed the SNESIM algorithm to simulate nonlinear geological structure,but it has the disadvantage of high computational burden (Rezaeeet al.2013).Straubhaaret al.(2011) improved SNESIM and further proposed the IMPALA algorithm that reduced the memory requirement to achieve higher computing efficiency.Whereas,all of these MPG methods can only simulate discontinuous variables,which are not suitable for the hydraulic conductivity field and other continuous variables.Mariethozet al.(2010) proposed the DS (Direct Sampling)algorithm,a pixel-based MPG method,which directly searches for the optimal estimation value of simulated points in the training images.Gravey and Mariethoz (2020) developed the QS algorithm(Quick Sampling) to find the best matching point quickly by calculating the mismatch map between estimating point and entire training image.It is found that QS is superior to the DS in computational efficiency.Also,the QS algorithm is much easier to implement because it has fewer parameters than DS (Gravey and Mariethoz,2020).

Hydraulic conductivity field directly affects the groundwater flow and pollutant transport.Nonetheless,based on the limited and unevenly distributed hard data,geostatistical methods are often inaccurate when estimating hydraulic conductivity fields,especially for complex non-Gaussian fields.Therefore,the data assimilation method can also be employed to derive the hydraulic conductivity field from observation data using inverse simulation from which the pollution source parameters and hydraulic conductivity field can be simultaneously solved.

Stochastic methods,such as the Markov chain Monte Carlo method (MCMC),Ensemble Kalman filter (EnKF),Ensemble Smoother (ES) and its improved algorithms,are commonly used to estimate groundwater parameters by inverse simulation.Compared with other algorithms,En-semble Smoother(Chen and Oliver,2012;Emerick and Reynolds,2013;Evensen and Van Leeuwen,2000),as a batchprocessing algorithm,does not require frequent operations of model files and shows its advantage in low computational cost,so its popularity has largely increased recently.Considering the application of ES and its improved algorithm to inverse modelling of groundwater flow,the current research can be divided into the following categories:(1) only assimilating the observed piezometric head data to estimate the Gaussian hydraulic conductivity field(JU Leiet al.2018) or non-Gaussian hydraulic conductivity fields (CAO Zhen-danet al.2018;LI Liang-pinget al.2018a;LI Liang-pinget al.2018b);(2) assimilating both observation data of piezometric head and contaminant concentration,to derive pollu-tion source parameters and gaussian hydraulic conductivity field simultaneously(ZHANG Jiang-Jianget al.2018;YANG Ai-linet al.2020;MO Shao-xinget al.2019);(3) only estimating hydraulic conductivity field by using both observation data of piezometric head and contaminant concentration (MO Shao-xinget al.2020;ZONG Cheng-yuanet al.2020).

So far,few studies have achieved simultaneous estimation of contaminant source parameters and non-Gaussian hydraulic conductivity field by assimilating observation data of piezometric head and contaminant concentration.For the first time,this paper proposed and introduced QS algorithm of MPG to characterize the non-gaussian hydraulic conductivity field,and further build a new data assimilation framework (QS-ILUES) based on the ILUES (an improvement of IES).ILUES is more efficient and suitable for the inverse estimation of the non-gaussian hydraulic conductivity field(ZHANG Jiang-jianget al.2018).Combining the advantages of ILUES and QS,the QS-ILUES framework achieves the simultaneous inverse estimation of contaminant source parameters and non-gaussian hydraulic conductivity field through assimilating both observation data of piezometric head and pollutant concentration.In this paper,referring to the research of CAO Zhen-danet al.(2018),the pilot point method is used to reduce the dimension of the entire hydraulic conductivity field for high-dimensional parameter inverse estimation(Ramaraoet al.1995),which not only ensures the accuracy of the simulation but also improves the computational efficiency.

1 ldentification model of groundwater contaminant

Identification of groundwater pollution source is an inverse groundwater problem,including groundwater pollution transport model and the inverse model in which direct method,optimization method,data assimilation method,etc.can be used.The objective of this research is to obtain the information of groundwater pollution sources such as location,concentration,discharge time,etc.and other aquifer parameters mainly including hydraulic conductivity,storage coefficient,dispersivity,effective porosity and so on.These data are then assimilated to minimize the deviation between measured and simulated values of observation.

1.1 Groundwater contaminant transport model

Common pollutant transport programs include MOC3D,MT3DMS,RT3D,FEMWATER,FEFLOW and so on,of which MT3DMS is the most popular and also used in this study.It should be used with MODFLOW for both the groundwater flow and transport simulation.

The subsurface flow equation of MODFLOW is as follows:

Where:it is assumed that the main direction of the hydraulic conductivity is consistent with the direction of the coordinate axis,andKx,Ky,Kz[LT-1]are the components of the hydraulic conductivity in the axis ofx,yandz.h[L]is piezometric head;W[T-1]represents groundwater source and sink;Ss[L-1]is the storage coefficient of the aquifer media;t[t]is time.

The groundwater solute transport equation of MT3DMS is as follows:

Where:θ[dimensionless]is the porosity of aquifer;ck[ML-3]is the concentration of solutek;Dij[L2T-1]is the tensor of hydrodynamic dispersion coefficient;vi[LT-1]is the actual flow velocity in the aquifer,and the mathematical relationship with Darcy flow velocityqiisvi=qi/θ;qs[T-1]is the flow rate of the term of source and sink per unit volume for the aquifer;cksis the concentration of solutekin the term of source and sink;∑Rn[ML-3T-1]is the sum of chemical reaction term.

1.2 QS-ILUES framework

As the latest multi-point geostatistics method,QS helps to reconstruct the complex geological patterns such as non-Gaussian parameter fields(Gravey and Mariethoz,2020).It calculates the mismatch degree between each simulation point and the whole training image to obtain estimated value that sampling from the best matching candidate points.Local updating strategy of sample set and simple iteration process are adopted to improve ES algorithm to solve non-Gaussian fields(ZHANG Jiang-Jianget al.2018).In this paper,a new parameter estimation method is proposed to achieve simultaneous inverse estimation of contaminant source parameters and hydraulic conductivity field by constructing data assimilation framework between QS and ILUES.Considering the high dimension and nonlinear characteristics of non-Gaussian hydraulic conductivity fields,the QS-ILUES framework applies pilot points to reduce the dimension of hydraulic conductivity fields.The estimated hydraulic conductivity at pilot points together with the limited available ones are combined and used as the conditioning data for QS to estimate the whole hydraulic conductivity field.The data-assimilation framework is shown in Fig.1.

(1) Based on the training images and a few available hydraulic conductivity values in the aquifer,QS is used to generate the initial sample field of hydraulic conductivity.Hydraulic conductivity values from pilot pointsPeand sampled pollution source parametersSeare put together to constitute the initial parameter matrix for the inverse modelling:

(2)miis put into the groundwater modelF(·)Equation (3),and the observation values of simulated piezometric head and contaminant concentration can be obtained.

Where:mis the model parameters,andeis the observation error.

(3) The difference between observation valuesdrealand simulated valuesdsimof piezometric head and contaminant concentration is calculated,and the parameter matrix is updated based on the following formula:

Where:mais the updated parameter matrix,andKgis the Kalman gain.

(4) The hydraulic conductivity at the pilot points updated with ILUES is also treated as conditioning data together with a few available hydraulic conductivity values for QS to generate the new hydraulic conductivity field,and an update process for the contaminant source parameters and hydraulic conductivity field is conducted.Repeat step (2)-(4) until the pre-set iteration times are completed.

Fig.1 The data assimilation framework of QS-ILUES

2 Case study

2.1 Problem overview

As shown in Fig.2a,a two-dimensional confined aquifer with the size of 150 m × 90 m is heterogeneous and anisotropic,and divided into 30 rows and 50 columns with a cell size of 3 m as the initial model grid size.The grid size is adjusted automatically in the process of inverse simulation.The contaminant source is set to be located at the center of the grid.The groundwater flow is assumed to be steady,and the aquifer thickness is 5 m.Besides,the upper and lower boundaries are assumed to be Neumann boundary,while the left and right boundaries are assumed to be Dirichlet boundary with specified head of 15 m along the left and 10 m along the right.There is no external source and sink of groundwater flow terms.No contaminant is introduced at the beginning of the simulation.The porosity of the aquifer is 0.30,the longitudinal dispersivity is 5.0 m,and the horizontal dispersivity is 1.0 m.

The aquifer contains two facies and the general distribution of hydraulic conductivity field shows clear non-Gaussian characteristics.Both of the two facies follow the lognormal distribution,with the mean value of (lnK) of 5.0 and 1.0,respectively,and the same variance (σ2lnK) of 0.50.In addition,the correlation length of the two fields is the same,which is 10 m along bothxandydirections.The variogram types of the two fields are both exponential.The true hydraulic conductivity field is shown in Fig.2a.

Point source contaminant occurred in the aquifer,and the release location of contaminant is shown in Fig.2a.The possible location and area of the contaminant source (marked by the rectangle area with red dotted line in Fig.2a are determined through site investigation and they are also treated as the prior information of the contaminant source parameters (Table 1) for inverse simulation.The piezometric head and contaminant concentration at the observation points are simulated in MODFLOW and MT3DMS program.10 stress periods are defined and each has a duration of 5 days.The contaminant only releases in the first eight stress periods (Table 1).20 observation wells are placed in the aquifer (Fig.2a),and observation values of the piezometric head and contaminant concentration are read every two days from the 26thday to the 50thday.The errors of all observed values obey normal distributionN(0,σ2),andσobeys uniform distributionU(1,2).

Fig.2 (a) the reference hydraulic conductivity field;(b) training image of facies;(c) training image of hydraulic conductivity field

Table 1 Reference values and prior ranges of contaminant source parameters

2.2 Results and discussion

Large size training images (TI) are usually selected for pattern searching when the multipoint geostatistics is used for the characterization of non-Gaussian hydraulic conductivity of aquifers(LI Liang-pinget al.2018;CAO Zhen-danet al.2018;ZONG Cheng-yuanet al.2020).In practical applications,large size training images may not be satisfied.In this study,the available data include:Training images of facies obtained from the previous site investigation (Fig.2b),hard data or field measurements including 20 observation wells and the measured hydraulic conductivity (Fig.2a).Table 2 shows the estimated values of the two hydraulic conductivity fields obtained by SGEMS program on the basis of hard data,from which two hydraulic conductivity fields are generated and filled into the training image of facies to acquire the training image of the hydraulic conductivity field (Fig.2c).The values in Table 2 show that the estimation is more accurate for the highpermeability facies.For the low-permeability facies,significant errors commonly exist in all statistics except the mean value.In addition,the training image of facies in Fig.2b shows large deviation compared with the reference aquifer(Fig.2a),which means that the training image of hydraulic conductivity contains the stacked errors and is very different from the reference hydraulic conductivity field.

When QS is used to characterize the hydraulic conductivity field,154 pilot points,as shown in Fig.3,are also set up in the aquifer in addition to the hard data from the 20 observation wells.The initial ensemble of 10 contaminant source parameters in this study was sampled from their prior distribution (Table 1).After several iterations,the mean values of the ensembles are regarded as the estimated value of the contaminant source parameters.The estimation of hydraulic conductivity values follows the same procedure.

RMSE (Root mean Square Error) and AES(Average Ensemble Spread) are used to evaluate the quality of the simulation results.Where the lower the RMSE value is,the higher the parameter estimation accuracy is.AES reflects the dispersive degree of the parameter ensemble relative to the real value.The smaller the value is,the more the ensemble is concentrated near the real value,and the better effect the inverse simulation has.The calculation formulas of the two indexes are as follows:

In the calculation formula of RMSE,VestandVactrepresent the estimated value and true value of the ith inverse parameter respectively,andNis the total number of parameters.In the AES calculation formula,(Vest)ijrepresents the estimated value in thej-th ensemble of thei-th parameter,and other parameters have the same meaning as RMSE.

Table 3 The simulation setting of the case study

Fig.3 The distribution of pilot points

2.2.1 Influence of the placement mode of pilot point

For Case 1 and 2,the size of ensemble and iteration times of the ILUES algorithms are set the same,the difference is that the 154 pilot points are placed in different ways (Table 3),and the specific locations are shown in Fig.3.

In Table 5,the values of RMSE and AES show that the estimation accuracy of Case 1 for the parameters of contaminant source is higher than that of Case 2,while the accuracy of hydraulic conductivity at pilot points and hydraulic conductivity field for Case 2 is better than that of Case 1.However,by analyzing the specific values of source parameters for Case 2 in Table 4,it is found that except for the estimation accuracy of contaminant release concentration in the first stress period is slightly worse where the difference between estimated and true value is 7.42,and the relative deviation is 12.4%,the estimated values are close to the true values in other stress periods,indicating that Case 2 also shows a reliable inverse estimation.In addition,in terms of the estimation and variance of the hydraulic conductivity field,the accuracy of Case 2 is significantly better than that of Case 1 (Fig.4b and Fig.4c).

From the above analysis,the random distribution of pilot points may help to improve the estimation accuracy of contaminant source parameters,but it is not as good as the even distribution of pilot points for the estimation accuracy of hydraulic conductivity.Therefore,the pilot points are placed in the mode of even distribution in the subsequent cases.

Table 4 Estimation values of contaminant source parameters

Table 5 Inverse estimation accuracy of source parameters,hydraulic conductivity at pilot point,hydraulic conductivity field in Case 1～3

Fig.4 a.The reference hydraulic conductivity field;b-d.The estimated hydraulic conductivity field of Case 1～3;e.The training image of hydraulic conductivity field;f-h.The variance field for the estimated hydraulic conductivity field of Case 1～3.

2.2.2 Influence of the ensemble size of ILUES

Both the contaminant source parameters and hydraulic conductivity field have been accurately estimated in Case 2,and the high permeability zone of the hydraulic conductivity field can be roughly delineated from the result.However,there is still a gap between the estimated and reference field in the high permeability area circled by the red dotted rectangle in Fig.4c.For ILUES algorithm,increasing the size of ensemble can improve the estimation accuracy of parameters to some extent(ZHANG Jiang-Jianget al.2018).Therefore,Case 3 is constructed,where the ensemble size is increased to 1 500,and the pilot points are placed as the mode of even distribution.The inverse estimation accuracy of the hydraulic conductivity at the pilot points and the hydraulic conductivity field are shown in Table 5.Among the RMSE and AES of the three cases,the value of Case 3 is the minimum,which indicates that increasing the size of ILUES ensemble can improve the estimation accuracy of hydraulic conductivity at the pilot points and also hydraulic conductivity field to some extent.As a result,in the estimation field of Case 3 (Fig. 4d),the high permeability channel in the region indicated by the red dotted rectangle is closer to the reference hydraulic conductivity field.

In terms of the RMSE and AES,the accuracy of Case 3 is the lowest for the contaminant source parameters.However,compared with previous studies (ZHANG Jiang-jianget al.2018;YANG Ai-linet al.2020;MO Shao-xinget al.2019),the numerical cases in this paper have a larger prior estimation range for given contaminant source parameters,which increases the difficulty of parameter inverse estimation.According to the specific parameter values in Table 4,except for the large estimation error of contaminant release concentration in the second and third stress periods(the difference between estimated and true value are 8.58 and 7.36,respectively,and the relative deviation are 14.8% and 13.4%),the estimated values of other parameters are quietly close to the true values.Therefore,the inverse estimation accuracy of contaminant source parameters is completely acceptable.

3 Conclusions

(1) In this study,the multi-point geostatistics method Quick Sampling (QS) is proposed and introduced and a new data assimilation framework QS-ILUES is then developed by combination with the Iterative Local Update Ensemble Smoother(ILUES) algorithm to achieve simultaneous inverse estimation of contaminant source parameters and non-Gaussian hydraulic conductivity fields.

(2) When QS-ILUES framework is applied,it can help to identify the source location and contaminant release concentration in each stress period from a larger prior range,and characterize the high permeability channel accurately.

(3) When applying the descending dimension method of the pilot points to estimate the hydraulic conductivity field,the even distribution of pilot points shows a higher accuracy than the random distribution.Further,increasing the ensemble size of ILUES algorithm with evenly distributed pilot points can also improve the inverse estimation accuracy of hydraulic conductivity field to some extent.

(4) Due to the difficulty in obtaining training images,most studies on non-Gaussian hydraulic conductivity field focus on verifying the applicability of various methods by increasing the number of case tests.In the subsequent studies of this paper,the QS-ILUES method framework will be considered to apply on large-scale site.

Ackonwledgements

This work was supported by the Belt and Road Special Foundation of the State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering (No.2019nkzd01) and National Natural Science Foundation of China (42077176).