Joint estimation of states and parameters of two-layer coastal aquifers based on ENKF

Groundwater resources are one of the most important sources of urban water supply, especially for the northern cities of China (Zheng et al. 2015; Chen et al. 2020). The exploitation and utilization of groundwater resources are increasing with the development of the economy in China, which leads to the continuous decline of groundwater levels and many other adverse environmental consequences, such as land subsidence (Chen et al. 2020), water quality deterioration (Purushotham et al. 2011), seawater intrusion (Ma et al. 2020), etc. In the past few decades, climate impacts, together with those of excessive human water use have changed the country's water availability structure (Khaki et al. 2018). At the same time, the spatial and temporal characterizations of groundwater were also altered (Chakraborty et al. 2020). Therefore, identifying the characterizations of groundwater has become one of the key points in scientifically evaluating and managing groundwater resources. Many methods can be used to predicting the dynamics of groundwater now, such as the evaluation of water balance (Portoghese et al. 2005), numerical simulation (finite-difference method (Ansarifar et al. 2020), finite-element method (Huang et al. 2020)), regression analysis (Chenini & Msaddek 2020), and spectral analysis (Kim et al. 2005). With the improvement of computing techniques, some calculation methods, such as the artificial neural network and wavelet transform (Rahman et al. 2020), have been applied to forecast the groundwater level gradually. At present, relatively traditional methods have been widely used in the analyses of groundwater dynamics. However, it was difficult to integrate real-time data into the model efficiently. The method of sequential data assimilation could utilize the monitored value in time, and update the parameters and variables based on historical estimations and the latest data (Tong et al. 2010). As one of the most popular methods of sequential data assimilation, the Ensemble Kalman Filter (ENKF) was first proposed in 1994 (Evensen 1994) and was refined by Burgers et al. (1998). ENKF extended the traditional linear Kalman Filtering (KF) (Zhou et al. 1991) to the non-linear field, and showed a more stabilized calculation than the Extended Kalman Filter (EKF) (Evensen 1992) in the non-linear field.

ENKF has previously been applied for large-scale nonlinear models in oceanography (Haugen & Evensen 2002) and hydrology (Margulis et al. 2002) successfully. The first application of ENKF for flow in porous media was petroleum engineering (Naevdal et al. 2005). Then, ENKF was introduced into hydrogeology gradually to acquire the hydraulic conductivity and the migration of pollutant plumes in the synthetic models (Huang et al. 2009) under different situations (Chen & Zhang 2006; Hendricks Franssen & Kinzelbach 2008), under various boundary conditions of the model (Tong et al. 2010). While some problems have arisen when operating the ENKF, such as filter inbreeding, unstable calculation processes, and the effect of observation data, some scholars (Hendricks Franssen & Kinzelbach 2008; Nan & Wu 2011; Schöniger et al. 2012; Cui & Wu 2013) have put forward the solutions for these deficiencies. Hendricks Franssen & Kinzelbach (2008) provided a solution for the filter inbreeding problem. Nan & Wu (2011) introduced distance-related weight to the ENKF which could stabilize the assimilation process using a small size ensemble. Schöniger et al. (2012) applied nonlinear, monotonic transformations to the observations rendering them Gaussian. Cui & Wu (2013) reported that the appropriate spatial and temporal distribution of observation wells with low density was superior to that of observation wells with high density. In recent years, ENKF was used to analyze pollutants in groundwater (Ross & Andersen 2018), identify pollutants in surface water (Wang et al. 2019), invert the parameters of unsaturated zones (Yu et al. 2019), analyze the concentration data of the density-dependent flow model in coastal areas (Yoon et al. 2020), validate the model of surface water and groundwater (Khaki et al. 2018), and recognize the relationship among the atmosphere, surface water and groundwater (Gelsinari et al. 2020).

