Groundwater flow and transport model of the Estonian Artesian Basin and its hydrological developments

Owing to such a complex situation, an integrative investigation and assessment of all Estonian water resources is needed for their optimum management. For that purpose, it is necessary to carry out a basin-wide three-dimensional and time-dependent groundwater modelling, coupled with the other components of the hydrologic cycle as well as with contaminant mass balance (NRC 2012; Singh 2012; Mulligan et al. 2014; Cai et al. 2015).

Groundwater flow and transport modelling are both basic directions of contemporary quantitative investigation of the subsurface water environment (Contoux et al. 2013; Zhu et al. 2013; Dell'Arciprete et al. 2014; Joyce et al. 2014; Turner et al. 2015). Furthermore, a number of specific researches are being developed towards sustainable groundwater use (Gleeson et al. 2012; Brown et al. 2015), inverse estimation of hydrogeological parameters (Doherty & Simmons 2013; Chiu 2014), application of tracers (Girmay et al. 2015; Mengistu et al. 2015), impacts of global change (Gorelick & Zheng 2015), optimisation of groundwater management (Yeh 2015), interaction between surface water and groundwater (Rinaldi et al. 2014; Sutton et al. 2015) and uncertainties connected with modelling (Refsgaard et al. 2010; Gupta & Nearing 2014; Kitanidis 2015).

A steady-state groundwater flow model of the entire Baltic Artesian Basin has been completed by Virbulis et al. (2013). Unfortunately, the model suffers from excessive generalisation of conductivity of bedrock layers, and, therefore, the simulation results should be considered as very rough estimations. On the other hand, an authoritative steady-state flow model of Quaternary aquifers, based on isotope-hydrochemical data, has been compiled by Mokrik et al. (2014) for the eastern part of Lithuania.

The first version of the time-dependent groundwater flow and transport model of the Estonian Artesian Basin (EAB) and its close surroundings was completed by Vallner (2003). The model was named Hydrogeological Model of Estonia (HME hereinafter) and its updated versions have been applied to various water problem inventory analyses and to scientific research (NGI 2004; MEE 2006; Gavrilova et al. 2010; Marandi & Vallner 2010; Vallner et al. 2015a, 2015b). The data and software of the model are constantly being upgraded to enhance the adequacy of simulations.

The goal of this paper is to characterise the HME, to review the main modelling results and to propose further developments for a holistic research of the Estonian hydrosphere.

The EAB was the main object of the completed regional groundwater modelling. The areal extent of the EAB, encompassing 70,000 km², is defined by the outer borders of the groundwater flow systems formed in the territory of Estonia (45,000 km² in area). The EAB includes the Quaternary cover (Q), a monocline of sedimentary Palaeozoic strata (Pz) dipping southward at a general rate of 2–4 m km⁻¹ and the underlying portion of the crystalline Proterozoic Baltic Shield (PR). The EAB can be considered as a northern structural element of the larger Baltic Artesian Basin. The latter is the sole European artesian basin singled out for the worldwide fundamental research programme Groundwater and Global PalaeoSignals (GGPS 2015).

Below come the Ordovician-Cambrian (O-Cm) and Cambrian-Vendian (Cm-V) aquifer systems, consisting mostly of sandstone with siltstone interbeds. These two systems are separated by the Cambrian aquitard (Cm_a). The Cambrian-Vendian aquifer system (CVAS), including in the eastern part of the study area the upper, Voronka aquifer (Vr), the lower, Gdov aquifer (Gd), and the intermediate, Kotlin aquitard (Kt), crops out along the northern coast of Estonia on the bottom of the Gulf of Finland. The CVAS is the most important source of drinking water for North Estonia.

Cracks and pores of the up to 100 m thick upper portion of the crystalline basement (PR) contain a certain amount of water, which participates in the basin-wide groundwater flow. The depth of the basement from the ground surface increases from 100 to 150 m on the shore of the Gulf of Finland to 500–800 m along the southern border of the model area.

The lateral hydraulic conductivity of sandstones, limestones and dolomites is usually between 2 and 8 m day^–1, the storage coefficient ranges mostly from 10^–5 to 10^–3 m⁻¹. The lateral conductivity of karstified carbonate rocks can reach 50 m day^–1 or even more. The vertical hydraulic conductivity of aquitards mostly varies in an interval of 10^–9–10^–2 m day^–1.

The upper strata of the EAB usually contain fresh water, but in deeper strata total dissolved solids (TDS) values increase southward, up to 22 g L⁻¹ (Vallner 1994). The content of fluoride exceeds the Estonian sanitary limit (1.5 mg L⁻¹) in S-O layers in Central and Western Estonia (Karro & Marandi 2003). Groundwater often comprises too much iron and requires special treatment before use. In Northeast Estonia, mining and oil shale processing have caused large-scale contamination of groundwater and surface water bodies by sulphides and phenols (Gavrilova et al. 2010; KM 2011; Vallner et al. 2015a).

During the last seven decades, intensive groundwater abstraction for water supply and mine dewatering totalling up to 600,000–800,000 m³ day⁻¹ has taken place in the model area. Therefore, basin-wide head depressions have formed in several aquifers, engendering the intrusion of unpotable water toward groundwater intakes (Vallner & Järvet 1998; MEE 2006; Marandi & Karro 2008; Gavrilova et al. 2010; Marandi & Vallner 2010). The groundwater public supply was about 120,000 m³ day⁻¹ and the water removal from quarries and underground mines was 460,000 m³ day⁻¹ in Estonia in 2014.

The HME involves all main aquifers, aquifer systems and aquitards marked by hydrogeological indices in the previous section. The units are represented by 13 model-layers in the basic version of the HME. Some versions of the HME dedicated to handling sophisticated problems of water supply or dewatering of mines comprise up to 22 model layers. Additional model layers were produced by splitting some principal hydrogeological units along the second-rate aquitards. The database used for modelling contains characteristics of about 16,000 boreholes and water wells, monitoring data of hydraulic head and groundwater quality. The results of about 1,000 time-drawdown pumping tests were used to characterise the hydraulic conductivity, the specific storage and the specific yield of the water-bearing strata.

The top boundary of the model coincides with the ground surface or the bottom of streams, lakes and the sea. The surface lying at a depth of 100 m beneath the upper surface of the crystalline basement acts as the supposedly impermeable bottom boundary of the model. The thickness of the water-bearing formation modelled varies from 100 m on the northern border of the model to 600–900 m along its southern border (Figure 3).

The northern border of the model area practically coincides with the updip limit of the EAB sedimentary bedrock on the bottom of the Gulf of Finland. The western border follows the general flow line of Cambrian water stretching northwards across contours of the hydraulic head. Albeit the eastern border of the EAB actually runs along the Narva River and the central part of Lake Peipsi and Lake Pskov, the border of the model area has been shifted eastward by 10–80 km. This allowed taking into consideration the entire catchment of the Narva River, which was essential in solving water management problems between Estonia and Russia.

The southern border of the model, between Estonia and Latvia, was established in light of the same aim.

The Dirichlet type of boundary condition prescribing the given value of the hydraulic head h(x, y, z, t) was mostly assigned to outer boundaries of the modelling domain according to data of hydrogeological mappings. The straight border lines applied did not affect the general authenticity of model outputs.

The modelling area was covered with an imaginary rectangular computational grid at spacing from 1,000 to 4,000 m for finite-difference discretisation of water-bearing layers in the regional scale. The initial and boundary conditions of the model were positioned and their numerical values assigned to the cells of the computational grid as input parameters of the HME.

The groundwater discharge in rivers R and pumping from layers P are the instrumentally measured terms in Equation (3); the direct seepage of groundwater to the sea M, and the subsurface flux V are calculated by Darcy's formula.

Net infiltration I is an important parameter for a reliable groundwater modelling. It is expressed by the term W in Equation (1). Due to positive value of net infiltration, divides of the groundwater table have formed supporting the groundwater flow to river network. The net infiltration is a main component of the groundwater inflow at the regional groundwater budget.

Direct measuring of net infiltration for wide areas by lysimeters is too expensive. It is usually possible only in restricted areas not exceeding a few square kilometres. Net infiltration of greater areas could be determined by variable formulas, elaborated on the basis of such indirect data as the air temperature, atmospheric humidity, wind speed, albedo, plant characteristics, etc. (Scanlon et al. 2002; Herrmann et al. 2009). Unfortunately, those kinds of calculations can often be quite unreliable. Their results are mainly acceptable for the same domains used for deducing empiric functions.

A collation of measured R and P values with M and V values calculated by Darcy's formula in Equation (1) shows that both the groundwater discharge in the channel network R and the groundwater pumping P make up about 90% of the sum of the right-hand-side in the equation. It enables a provisional estimation of net infiltration I with a relative error of about 10% through the sum of the base flow and groundwater pumping. The I values preliminarily estimated were put into the computer model of the HME. All terms of Equation (3) are inputs or outputs of the regional model. Their adjustment was continued in the course of the model calibration.

Water exchange between aquifers and large surface waterbodies (the sea, Lake Peipsi and Lake Pskov) has also been simulated by the River Package of Visual MODFLOW Classic for the HME.

The groundwater outflow from Ordovician layers along the North Estonian Cliff, and the discharge of the Silurian-Ordovician aquifer system through numerous springs were modelled by modified Cauchy boundary conditions represented by the Drain Package in Visual MODFLOW Classic.

The initial condition needed for proceeding of time-dependent simulations was established by deactivating all groundwater intakes put into the preliminary calibrated steady-state model. A corresponding distribution of the heads h(x, y, z) was simulated and considered as a mathematically correct description of the initial condition for the transient transport model. The versions of the HME obtained for the steady-state flow regime (∂h/∂t = 0) describe the three-dimensional distribution of hydraulic head in the modelling domain at t = 0. This moment of the computational time corresponds to the year 1910, marking the beginning of intensive groundwater abstraction in Estonia.

The subprogram MT3DMS, coupled with the code Visual MODFLOW Classic, was also activated in the HME to study both local and basin-wide problems of groundwater transport. For that purpose, the geometry of the transport boundaries was defined and the concentrations of groundwater ingredients, layers’ dispersivity, characteristics of chemical reactions and other necessary parameters were input into the HME. The procedures will be more closely described hereinafter in sections ‘Palaeohydrological reconstructions’ and ‘Environmental risks of oil shale semi-coke’.

The long-term groundwater discharge in streams has been estimated based on observations carried out during several decades at more than 100 hydrological gauging stations distributed quite evenly all over the modelling area (Figure 2). The observation series have been perfectly studied by methods of mathematical statistics and the mean long-term low base flow during 30 days Q has been defined for gauging stations (Protas'eva & Eipre 1972). Apart from gauging stations, many irregular measurements of the low flow have been made in approximately 1,000 stream cross sections. The obtained sporadic low flow data have been modified to an average long-time low base flow by statistical methods based on regular observations at gauging stations. Detailed maps of the low base flow at a scale of 1:200,000 for the modelling area were compiled by L. Vallner (ETAGI 1970; IGANE 1975).

To estimate the net infiltration I, the study area was divided into a number of calculation zones, which coincided with catchments of hydrological gauging stations. The I value was found by Equation (3) for every calculation zone and the corresponding map of the net infiltration was completed (Vallner 1976, 1997). The map was used for the input of groundwater recharge data into the HME.

The value of the parameter R_c needed for calculating the conductance by Equation (5) was obtained from the aforementioned base flow map and the heads H and H_r were estimated on the basis of the hydrogeological mapping carried out in the course of an official geological survey. In this way, a feasible calculation result for the riverbed conductance was obtained.

As mentioned above, large lakes and the sea were modelled as rivers; the conductance of the bed material was estimated for those as well as for real rivers. Instead of the value of R_c, the direct seepage of groundwater to a large surface waterbody M was used in Equation (5). The value of M was calculated by Darcy's formula and checked by Equation (3).

The pumping data were obtained from state institutions monitoring the groundwater use. All significant groundwater intakes with their abstraction rates were taken into account and incorporated into the HME. The results of about 1,000 time-drawdown and distance-drawdown pumping tests were used to characterise the hydraulic conductivity, the specific storage and specific yield of water-bearing layers.

The data on groundwater heads and water quality were obtained from state institutions dealing with the monitoring of the water environment.

The trustworthiness of the HME has been assessed by comparing the measured and calculated data. If too large a deviation was recorded, the model parameters were corrected until the deviation no longer exceeded the permitted range. In general, the calibrated version of the HME was accepted if the correlation coefficient between the measured and modelled data was more than 0.8. The calibration process, besides its checking and correcting functions, was also used for determination of the model's parameters, whose direct experimental estimation was very difficult, if not impossible, in practice. For that purpose, methods of inverse modelling were used (Schwede & Cirpka 2010). The adequacy of such corrections to model parameters was verified by secondary model calibrations.

The current HME flow model has been calibrated against two different sets of calibration targets. One of the sets represented the statistically established mean long-term value of the lowest base flow during 30 days Q_m in the hydrological gauging stations, and the other one corresponded to the mean long-term measured elevations of the groundwater table and heads during low base flow periods in the representative monitoring points. The Q_m-target was chosen to exclude the possibility that a portion of the surface runoff could be involved at calibration. On the other hand, it was presumed that the conductivity of layers calibrated for the low flow period would be the same in medium and high flow periods.

The HME is an indispensable tool to investigate the regional-scale problems of the Estonian water environment. The model provides scientifically well-founded authentic data ensembles of the EAB compound elements. It renders it possible to quantify the time-dependent three-dimensional hydraulic and hydrogeochemical processes occurring in the HME, as well as their mutual relationship with the surface waterbodies and climatic factors. Flows between the adjacent water-bearing layers, the channel network and the sea (Figure 3), and changes in water quality can be determined both in natural and artificial conditions. The model enables simulation of the drawdown caused by groundwater abstraction from dissimilar layers and in different places. It supports the compilation of long-term plans for optimum water management and protection. Some examples of the regional characteristics of the EAB simulated by the HME are discussed below.

The hydrogeological setting of HME is appropriate for assembling a regional water balance system. Seas and large lakes from three sides surround the water-bearing formation here. A comparatively dense network of gauging stations on rivers controls the groundwater runoff. All this helps bring about an adequate specification of boundary conditions for a correct determination of main components of the regional water budget. A concise water budget of the EAB itemised by main hydrogeological units, is shown in Table 1. The main budget units can be articulated to a sequence of hierarchic subunits to consider the local budget problems. A similar water budget has also been completed for the entire modelling area.

Determination of volumes of hydrogeological units by means of the code ArcGIS enabled estimation of the approximate amount of groundwater in layers and assessment of the duration of its total exchange (Table 2).

Local water budgets can be completed independently from the EAB, but despite their virtual correctness, there is no guarantee that they always match a higher-rank budget system. A discordance between neighbouring budget units shows that the values of some flow components calculated are not correct. The problem can be overcome by completing a hierarchic water budget system, where the budget units of the lower rank frame into units of the higher rank; the latter, in turn, are suited for a general water budget system based on the HME. A regional water budget system containing subsystems reciprocally balanced minimises possible errors at establishing groundwater flow rates. It is very important for determination of an optimum exploitation regime of aquifer systems (Custodio 2002; Devlin & Sophocleous 2005; Bredehoeft & Durbin 2009).

The evaporation from groundwater table for the EAB can also be estimated by the HME.

The practical concurrence of measured and simulated values of the base flow proved by results of calibration (Figure 5) show that the HME could be used as the base flow model of the EAB. It enables the areal estimating of the base flow intensity, the quantifying of the base flow and a computerised coupling of base flow data with the hydraulic parameters of water bearing layers and climatic elements. The seasonal dynamics of the base flow and its share in the total surface runoff can be determined.

As pointed out above, the BAB along with the EAB are outstanding objects of the worldwide scientific research on palaeoclimatology (GGPS 2015). The isotope composition of groundwater indicates that the groundwater of the EAB has been strongly influenced by the advancing and retreating ice sheet during the last continental glaciation (Vaikmäe et al. 2001, 2008). However, the time-dependent mechanism of such an influence is still disputable, as well as the spatial extent of palaeowater that has intruded or remained in the layers.

The time-dependent formation and spatial distribution of groundwater's genetic types during the last 22 ka were reconstructed by a series of simulations carried out by means of the HME (Vallner et al. 2015b). For that purpose, it was assumed that the time-dependent spatial distribution of ¹⁸O concentration in groundwater could be simulated by code MT3DMS as a process of molecular diffusion. MT3DMS does not accept negative values of boundary conditions except a case connected with selection of model outputs desired (Zheng & Wang 1999). Therefore, the negative values of δ¹⁸O were recalculated to the respective positive values of ¹⁸O mass concentrations having dimension mg L⁻¹. Representative groundwater heads induced by the load of the ice sheet as well as the presumed values of ¹⁸O mass concentrations of glacial melt water, meteoric water and brines were put into the HME as initial and boundary conditions.

According to initial conditions of the time-dependent δ¹⁸O transport model, the upper layers of the model area down to the O-Cm aquifer system contained meteoric water with the δ¹⁸O concentration equal to −10‰ at 22 ka BP. At the same time, the underlying layers were filled with brine having the δ¹⁸O value of −6‰.

Transport simulations of ¹⁸O as a conservative tracer were performed for all the main stages of the last glaciation and deglaciation. The advection term of Equation (2) was solved by implicit finite-difference method using upstream weighting at a Courant number of 0.75 (Zheng 2010). Since the MT3DMS calculations are based on a steady-state flow model, the period under consideration was modelled by seven single steady-state groundwater flow complementary models coupled with MT3DMS transport code representing the phases of glaciation, deglaciation, and the evolution of the Baltic Sea. Seven sequential simulations were performed in which the distribution of ¹⁸O concentration simulated by the previous model was imported into the next. This way, all complementary models were functionally connected to each other and, together, they formed an integrated calculation unit suitable for interpreting the successive palaeogeographical events.

The simulations of the ¹⁸O transport were based on general conception about the development of the last glaciation and its degradation in Northeast Europe, presented by Kalm et al. (2011) and Siegert (2001). The modelling showed that groundwater heads and flow directions were completely rearranged because of the changing ice load in the EAB. Between 22 and 18 ka BP, the continental ice advanced in southern and southeastern directions, covering the EAB completely.

At 18 ka BP deglaciation began and continental ice perished from the area of the EAB at 12 ka BP. After that, the ice lakes and the Baltic Sea flooded the western lowlands of the present-day Estonia until 2 ka BP. During deglaciation and flooding of the model area, the meteoric water began to reinfiltrate, intruding the layers filled earlier with glacial meltwater, thus producing a complicated mixture of waters of different genetic types. All main stages of glaciation, deglaciation and evolution of the Baltic Sea were simulated by the HME for the time interval from 22 until 0.1 ka BP.

According to the requirements of the European Commission, the long-term average groundwater abstraction must not exceed the available groundwater resources, which express the long-term natural average rate of the overall recharge of a groundwater body (Water Directive 2000). If the real abstraction from a groundwater body exceeds its available resources, both quantity and quality of the surface water hydraulically connected with this groundwater body can suffer. Alteration of groundwater flow direction must not cause the intrusion of saline water or other harmful substances into a groundwater body. To check the fulfilment of these criteria, the available resources of all main groundwater bodies of Estonia were simulated by HME. A detailed analysis of the regional groundwater balance was completed, which proved that the status of deep aquifer systems in North Estonia was not sustainable (MEE 2006; Gavrilova et al. 2010).

The results of the study show that oil shale mining is a major pressure on water resources in Northeast Estonia. Because of mine dewatering, the groundwater table has lowered along the perimeter of the mines by 10–70 m to the roof level of the goaf. Water table depressions extend up to 6–7 km from the operating mines. The amount of water pumped from oil shale mines exceeds the available resources of the Ordovician groundwater body by five times. The mine water is pumped into the channel network without any chemical purification. Therefore, sulphates, petroleum products and phenols in the region of the oil shale industry have contaminated the stream water.

Transport simulations (Marandi & Vallner 2010) investigated the reason for the sporadic increase in the salinity over the permitted limit in the CVAS exploited in Tallinn, the capital of Estonia. The study area comprising 9 km² was functionally coupled with the HME by the local refining of the HME's computational grid to the spacing of 250 × 250 m. It enhanced the authenticity of the local model. Taking into account the density-dependent transport of solutes based on codes MT3DMS and SEAWAT-2000 (Weixing & Langevin 2002; Zheng 2010), it was proved by simulations that the observed increase in the TDS was caused by upconing of salty brine from the crystalline basement due to overexploitation of the CVAS. Some measures were proposed to reduce or avoid this process by proper design of production wells.

The Kohtla-Järve oil shale semi-coke and ash landfill, which is likely the largest of its kind in the world, was started in 1938. No artificial barrier impeding the vertical movement of landfill leachate was arranged on the dumping site. Approximately 10⁸ t or 8 × 10⁷ m³ of waste were stored in an area of 2.14 km² until 2010. According to the EU Landfill Directive (1999), Estonia's Ministry of the Environment initiated the implementation of measures to alleviate the presumptive contamination impact of the landfill in 2003. It was decided to close approximately half of the landfill, where depositing of waste had mostly ceased decades before, to smooth the ground surface of the closed area and to cap it partially with a presumably waterproof liner. It was declared without any scientific argumentation that the planned remediation actions would significantly lessen the groundwater contamination.

The work started in 2010; however, the earth-moving operations caused self-ignition of the landfill in 2012 because of intensive oxidation of the organic matter contained in semi-coke. The active surface fire was extinguished, but even to date the temperature sporadically reaches 400 °C in deeper layers of the landfill.

Owing to the situation, the Estonian public raised the question whether the capping action of the Kohtla-Järve landfill site had been scientifically and economically sufficiently grounded, or had it been an incompetent project and a large-scale misuse of international funds (Haravee 2012; Niitra 2012; Strandberg 2012)?

This problem was investigated by means of HME (NGI 2004; KM 2011; Vallner et al. 2015a). A local groundwater flow and transport model of the landfill's area was completed by a sectional refining of the computational grid of the HME to the spacing of 125 × 125 m. Such linking contributed to motivation of boundary conditions of the local model and enabled checking its matching with the regional water budget system. The size of the rectangular study area was 5.3 × 5.0 km covering 26.5 km².

The TDS was considered as a general contamination indicator and as a conservative tracer not participating in chemical reactions. Concentration of other non-reactive groundwater ingredients might presumably change proportionally to the TDS at groundwater transport. On the background of the initial TDS equal to 500 mg L⁻¹ for the study area, the TDS values exceeding this threshold were handled as evident signs of groundwater contamination. The content of phenols was the most suitable parameter for characterising the transport and fate of organic contaminants subjected to processes of sorption and biological degradation in the study area. Leaching of phenols with other organic contaminants from the dump deposits of the landfill site started in 1938.

The transport models were calibrated against the observed concentrations of groundwater ingredients in the monitoring points. Estimation of parameters for transport models was quite complicated. For instance, no special tests were done for experimental determination of the porosity and hydrodynamic dispersion of the groundwater environment in the study area of the Kohtla-Järve oil shale semi-coke landfill. The sole data for the construction of an adequate transport model were obtained from water quality monitoring performed in 1996–2010 (MAVES 2007, 2008; KM 2011). Therefore, the simulation of the subsurface contamination plume was mainly based on presumed data; later, the simulated concentration data were compared with the observed concentration data at monitoring points by model calibration (Vallner et al. 2015a). The model conductivity and transport parameters were adjusted during calibration. Afterwards, the spatial time-dependent spreading and fate of contaminants were predicted, based on model parameters specified by the calibration. In such a way, the calibration of the transport model was changed from an adjusting and checking operation into a principal method for establishing model parameters.

The hydrogeochemical parameters of the phenol model specified by its inverse calibration are as follows: the longitudinal, horizontal and vertical dispersivity are 50 m, 5 m and 0.5 m, respectively; the coefficient of molecular diffusion is 0.1 m² day⁻¹, the distribution coefficient is 1.1 × 10⁻⁶L mg⁻¹; first-order decay rate for the dissolved phase is 10⁻³ day⁻¹; and the same for the sorbed (solid) phase is 10⁻⁴ day⁻¹, and the average bulk density of model layers is 1.7 × 10⁶ mg L⁻¹.

Based on the transport modelling described above, it was established that the landfill has been a serious source of contamination for the subsurface throughout the past century, but the plumes of organic and inorganic contaminants have currently stabilised in the groundwater. The size of the subsurface contamination plume, which covers an area of 8.3 km², will increase during the next century by only a certain percentage because of biodegradation of organic contaminants and due to attenuation of contaminated fluxes (Figure 12). The expansion of the pollution plume is not essential from the viewpoint of practical environmental protection.

Current emission of gases from the landfill into the atmosphere is the most significant hazard presently, since every inhabitant of the Kohtla-Järve region breathing the contaminated air suffers from it to some extent. A lesser environmental hazard is contamination of water in the channel network of the Kohtla River by leachates of the landfill's dump deposits, which are toxic for the aquatic biota.

Thus, the groundwater transport simulations showed that capping of the Kohtla-Järve landfill has not been sufficiently substantiated. The landfill liner used is not able to restrict the expansion of pollutants in groundwater. All expenditures for accomplishing this landfill capping financed by Estonian and EC taxpayers have been useless.

The EU Landfill Directive (1999) requires that the base and sides of a landfill for hazardous waste shall consist of a mineral layer for which the hydraulic conductivity K does not exceed 1.0 × 10⁻⁹ m s⁻¹ and the thickness is at least 5 m. The simulations based on the HME demonstrated that fulfilling these criteria does not prevent the outflow of leachate from the landfill.

The HME is continuously upgraded to enhance the adequacy of the simulation results and to expand the spectrum of research topics under consideration.

The geometry of the model layers has to adjust to the continuously refined spacing of the virtual computational grid. It is especially necessary for constructing additional local models hydraulically connected with the regional model. The account of the impact of collateral aquitards or semi-pervious layers is needed to solve many local problems of water supply and environmental protection. Therefore, the basic units of the hydrogeological stratigraphy should be split into finer seams. In Estonian databases, there are enough records for such corrections.

Recently, the geometry of the HME's top layer was significantly elaborated by Porman (2015). The airborne LIDAR data of the State Digital Elevation Model with the vertical accuracy of ±0.34 m were used to adjust the elevation of the ground surface of the HME. The modern data on seabed elevation were also used. The most difficult task was to collect thousands of borehole logs in many different formats scattered between several archives. However, a digital database on the geometry of the Quaternary cover was created for the first time and its data were put into the HME. Since the geometry of the Quaternary deposits was corrected, the riverbed thickness and the river stage also needed adjustment in places. Therefore, a new series of the HME calibrations was conducted.

The hydraulic and transport parameters of layers also need specifying, but, unfortunately, new comprehensive field experiments are rare. To compensate for the lack of adequate hydraulic and hydrogeochemical parameters, special series of inverse modelling and calibration are performed. The results are used for upgrading the model.

Usable information can be obtained from various groundwater head and water quality observations, which are continued or started in new places or include additional layers. Therefore, the records of isotopic contents are of special interest, insofar as these data can be used as tracers of the water movement.

To date, a renewal of the basic software of the HME is taking place. Most of the existing versions of the HML based on Visual MODFLOW Classic will be transferred to Visual MODFLOW Flex 2014.2 (Schlumberger 2015b). The computer code generating an unstructured computational grid facilitates the detailed investigation of detached objects on the background of the regional model. The pinching out of layers can be taken into account with a satisfying adequacy. These features will make the modelling more flexible and increase the authenticity of model versions.

The HML based on Visual MODFLOW Flex 2014.2 will be coupled with MIKE SHE (DHI 2014) or with another effective computer code aimed at a sophisticated hydrological modelling. Such a system connected with regional databases can be used for the analysis, planning and management of a wide range of water resource and environmental problems (Refsgaard et al. 2010). Interconnection between channel water, groundwater and the sea, the conjunctive use of groundwater and surface water, wetland management and restoration, the land use and climate can be studied in their mutual relationship. The main goal of the hydrological development of the HME is to create an integrated and universal computer model of the Estonian hydrosphere.

The sophisticated research and management problems of the EAB including the territory of Estonia should be considered by means of a holistic modelling of the water environment. The regional groundwater flow and transport model encompassing the entire EAB, and border districts of Russia and Latvia has been completed based on the code Visual MODFLOW Classic. The model covers about 88,000 km² and involves all main aquifers and aquitards from ground surface to as low as the provisionally impermeable part of the crystalline basement.

The model development is based on quantifying of the relationship between groundwater fluxes and surface waterbodies. Detailed base flow data specified by additional field measurements in conjunction with profound water budget investigations were used to adjust the model input parameters. The model has been properly calibrated against the observed values of groundwater heads and statistically assessed rates of the low base flow.

The model is a mighty and trustworthy tool for the advanced study of water problems. The main achievements of the EAB modelling are as follows:

• The time-dependent three-dimensional distribution of groundwater heads and main hydrogeological parameters of layers have been determined.

• The direction, velocity and rate of subsurface fluxes have been estimated.

• The itemised regional water budgets, durations of groundwater exchange, and base flow rates have been simulated.

• The principal problems of the sustainable management of the water environment have been elaborated.

• The palaeohydrological situation during the last continental glaciation and deglaciation has been reconstructed proceeding from three-dimensional transport simulations of ¹⁸O concentration in groundwater – this is a scientific innovation for tracer application in hydrogeological methodology.

• It has been proved that hydraulic criteria of an effective geological barrier for a landfill of hazardous waste posed by the EU Landfill Directive are too mild and need revision – this finding should be taken into consideration in all EU member states.

The model is continuously upgraded to enhance the adequacy of the simulation results. The main goal of the further hydrological development of the current model is to create an integrated and universal computer model of the entire Estonian hydrosphere.

This paper is a contribution to the research project IUT19-22 ‘Groundwater flow history, global palaeoclimate signals and anthropogenic influence in the Baltic Artesian Basin: a synthesis of numerical models and hydrogeochemical data’, financed by the Estonian Research Council, and to the international programme ‘Groundwater and Global PalaeoSignals’. The authors are very grateful to the institutions for the assistance rendered. This paper benefited greatly from suggestions made by two anonymous reviewer.