However, the above-mentioned studies about ENKF were usually based on virtually synthetic models and rarely applied in the real world. The ENKF has not been integrated into the two-layer confined model in the coastal area so far. Thus, conducting a study on the complex groundwater model in the real world is needed (e.g. the groundwater model of the coastal area subjected to strong tidal effects). This study aims to evaluate the capability of ENKF in forecasting groundwater levels and inverting the hydraulic conductivity of the heterogeneous field in the relatively complex study area. The advantages and drawbacks of ENKF were also comparatively analyzed. In this study, we introduced ENKF into a petrochemical project in the coastal area of Tianjin city, took the two-layer confined flow model as a prediction function, assimilated the measured groundwater levels into the model, forecast groundwater levels, and inverted the hydraulic conductivity to obtain better prediction results. This paper was structured as follows. Firstly, the operation of ENKF, site description, calculations schemes, and evaluating indexes about the capability of ENKF were described. Then the results of forecasting groundwater levels and second inverse K by ENKF were presented. To show the efficiency of ENKF in forecasting aspects, comparative analyses were depicted as discussions. Finally, conclusions and guidelines for future work were given.

The ENKF is a widely used sequential data-assimilation method based on a Monte Carlo approximation that can update model parameters as well as model variables included in the state vector using various types of serially dynamic observations. ENKF has a better performance than the standard Kalman Filter (KF) as KF needs to compute and propagate the error covariance matrix explicitly in time, which may pose a significant computation burden for large or nonlinear problems. In contrast, ENKF abandoned the forecasting of the error covariance matrix in the model and directly carried out the multiple integrations based on ensembles via the Monte Carlo method, which means it can be calculated based on the continuously updated ensemble of realizations (Chen & Zhang 2006). Therefore, one of the reasons why ENKF widely used was that the number of calculations was greatly reduced. The main procedure for ENKF is as follows:

• (1)
  Generation of the initial ensemble

  

  (1)

  where X is an ensemble which concluded some samples; is a sample in the ensemble X; m is the size of samples;
• (2)
  Forecasting the parameter and state vector

  

  (2)

  where is the forecasting value of the sample k at time step j, j stands for the time step, k stands for the number of the sample; is the analyzed value of the sample k at time step j− 1, the superscripts f and a indicate the forecast and assimilation procedure, respectively; F is a forecast operator, representing the flow equations for our study.
• (3)
  Updating the parameter and state vector

  

  (3)

  

  (4)

  where is the analyzed value of the sample k at time step j; denotes the Kalman gain; is the observation vector (monitored data) at the time step j; is the observation operator which represents the relationship between the state vector and the observation vector; denotes the state error covariance matrix; is the error covariance matrix of the observations.
• (4)

  The assimilated value of ENKF

• Finally, averaged valuesof all samples were set as the assimilated value at time step j:

  

  (5)

• The four calculated steps mentioned above were the main procedure of ENKF at time step j, then the forecast model will run until new observations become available. so values at next time step j+ 1 could be forecast by returning the 2nd calculated step, namely Equation (2) and the forecasting values could fit the monitored data at the next time step better by assimilating the monitored data at the last time step into the ENKF model (Shen 2014). The above steps were realized by an interactive operation between MATLAB and GMS in this study.

The study area located in the intertidal zone which was reclaimed by the silty sand was about 15 km². According to the hydrogeological survey, the aquifer in the intertidal zone was composed of two confined layers, and the depth of groundwater level was relatively shallow due to the influence of overlying additional pressure and the surrounded seawater. Figure 1 shows the location and the specific environment of the study area. Q1, Q2, Q3, Q4, Q5, Q7, Q8, W4, and W5 were observation wells. Figure 1(b) presents the partitioned blocks of K in layer 2, which was divided into nine blocks. This was the first inversion of K. The first inversion was obtained by using the method of trial-and-error based on the two-period groundwater levels. That is to say, the first inversion K was obtained through changing the value of K and making the predicted groundwater levels fit to monitored levels. The process of changing the values and blocks of K of the groundwater model was the implementation of trial-and-error method. The K of nine blocks in layer 2 were named by observation wells, i.e., Q1-K2, Q2-K2, Q3-K2, Q4-K2, and so on. In each block, as shown in Figure 1(b), there was a notation like ‘W5-K2(2.8)’ which meant the first inverse value of the hydraulic conductivity was 2.8 m⁻¹ d through the method of trial-and-error in the W5 block. Since the hydraulic conductivity of layer 1 was not partitioned (not displayed in the figure), there were 10 inverse blocks of the hydraulic conductivity which would be inverted by ENKF. The first inverse value of partitioned K was utilized as the initial value in ENKF which means the 10 blocks values of K would be inverted by ENKF secondly.

The study area was discretized into squares as shown in Figure 2, the distance between two sides of an element was 10 m, and the observation points were refined, with a total of 57,956 elements and 123,717 nodes. The model was divided into two layers vertically, and the single-layer was divided into 28,978 elements. A set of initial samples conforming to Gaussian distribution was generated by superimposing the water level value on the 0th day with white noise (the mean value of 0 m and a standard deviation of 0.05 m (Hu 2018)), to assimilate and forecast groundwater levels of 9 observation points for 29 days in this model. Then the same white noise used in water level was superimposed (i.e., the mean value was 0 and the standard deviation was 0.05) with logarithm value of the initial hydraulic conductivity (lgK) which was the first step before the second inversion of K through ENKF. As the model was built from a real field case, the intervals between the monitored time were different, thus the assimilated time intervals of ENKF were also different. According to the monitored time in the local hydrogeological investigation, the assimilated time intervals were set 1.75 d (two), 1.79 d (six), 1.83 d (three), 1.84 d (three), 1.87 d (one), and 1.88 d (one), with a total of 16 intervals summarized 29 days.

The essence of ENKF was a calculation based on the statistics of samples. In theory, the larger samples stand for the higher efficiency of the calculation, but the relatively large samples in the ensemble would increase the computing cost. Some studies (Evensen 1994; Naevdal et al. 2005) have shown 100 samples used in ENKF could achieve good forecasting groundwater levels which fitted to the monitored data well, the accuracy of the forecasting results could satisfy the need for a general application in the real engineering project. For example, the result of forecasting groundwater levels by ENKF could be used for the groundwater environmental impact assessment for some petroleum projects and the early warning of the groundwater levels for some mining projects. Combined with previous studies and considering the area of this study, the 100 samples were adopted, that is, m = 100 (k = 1, 2, 3 … 100). On the other hand, considering the accuracy of water level probes and the location of study area is in the intertidal zone, the observation error in the groundwater level was set as 0.01 m and the groundwater level was monitored at each assimilated step during the operation of ENKF (Hu 2018). The observation error of K was set as 5% of their value and it was only observed in the first time step (i.e., 1.84 d). The relative observation error of 5% was adopted to study the soil moisture in some previous studies (Shi et al. 2015; Hu 2018). In this study, the relative observation error of 5% for K was employed to study whether the same observation error used in both layers could present the same convergent results in ENKF. In reality, it is difficult to define the observation error of K, the value of relative error (5%) was just to test the performance of ENKF in the two-layer confined aquifers.

RMSE was used to evaluate the overall capability of ENKF, the smaller RMSE often means the better performance of ENKF. SPREAD was used to evaluate the convergence of ENKF, the rapidly decreased values of SPREAD indicated the fast convergence of ENKF. When E was more than 0 demonstrated the forecasting value by ENKF fitted the monitored data better than the initial simulation, vice versa.

In addition, comparative analyses among ENKF, the initial model, and the modified model were done to see whether ENKF, to some extent, could improve the accuracy of forecasting groundwater levels. In terms of the comparative analyses between the forecasting values by ENKF and the initial simulation, the groundwater levels at some observation points were chosen to compare directly in the following sections and the second inverse K by ENKF was introduced to the initial model to get the modified model. The forecasting error between the modified model and the initial model was analyzed. The analytical method was a calculation of relative errors that subtracted forecasting errors of the modified model from the initial model. If the relative errors were calculated above 0, indicated that there was a certain improvement in the accuracy of forecasting water level, otherwise, no obvious improvement existed in terms of the modified model. The logic diagram of this study can be seen in Figure 3.

Results of forecasting groundwater levels through ENKF presented a good fit for the real-field monitored values as shown in Figure 4. SPREAD dropped rapidly in an exponential form, demonstrated that the assimilation of ENKF converged quickly. RMSE was small in the assimilated process, it decreased to 0.16471 at the final time step (i.e., the 16th time step), and the relatively small value of RMSE indicated that forecasting values were generally consistent with the real-field monitored values. However, RMSE showed a fluctuated decline process, which suggested that ENKF did not perform well at some assimilated time steps (mainly the odd time steps of assimilation). This may be related to the location of the study area, ENKF operated well at some time steps (even time steps) because of the well-fluctuated groundwater levels but performed worse at some time steps due to the relatively weak fluctuations in monitored values. This phenomenon showed consistency with the previous studies (Cui & Wu 2013; Song et al. 2014). The E (efficiency of calculation in ENKF) at the final assimilated time step was calculated using Equation (11), according to forecasting values of ENKF and the initial simulation which did not assimilate the real-field monitored values. The result of E was 68.55%, which showed ENKF performed better than the initial simulation in fitting the monitored data to a certain extent.

The absolute errors at all observation points were calculated by subtracting real-field monitored values from forecasting values, as shown in Figure 5. Among groundwater levels of nine observation points forecast by ENKF, seven observation points had relatively good results with the absolute errors were less than 0.1 m. The absolute error of Q8 was close to 0 m, presenting the best result of seven observations, and ENKF just presented a relatively good forecasting result at Q1, from Figure 5 the absolute error of Q1 fluctuated in the process of ENKF. However, the forecasting results of Q5 and W5 showed a poor performance, which could be explained by the local hydrogeological conditions and the intertidal zone. Q5, located in the middle of the study area, was less affected by tides and the hydraulic conductivity was small, these two factors led to the weak fluctuations in the groundwater level. Since the groundwater level at Q5 provided less information, ENKF could not forecast the slight fluctuations in groundwater level. Although the absolute error of the groundwater level at Q5 was large (in Figure 5), the error was decreased gradually with the operation of assimilation, which indicated that ENKF needed more long time monitored data at Q5 to forecast water level. The groundwater level at W5 fluctuated strongly because it was located at the tidal boundary (Figure 1). Although the fluctuations in groundwater levels could make ENKF perform well (Song et al. 2014), the inevitable phenomenon was that the monitored data associated with much more noise from the real field could make ENKF perform variously between different assimilated steps. Thus ENKF could perform well when monitored data of real field were subjected to a range of interferences, while for the monitored data with relatively strong or weak interferences i.e., the intensity of interferences was out of the range accepted by ENKF, and ENKF would produce a bad result.

The 10 blocks values of hydraulic conductivity (K) in the two-layer confined aquifer were inverted by ENKF secondly of which the initial values were obtained by the first inversion based on the two-period water levels. The second inverse values of K at the final assimilated step, the initial values of K, and the improvement degrees are shown in Table 1. The improvement degrees were calculated by the differences between the second inverse values and the initial values dividing the initial values. It can be seen from Table 1 that the second inverse K of layer 2 smoothly converged to a specific value within 16 assimilation steps, with little difference from initial values. That means when the relative observation error of K was 5%, the assimilation process was convergent for layer 2, while the second inverse K of layer 1 was quite different from the initial value, and even went beyond the initial scope of interferences. To further explain the differences of inverse K between the two layers in the same model when taking the same observation error, the convergent details of 100 samples in the assimilation process are drawn in Figure 6. To avoid repetition, the convergence of Q1-K2 (hydraulic conductivity of Q1 block in layer 2) was selected to demonstrate the capability of ENKF in layer 2. Figure 6(a) indicates that Q1-K2 had a great convergence at time step 2, the values of samples kept approaching the initial values with the process of assimilation. There was no obvious change in the mean value of samples after the 8th assimilation step. While the divergence appeared in the assimilated process of ENKF for K1, values of samples could not converge as shown in Figure 6(b). The reason why the divergence appeared was that the noise of model was too large to make the model inaccurate or the observation error was too big to make the Kalman gain matrix small, which made the innovation account for a large proportion and lead to the inability to fit the real-field monitored values (Evensen 2009; Dan et al. 2016). Because of the same numerical simulation adopted in this paper, the system noise was considered as zero, so we tried to make the assimilated process of K1 converge by reducing the observation error and increasing the proportion of the original value. When relative observation error 4% was adopted for the layer 1 in ENKF, although the mean value of K1 did not converge to a specific value within 16 assimilation steps, it was obvious that the convergence of samples appeared as seen in Figure 6(c). While for the relative observation error 3%, ENKF performed better and converged to a specific value in a short time as demonstrated by Figure 6(d). The observation error had a great influence on the assimilation process as concluded from the above. K2 (the nine blocks of K in layer 2) could converge with a larger observation error in the ENKF system, while for K1 (the K in layer 1), a smaller observation error should be applied to make the assimilation converge. It is interesting that ENKF shows a huge difference for the same parameter in the same model. The reason may because the model is two-layer and the location of study area is an intertidal zone. Anthropogenic activities and tidal influences are two main factors that can strongly affect the inversion of K, especially for layer 1 in this study. As described above, the K of layer 1 was not partitioned and the initial inversion based on the two-period groundwater levels could have a larger error compared with the K in layer 2. While for the operation of ENKF, the same observation error of 5% was adopted in the two layers, which led to the divergence of ENKF for K1.

The above analysis showed that even if the same error was adopted in the same model for the same parameter, the performance of ENKF could be different. Especially for the real monitored data subjected to much more external interferences, the inversion of K had a larger error based on the monitored data. If the inverse K with a large error and the inverse K with a small error were inverted by ENKF simultaneously, the possibility of filter divergence was greater. Thus, observation errors of the same parameter should be considered different in the real-field model, especially for the two or more layers of groundwater models.

The forecasting groundwater levels by ENKF that fitted the real monitored levels have been analyzed in the above sections. To further demonstrate the positive and negative effects of ENKF, forecasting levels by ENKF and the initial simulation without any modification were compared as displayed in the following sections. The results showed that the forecasting by ENKF performed worse than the initial simulation at W5 and Q5, while performances of ENKF at the other seven observation points were better. The comparative analysis of the two methods is displayed in Figure 7, there were analyses of W5, Q5, Q4, and Q8, other observation points (Q1, Q2, Q3, Q7, and W4) were not reported here as they were similar to Q4 and Q8. As detailed above, since the intensities of the interferences at W5, Q5 were out of the range which could be accepted by ENKF, the efficiency of forecasting by ENKF was lower than the initial simulation. By the initial simulation, W5 showed a well-fluctuated curve which fitted the monitored data well due to its special locations, while for the other observation points, whether using ENKF or the initial simulation, all showed weak fluctuations as indicated by Figure 7(b)–7(d). In terms of the differences between values of forecasting and real monitored, ENKF performed better at the other seven observation points, which can be seen from Figure 7(c) and 7(d), the forecasting groundwater levels by ENKF at Q4 and Q8 were closer to the real monitored curves than the initial simulation. Both ENKF and the initial simulation have advantages and disadvantages in forecasting the groundwater level, but they could be complementary to improve the accuracy of the model. Thus, if the accuracy of forecasting groundwater levels was required to be higher, the combination of two methods could be considered to increase the forecasting accuracy of the flow field.

As mentioned before, the modified model was obtained by introducing the second inverse K presented in Table 1 into the initial simulation. Relative errors of nine observation points are displayed in Figure 8. Figure 8(a) shows that the relative errors of Q3, Q4, Q5, W4, and W5 were more than 0, which meant the modified model improved the accuracy in forecasting groundwater levels. Although the relative error of W5 showed a fluctuation revealing there were some unstable forecasting processes at some time steps in the modified simulation, this was likely to be accounted for by the special position of W5. From the values of relative errors at W5, the modified model improved the accuracy of forecasting distinctly compared with other observation wells. Q5 was located in the middle of the study area, the weak fluctuations made ENKF perform poorly, while the improvement degree of hydraulic conductivity at Q5 block (Q5-K2) was the highest among all inverse blocks from the Table 1. The relative large improvement degree in the hydraulic conductivity made Q5 present a large relative error that meant that the modified model was greater than the initial model. This means forecasting accuracy of groundwater levels in the middle of the study area could be enhanced by the large improvement degree of K through ENKF. Figure 8(a) shows the higher the improvement degree of hydraulic conductivity (mainly for the hydraulic conductivity of layer 2, K2), the greater the accuracy of groundwater level forecasting by the modified model. The relative errors calculated at Q2, Q7, and Q8 showed fluctuations as shown in Figure 8(b) which meant that the modified model forecasted the groundwater level unstably at the three observation points. For Q1, the modified model showed a poor performance within the 16-time steps. Although the performance of the modified model was not good at the four observation points, the trend that the higher improvement degree of K2 could enhance the accuracy of groundwater forecasting still appears in Figure 8(b). Figure 8 demonstrates that not every observation point displayed the good result when forecast by the modified model, which may be related to the fluctuation of groundwater level and the improvement degree of second inverse hydraulic conductivity. From Table 1, the improvement degrees of the second inverse hydraulic conductivity were small except for Q5-K2. Even if the improvement degrees of W5-K2, W4-K2, Q4-K2, and Q3-K2 were small, the modified model could also perform well in forecasting water level due to the strong fluctuations existed in the real-field monitored data. While Q1, Q2, Q7, and Q8 were located in the relatively middle of the model far away from the tidal boundary, monitored groundwater levels at these observation points showed a weaker fluctuation compared with W5, W4, Q4, and Q3. Thus, under the conditions of a small improvement degree of K2, the modified model showed fluctuating results at Q1, Q2, Q7, and Q8.

From the comparative analysis between the modified simulation and the initial simulation, if the observation points had strongly fluctuating groundwater levels, even a small improvement degree of K in the modified model could perform better than the initial model. For the weakly fluctuating groundwater levels, a large improvement degree of K needed to enhance the efficiency of forecasting by the modified model. This means that the modified model with a large improvement degree of hydraulic conductivity can compensate for the deficiency caused by weak fluctuation in groundwater levels and then improve the accuracy of forecasting waters.

In this study, an assimilated system based on ENKF and 2D confined flow models of two layers was constructed. Groundwater levels in the study area were forecast gradually with different time intervals by assimilating the real-field monitored data and 10 blocks of hydraulic conductivity were inverted by ENKF secondly. Comparative analyses among ENKF, the modified model according to the second inverse K, and the initial model were discussed in this study. The key findings were as follows:

• (1)

  ENKF could utilize all the real-field monitored data efficiently, the performance was very acceptable when assimilating and forecasting groundwater levels. However, when monitored data at some observation points were beyond the intensity of interferences accepted by ENKF, the performance was not good.

• (2)

  Results of the second inverse hydraulic conductivity demonstrated that even if the same parameter was used in the same model, the relative observation errors of the same parameter could be different, especially for the two or more layer groundwater model. To improve the accuracy of forecasting, different relative observation errors of the same parameter should be considered in ENKF for the complex groundwater model.

• (3)

  The combination of ENKF and the initial model could be a considered a method to improve the accuracy of forecasting. For example, ENKF can be used at almost all observation points, and the initial model could be a supplementary method to improve the accuracy of some observation points at which ENKF performed poorly.

• (4)

  The large improvement degree of hydraulic conductivity could enhance the efficiency of forecasting especially for the observation points located in the middle of the model which could compensate for the deficiency caused by weak fluctuation in groundwater levels

ENKF was successfully introduced into the real-field confined flow model of two layers which was located in the coastal area of Tianjin city and performed well on the whole. This study demonstrated how to improve the accuracy of forecasting groundwater levels by combining ENKF and initial model or using a modified model with the second inverse K by ENKF. These methods could be used to improve the accuracy of forecasting groundwater levels, which could provide the relevant reference for urban management of groundwater resources and give early warning of water levels, especially for the coastal areas subjected to seawater intrusion. While the range of interferences accepted by ENKF was not analyzed quantitatively in this study, this aspect represented the direction for further studying the forecasting by ENKF.

The authors sincerely thank Dr Yang at the Laboratory of Eco-hydrology and Water-soil Conservation for her help. The authors acknowledge the valuable comments from two anonymous reviewers and editors.

No potential conflict of interest was reported by the authors.

All relevant data are included in the paper or its Supplementary Information.