The sophisticated research and management problems of the Estonian water environment are considered by means of holistic modelling. The model elaborated is based on the computer code Visual MODFLOW Classic. The model encompasses the entire Estonian Artesian Basin (EAB) and the border districts of Russia and Latvia. It involves all main aquifers and aquitards on an area of 88,000 km2. The main hydrogeological and hydrological characteristics of the study area, including the time-dependent three-dimensional distribution of groundwater heads, the direction, velocity and rate of subsurface fluxes, itemised water budgets, volumes of hydrogeological units, and durations of groundwater exchange have been determined by modelling. The palaeohydrological situation during the last continental glaciation of the EAB was reconstructed and the principal problems of the sustainable management of water environment were elucidated. The model has been used to simulate the local and cumulative rates of the base flow.

INTRODUCTION

The state of groundwater is complicated and causes anxiety in Estonia (Figure 1). Upper aquifers are suffering from industrial, agricultural and former military pollution. In deeper aquifers, a number of large depression cones have formed due to intensive withdrawal from wells, which induces harmful intrusions from artificial pollution sources, sea and natural bodies of unpotable groundwater towards groundwater intakes (Vallner 1996; Vallner & Järvet 1998; Karro & Marandi 2003; MEE 2006; KM 2011; Vallner et al. 2015a). The state of the water environment has seriously worsened in Northeast Estonia due to oil shale mining and processing. The natural quality of groundwater is poor in several main aquifers.

Figure 1

Location of the study area.

Figure 1

Location of the study area.

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.

HYDROGEOLOGICAL SETTING

The model area, about 88,000 km2 in total, includes the territory of Estonia with the surrounding portions of the Baltic Sea, the Gulf of Finland, Lake Peipsi with Lake Pskov and border regions of Russia and Latvia (Figures 1 and 2). The average absolute height of the flat terrestrial part of the area is about 50 m, and only a few places are 150–250 m above sea level. The Baltic Sea, with the Gulf of Finland, is the main drainage basin, Lake Peipsi with an elevation of 30 m a.s.l. being the second largest. The climate is moderately cool and humid. Average annual precipitation is 500–750 mm. The total surface runoff from the territory is about 270 mm year−1.

Figure 2

Model area and EAB.

Figure 2

Model area and EAB.

The EAB was the main object of the completed regional groundwater modelling. The areal extent of the EAB, encompassing 70,000 km2, is defined by the outer borders of the groundwater flow systems formed in the territory of Estonia (45,000 km2 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−1 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).

Quaternary deposits (Q), consisting predominantly of glacial till and glaciolacustrine sandy loam, form the uppermost aquifer system in the model area (Perens & Vallner 1997). Their thickness usually ranges from a couple of metres to 30 m, occasionally reaching 100–150 m. In the southern and eastern parts of the model area, Quaternary deposits cover Devonian (D) aquifer systems (D3, D2, D2-1), including the aquitard D2Nr, which are represented chiefly by sandstone and siltstone (Figure 3). In North and Central Estonia and on the West Estonian islands, the Quaternary deposits lie on the outcrop of the Silurian-Ordovician aquifer system (S-O) consisting of limestones and dolomites. Frequently, the upper part of the formation is heavily karstified and cavernous. Deeper than 100 m from the bedrock surface the fissures are almost closed in carbonate strata, which turn into the Silurian-Ordovician aquitard (S-Or) of regional extent.

Figure 3

Hydrogeological cross-section A–B (section line in Figure 2).

Figure 3

Hydrogeological cross-section A–B (section line in Figure 2).

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 (Cma). 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−1. 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−1 (Vallner 1994). The content of fluoride exceeds the Estonian sanitary limit (1.5 mg L−1) 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 m3 day−1 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 m3 day−1 and the water removal from quarries and underground mines was 460,000 m3 day−1 in Estonia in 2014.

CONCEPTUAL MODEL

Computer codes

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 current version of the three-dimensional transient groundwater flow and transport model HME is based on the computer code Visual MODFLOW Classic (Schlumberger 2015a). In this code, the distribution of hydraulic head in a porous heterogeneous domain containing groundwater is described by the governing partial differential equation (Harbaugh 2005; Fetter 2014): 
formula
1
where Kxx, Kyy and Kzz are values of hydraulic conductivity along the x, y and z coordinate axes [LT−1]; h is the hydraulic head [L]; W is volumetric flux per unit volume and represents sources and/or sinks of water [T−1]; Ss is the specific storage of the porous material [L−1]; and t is time [T].
Visual MODFLOW Classic is coupled with the codes SEAWAT and MT3DMS; the former enables simulation of the head h in a variable-density groundwater system (Weixing & Langevin 2002), the latter is a solute transport model for handling the broad range of equilibrium and kinetically controlled geochemical processes (Zheng & Wang 1999; Zheng 2010). The MT3DMS code is based on the fundamental differential equation describing the transport of groundwater ingredients in three-dimensional transient groundwater flow system (Zheng & Bennett 2002): 
formula
2
where Ck is the dissolved concentration of a groundwater ingredient with dimension [ML−3]; Θ is the porosity of the subsurface medium (dimensionless); xi is the distance along the respective Cartesian coordinate axis [L]; Dij is the hydrodynamic dispersion coefficient tensor [L2T−1]; νi is the linear pore water velocity [LT−1]; qs is the volumetric flow rate per unit volume of aquifer representing fluid sources (positive) or sinks (negative) [T−1]; Csk is the concentration of the source or sink flux [ML−3]; and ∑Rn is the chemical reaction term.

Boundary conditions and their parameters

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 most complicated problem at defining the boundary conditions was determination of the intensity of the groundwater areal recharge, which was simulated in Visual MODFLOW by the Recharge Package. To solve this issue, a fundamental equation of the regional groundwater budget was completed as follows: 
formula
3
where I is the net infiltration or the total groundwater recharge minus evaporation from the zone of saturation or capillary fringe; R is the groundwater discharge in rivers (the base flow); P is the pumping from layers; M is the direct seepage of groundwater to the sea; V is the subsurface exchange of groundwater between the HME and the surrounding region; and S is the storage change.

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.

To simulate the flux between a river and a groundwater system by the River Package based on the Cauchy boundary condition, the Visual MODFLOW Classic requires the input into the model of a special parameter C, called the conductance, and defined by the equation: 
formula
4
where K is the vertical hydraulic conductivity of the riverbed; L is the length of the river through a cell; W is the width of the river in the computational cell; and M is the thickness of the riverbed.
In general, a direct estimation of the vertical hydraulic conductivity K of a riverbed separating the surface waterbody from the underlying groundwater system is quite troublesome. Usually, there is not enough field test data for correct calculations. Therefore, the conductance C has been defined for the HME, proceeding from the equation: 
formula
5
where Rc is the groundwater discharge into the river in the computational cell under consideration; H is the average value of the groundwater head in the same cell; and Hc is the corresponding mean river stage elevation.

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’.

DETERMINATION OF FLOW MODEL INPUT PARAMETERS

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 Rc needed for calculating the conductance by Equation (5) was obtained from the aforementioned base flow map and the heads H and Hr 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 Rc, 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.

FLOW MODEL CALIBRATION

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 Qm 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 Qm-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.

Calibrations were completed for several different periods of extreme precipitation amounts or groundwater abstraction. Data sets of 1976, 1990, 1998, 2006 and 2013 were used (Figure 4). Using different calibration targets and different data sets has enhanced the reliability of the modelling results, especially the adequacy of budget calculations. The manual trial-and-error adjustment of model parameters was used at calibration to achieve the optimum match between simulated parameters and calibration targets. It was closely observed that the corrections to model parameters introduced by calibration would remain within the ranges of the hydrogeological reality.

Figure 4

Total pumpage from CVAS (curve 1) and change of the potentiometric surface in observation wells 705 (Tallinn, curve 2), and 864 (at Kohtla-Järve, curve 3).

Figure 4

Total pumpage from CVAS (curve 1) and change of the potentiometric surface in observation wells 705 (Tallinn, curve 2), and 864 (at Kohtla-Järve, curve 3).

The conductance C of a river cell was mostly modified at model calibration relying on the base flow data. The river stage Hc and the groundwater head H were also corrected when necessary. When adjustments did not help, the given value of the net infiltration I in the catchments was altered. Calibrations were continued until the correlation coefficient between measured and simulated data exceeded 0.9 (Figure 5). Simulated elevations of the groundwater table and heads were rechecked against statistical assessments of the field data series, until the difference between computed and measured values was less than 3.5 m in the regional model. While calibrating local models, the allowed discrepancy between measured and simulated groundwater heads was decreased to 1.0 m. The correlation coefficient between measured and simulated groundwater heads exceeded 0.9 and practically all calibration points fell into the 95% interval of the statistical exception at main groundwater intakes (Figure 6).

Figure 5

Relationship between statistically assessed low base flow Qm and simulated base flow Qs for river gauging stations.

Figure 5

Relationship between statistically assessed low base flow Qm and simulated base flow Qs for river gauging stations.

Figure 6

Relationship between measured and simulated mean annual groundwater heads for main groundwater intakes in 1991.

Figure 6

Relationship between measured and simulated mean annual groundwater heads for main groundwater intakes in 1991.

MODEL APPLICATIONS

General characteristics of the EAB

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.

Table 1

Groundwater flows of the EAB in predevelopment conditions (thousands of m3 day−1)

Hydrogeological unit Inflow (Qin)
 
Outflow (Qout)
 
Water exchangea 
Lateral From above From below Lateral Up Down 
4.1 6,809.1b 4,345.8 15.2 6,255.6c 4,888.2 11,159 
     762.1d   
D3 3.1 194.2 <0.1e 16.6 <0.1 180.6 197.2 
D3-2 77.2 2,562.1 0.2 464.2 1,943.4 231.9 2,639.5 
D2Nr 20.9 374.6 251.5 1.8 364.3 281 647.1 
D2-1 48.5 281.7 101 22.5 330.4 78.3 431.2 
S-O 100 2,592.9 234.5 74.6 2,679.3 173.5 2,927.4 
S-Or 0.1 90 33.1 0.1 38.3 84.9 123.3 
O-Cm 15.5 84.9 12.8 25.9 64.8 22.5 113.2 
Cmr <0.1 21.7 12.1 0.2 12.4 21.2 33.8 
Vr 20.1 12.8 2.6 24.8 14.5 41.9 
Kt 0.3 8.9 7.4 0.3 7.3 8.9 16.5 
Gd 15.2 15.5 3.9 6.3 25.2 3.1 34.6 
PR 1.5 3.8 0.1 5.2 5.3 
Hydrogeological unit Inflow (Qin)
 
Outflow (Qout)
 
Water exchangea 
Lateral From above From below Lateral Up Down 
4.1 6,809.1b 4,345.8 15.2 6,255.6c 4,888.2 11,159 
     762.1d   
D3 3.1 194.2 <0.1e 16.6 <0.1 180.6 197.2 
D3-2 77.2 2,562.1 0.2 464.2 1,943.4 231.9 2,639.5 
D2Nr 20.9 374.6 251.5 1.8 364.3 281 647.1 
D2-1 48.5 281.7 101 22.5 330.4 78.3 431.2 
S-O 100 2,592.9 234.5 74.6 2,679.3 173.5 2,927.4 
S-Or 0.1 90 33.1 0.1 38.3 84.9 123.3 
O-Cm 15.5 84.9 12.8 25.9 64.8 22.5 113.2 
Cmr <0.1 21.7 12.1 0.2 12.4 21.2 33.8 
Vr 20.1 12.8 2.6 24.8 14.5 41.9 
Kt 0.3 8.9 7.4 0.3 7.3 8.9 16.5 
Gd 15.2 15.5 3.9 6.3 25.2 3.1 34.6 
PR 1.5 3.8 0.1 5.2 5.3 

aWater exchange is inflow into a hydrogeological unit or outflow from it (Qin = Qout).

bNet infiltration.

cDirectly to channel network and sea.

dInto Quaternary deposits from the bedrock layer in the sea area.

eFlow rates less than 100 m3 day−1 have not been accounted at summation.

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).

Table 2

Groundwater amounts and exchange durations in the EAB

Hydrogeological unit Volume (km3Effective porosity Groundwater amount A (km3Annual water exchange E (km3 year−1Duration of total water exchange T = A/E (years) 
1,163 0.3 349.0 4.073 86 
D3 15 0.03 0.4 0.072 
D3-2 1,633 0.15 244.9 0.963 254 
D2Nr 1,347 0.01 13.5 0.236 57 
D2-1 1,331 0.15 199.7 0.157 1,269 
S-O 4,366 0.03 131.0 1.069 123 
S-Or 10,658 0.01 106.6 0.045 2,368 
O-Cm 4,371 0.05 218.5 0.041 5,289 
Cmr 1,598 0.01 16.0 0.012 1,295 
Vr 740 0.15 111.0 0.015 7,259 
Kt 414 0.15 62.2 0.006 10,323 
Gd 1,751 0.15 262.6 0.013 20,792 
PR 7,157 0.01 71.6 0.002 36,994 
Total 36,543  1,786.9   
Hydrogeological unit Volume (km3Effective porosity Groundwater amount A (km3Annual water exchange E (km3 year−1Duration of total water exchange T = A/E (years) 
1,163 0.3 349.0 4.073 86 
D3 15 0.03 0.4 0.072 
D3-2 1,633 0.15 244.9 0.963 254 
D2Nr 1,347 0.01 13.5 0.236 57 
D2-1 1,331 0.15 199.7 0.157 1,269 
S-O 4,366 0.03 131.0 1.069 123 
S-Or 10,658 0.01 106.6 0.045 2,368 
O-Cm 4,371 0.05 218.5 0.041 5,289 
Cmr 1,598 0.01 16.0 0.012 1,295 
Vr 740 0.15 111.0 0.015 7,259 
Kt 414 0.15 62.2 0.006 10,323 
Gd 1,751 0.15 262.6 0.013 20,792 
PR 7,157 0.01 71.6 0.002 36,994 
Total 36,543  1,786.9   

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 net infiltration is an essential component in the regional water budget. It has been determined by HME for the EAB as well as for the entire modelling area (Figure 7). The approximate local values of the net infiltration based on assessment by Equation (3) were put into the model. They were thoroughly checked and adjusted by calibration and inverse modelling and a general distribution of the intensity of the net infiltration I was simulated. It can be used to estimate the total evaporation E by the formula: 
formula
6
where N is the precipitation and G the total surface runoff.
Figure 7

Net infiltration (mm year–1) in the model area.

Figure 7

Net infiltration (mm year–1) in the model area.

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

The profoundly calibrated HME provides a unique chance to estimate the vertical hydraulic conductivity of strong aquitards. To carry out such an operation by pumping or other field tests would be extraordinarily troublesome and expensive. However, the vertical conductivity of aquitards can be calculated by trial and error, using its presumed values until the model calibration criteria have been met. In that way, the vertical hydraulic conductivity of the regional Cambrian aquitard was determined for the modelling area (Figure 8). The parameter is relevant to quantify the environmental risks connected with the vertical movement of contaminants and unpotable water in the aquifer systems of Estonia.

Figure 8

Vertical hydraulic conductivity (m day–1) of the regional Cambrian aquitard.

Figure 8

Vertical hydraulic conductivity (m day–1) of the regional Cambrian aquitard.

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.

Palaeohydrological reconstructions

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 18O 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 δ18O were recalculated to the respective positive values of 18O mass concentrations having dimension mg L−1. Representative groundwater heads induced by the load of the ice sheet as well as the presumed values of 18O 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 δ18O transport model, the upper layers of the model area down to the O-Cm aquifer system contained meteoric water with the δ18O concentration equal to −10‰ at 22 ka BP. At the same time, the underlying layers were filled with brine having the δ18O value of −6‰.

Transport simulations of 18O 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 18O 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 adequacy of scientific hypotheses on 18O concentrations used at simulations and the whole methodology of modelling were proved by the correlation coefficient equal to 0.86 between calculated and observed values of 18O determined for 185 monitoring points over the territory of the EAB (Figure 9).

Figure 9

Relationship between observed and simulated δ18O values.

Figure 9

Relationship between observed and simulated δ18O values.

The simulations of the 18O 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.

During this time, the infiltration of meteoric water characterised by δ18O equal to −10‰ was replaced by vertical intrusion of ice meltwater into subsoil (Figure 10). The meltwater with the δ18O value of −23‰ came into being because of an upward heat transfer and due to frictional heating of the moving glacier at the base of the ice sheet. The discharge of basinal brines with the δ18O value of −6‰ and meteoric water from the layers of the EAB into the depressions of the Baltic Sea and the Gulf of Finland was replaced with lateral encroachment of glacial meltwater into bedrock layers. The CVAS was mostly filled with glacial meltwater characterised by the δ18O value of −23‰. Deep groundwater, including basinal brines, partially raised upward being drained by network of subglacial currents.

Figure 10

Simulated δ18O values and groundwater genetic types in the cross-section C–D at 18 ka BP (section line in Figure 2).

Figure 10

Simulated δ18O values and groundwater genetic types in the cross-section C–D at 18 ka BP (section line in Figure 2).

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.

Optimum groundwater abstraction

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.

Over the last 60 years, groundwater abstraction has been up to three times greater than the available resources of deep aquifers. The result of excessive pumping allowed by the Ministry of the Environment was a basin-wide depression of groundwater head (Figure 11) inducing the saline seawater intrusion into the aquifers (Karro et al. 2004; Marandi & Karro 2008). In addition, in the past few decades, contamination of groundwater mainly by radium isotopes 226Ra and 228Ra has become evident in the Gdov aquifer. It may be caused by the upconing of radioactive nuclides from the crystalline basement due to heavy pumping from the overlying aquifers.

Figure 11

Simulated groundwater hydraulic heads in the Cambrian–Vendian aquifer system.

Figure 11

Simulated groundwater hydraulic heads in the Cambrian–Vendian aquifer system.

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 km2 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.

Environmental risks of oil shale semi-coke

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 108 t or 8 × 107 m3 of waste were stored in an area of 2.14 km2 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 km2.

The effect of complicated geochemical processes at oil shale waste disposal sites on the water environment of the Kohtla-Järve area was interpreted and predicted by groundwater transport simulations for the timespan from 1938 until 2113 (Figure 12). For that purpose, two detached time-dependent models were completed: one for modelling of conservative spreading of TDS in groundwater, and the other, to study the reactive transport of phenols.

Figure 12

Phenols’ concentration in a cross section through the oil shale semi-coke landfill in Kohtla-Järve.

Figure 12

Phenols’ concentration in a cross section through the oil shale semi-coke landfill in Kohtla-Järve.

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−1 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 m2 day−1, the distribution coefficient is 1.1 × 10−6L mg−1; first-order decay rate for the dissolved phase is 10−3 day−1; and the same for the sorbed (solid) phase is 10−4 day−1, and the average bulk density of model layers is 1.7 × 106 mg L−1.

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 km2, 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−9 m s−1 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.

MODEL DEVELOPMENT

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.

CONCLUSIONS

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 km2 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 18O 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.

ACKNOWLEDGEMENTS

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.

REFERENCES

REFERENCES
Brown
C. M.
Lund
J. R.
Cai
X.
Reed
P. M.
Zagona
E. A.
Ostfeld
A.
Hall
J.
Characklis
G. W.
Yu
W.
Brekke
L.
2015
The future of water resources systems analysis: toward a scientific framework for sustainable water management
.
Water Resour. Res.
51
,
6110
6124
.
Dell'Arciprete
D.
Vassena
C.
Baratelli
F.
Giudici
M.
Bersezio
R.
Felletti
F.
2014
Connectivity and single/dual domain transport models: tests on a point-bar/channel aquifer analogue
.
Hydrogeol. J.
22
,
761
778
.
DHI
2014
(
accessed 11 April 2015
).
ETAGI
1970
Intensiivse veevahetuse tsooni põhjavete miinimumäravool ja piesomeetriline režiim Eestis [Groundwater minimum discharge to the channel network from the active water exchange zone, and fluctuation of groundwater heads in Estonia]
.
Eesti NSV Teaduste Akadeemia Geoloogia Instituut (ETAGI)
,
Tallinn
,
Estonia
.
Fetter
C. W.
Jr
2014
Applied Hydrogeology
,
4th edn
.
Pearson New International Edition
,
Harlow
,
UK
.
GGPS
2015
Groundwater and global palaeoclimatic signals. International scientific project (GGPS) supported by UNESCO IGCP and INQUA TERPRO Committee. http://www.gw-gps.com/
(
accessed 11 January 2016
).
Gleeson
T.
Alley
W. M.
Allen
D. M.
Sophocleous
M. A.
Zhou
Y.
Taniguchi
M.
VanderSteen
J.
2012
Towards sustainable groundwater use: setting long-term goals, backcasting, and managing adaptively
.
Ground Water
50
,
19
26
.
Gorelick
S. M.
Zheng
C.
2015
Global change and the groundwater management challenge
.
Water Resour. Res.
51
,
3031
3051
.
Haravee
H.
2012
Võim ja raha sõitsid rahvusvaheliselt tunnustatud metsateadlase elutööst traktoriga mürinal üle [The might and money overrun the lifework of a forest scientist by caterpillars without any remorse]. Õhtuleht, 18 September 2012
.
Harbaugh
A. W.
2005
MODFLOW-2005: The U.S. Geological Survey Modular GroundWater Model – the Ground-WaterFlow Process. Techniques and Methods 6-A16
.
US Geological Survey
,
Reston, VA
,
USA
.
IGANE
1975
Harakteristika estestvennyh resursov podzemnyh vod Éstonii [Specification of the groundwater natural recharge in Estonia]
.
Institut geologii Akademii nauk Éstonskoj SSR (IGANE)
,
Tallin
,
Estonia
.
Kalm
V.
Raukas
A.
Rattas
M.
Lasberg
K.
2011
Pleistocene glaciations in Estonia
. In:
Developments in Quaternary Science
(
Ehlers
J.
Gibbard
P. L.
Hughes
P. D.
, eds).
15
,
Elsevier B. V.
,
Amsterdam
,
The Netherlands
, pp.
95
104
.
Karro
E.
Marandi
A.
2003
Mapping of potentially hazardous elements of the Cambrian-Vendian aquifer system, northern Estonia
.
Bull. Geol. Soc. Finland
75
,
17
27
.
Karro
E.
Marandi
A.
Vaikmäe
R.
2004
The origin of increased salinity in the Cambrian-Vendian aquifer system on the Kopli Peninsula, northern Estonia
.
Hydrogeol. J.
12
,
424
435
.
KM
2011
Jätkusuutlik põhjaveeseire süsteem Ida-Viru maakonnas. Norra finantsmehhanismi projekt 52/2006-EE0010. Projekti lõpukoolituse ja lõpuseminari materjalid [Sustainable groundwater monitoring system of East-Viru county, Estonia. Project 52/2006-EE0010 of the Norwegian Financial Mechanism. Proceedings of the summary workshop and seminar]. Keskkonnaministeerium (KM), Tartu Ülikool, Tallinna Tehnikaülikool, Eesti Geoloogiakeskus
,
Tartu
.
MAVES
2007
Tööstusjäätmete ja poolkoksi ladestuspaikade sulgemise ettevalmistus Kohtla-Järvel ja Kiviõlis. 2003/EE/16/P/PA/012. Keskkonnamõju hindamise aruanne [MAVES 2007 Arrangement to closing the depositing sites of industrial wastes and semi-coke in Kohtla-Järve and Kiviõli]
.
AS Maves (MAVES)
,
Tallinn
,
Estonia
.
MAVES
2008
Kirde-Eesti tööstuspiirkonna põhjavee orgaaniliste ühendite seire 2008. aastal [Monitoring of groundwater organic compounds in the industrial region of the northwest Estonia in 2008]
.
AS Maves (MAVES)
,
Tallinn
,
Estonia
. ).
MEE
2006
Harju sub-river basin district water management plan
.
Ministry of Environment of the Republic of Estonia (MEE), Grontmij Nederland bv
,
Tallinn
,
Amsterdam
.
Mulligan
K. B.
Brown
C.
Yang
Y.-C. E.
Ahlfeld
D. P.
2014
Assessing groundwater policy with coupled economic-groundwater hydrologic modelling
.
Water Resour. Res.
50
,
2257
2275
.
NGI
2004
Estonia, the oil shale industry
.
Risk based environmental site assessment of landfills. Norwegian Geotechnical Institute (NGI)
,
Oslo, Norway; Ministry of Environment of the Republic of Estonia
,
Tallinn
,
Estonia
.
Niitra
N.
2012
Teadlased teevad Kohtla-Järve tuhamäes sonkimise maatasa [Scientists are very critical about ongoing grabbing of the Kohtla-Järve landfill]. Postimees, 31 August 2012, nr. 6587
.
NRC
2012
Committee on the Challenges and Opportunities in the Hydrologic Sciences, Challenges and Opportunities in the Hydrologic Sciences
.
National Research Council
,
Washington, DC
,
USA
.
Perens
R.
Vallner
L.
1997
Water-bearing formation
. In:
Geology and Mineral Resources of Estonia
(
Raukas
A.
Teedumäe
A.
, eds).
Estonian Academy Publishers
,
Tallinn
,
Estonia
, pp.
137
145
.
Porman
A.
2015
Digital model of the shape of the Quaternary cover of the Estonian Artesian Basin
.
Master's Thesis
,
Department of Geography, University of Tartu
,
Tartu
,
Estonia
.
Protas'eva
M. S.
Eipre
T. F.
(eds)
1972
Resursy poverhnostnyh vod SSSR. T. 4, Pribaltijskij rajon. Estonia
[Surface water resources of Soviet Union. Vol. 4, Baltic Region, Estonia]. Gidrometeoizdat
,
Leningrad
.
Refsgaard
J. C.
Højberg
A. L.
Møller
I.
Hansen
M.
Søndergaard
V.
2010
Groundwater modeling in integrated water resources management – visions for 2020
.
Ground Water
48
,
633
648
.
Schlumberger
2015a
Visual MODFLOW Classic
.
Schlumberger Water Services
,
Kitchener, Ontario
,
Canada
.
Schlumberger
2015b
Visual MODFLOW Flex
.
Schlumberger Water Services
,
Kitchener, Ontario
,
Canada
.
Siegert
M. J.
2001
Ice Sheets and Late Quaternary Environmental Change
.
John Wiley & Sons
,
Chichester
,
UK
.
Strandberg
M.
2012
Põlevkivirämpsu põlengu varjamine meenutab Tšernobõli ja Fukushimat [Suppression of oil shale residuals burning recalls Černobyl and Fukushima]. Postimees, 28 July 2012
.
Vaikmäe
R.
Vallner
L.
Loosli
H. H.
Blaser
P. C.
Julliard-Tardent
M.
2001
Paleogroundwater of glacial origin in the Cambrian-Vendian aquifer of northern Estonia
. In:
Paleowaters in Coastal Europe: Evolution of Groundwater since the Late Pleistocene
(
Edmunds
W. M.
Milne
C. J.
, eds).
Geological Society, Special Publications
,
189
,
London
,
UK
, pp.
17
27
.
Vaikmäe
R.
Kaup
E.
Marandi
A.
Martma
T.
Raidla
V.
Vallner
L.
2008
The Cambrian-Vendian aquifer, Estonia
. In:
The Natural Baseline Quality of Groundwater
(
Edmunds
W. M.
Shand
P.
, eds).
Blackwell Publishing
,
Oxford
,
UK
, pp.
175
189
.
Vallner
L. K.
1976
Metodika i rezul'taty rasčeta podzemnogo pitanija rek Èstonii. V kn.: Trudy IV Vsesojuznogo gidrologičeskogo s’ezda. T. 8. Vzaimodejstvie poverhnostnyh i podzemnyh vod (A. A. Sokolov, red.) [Methodology and results of estimation of groundwater discharge in streams of Estonia. In: Proceedings of the 4th Hydrologic Congress of Soviet Union. Vol. 8. Relationship between surface water and groundwater (A. A. Sokolov, ed.)]. Gidrometeoizdat, Leningrad, SSSR, s
.
104
114
.
Vallner
L.
1994
Problems of groundwater quality management in Estonia
. In:
Groundwater Quality Management
(
Kovar
K.
Soveri
J.
, eds).
IAHS Publ. No. 220
,
Wallingford
,
UK
, pp.
449
460
.
Vallner
L.
1996
Groundwater
. In:
Estonian Environment. Past, Present and Future
(
Raukas
A.
, ed.).
Ministry of the Environment of Estonia
,
Tallinn
,
Estonia
, pp.
60
71
.
Vallner
L.
1997
Groundwater flow
. In:
Geology and Mineral Resources of Estonia
(
Raukas
A.
Teedumäe
A.
, eds).
Estonian Academy Publishers
,
Tallinn
,
Estonia
, pp.
137
152
.
Vallner
L.
2003
Hydrogeological model of Estonia and its applications
.
Proc. Estonian Acad. Sci. Geol.
52
,
179
192
.
Vallner
L.
Järvet
A.
1998
Groundwater recharge and extraction in Estonia
. In:
XX Nordic Hydrological Conference
(
Kajander
J.
, ed.).
The Nordic Coordinating Committee for Hydrology
,
Helsinki
,
Finland
, pp.
283
293
.
Vallner
L.
Ivask
J.
Marandi
A.
Vaikmäe
R.
2015b
Simulation of 18O concentration in the Estonian Artesian Basin
. In:
4th Annual Meeting of G@GPS IGCP 618 Project. Paleogroundwater from Past and Present Glaciated Areas. Estonia, 5–9 July 2015
(
Pärn
J.
Raidla
V.
Vaikmäe
R.
Raukas
A.
Bauert
H.
, eds).
Tallinn University of Technology
,
Tallinn
,
Estonia
, pp.
32
34
.
Virbulis
J.
Bethers
U.
Saks
T.
Sennikovs
J.
Timuhins
A.
2013
Hydrogeological model of the Baltic Artesian Basin
.
Hydrogeol. J.
21
,
845
862
.
Water Directive
2000
Directive 2000/60/EC of the European Parliament and of the Council establishing a framework for the Community action in the field of water policy (Water Directive 2000). http://ec.europa.eu/environment/water/water-framework/index_en.html
.
Weixing
G.
Langevin
C. D.
2002
User’s Guide to SEAWAT: A Computer Program for Simulation of Three-Dimensional Variable-density Ground-water Flow
.
Techniques of water-resources investications 6-A7. US Geological Survey
,
Tallahassee, FL
,
USA
.
Zheng
C.
2010
MT3DMS v5.3: a modular three-dimensional multispecies transport model for simulation of advection, dispersion and chemical reactions of contaminants in groundwater systems. Supplemental User's Guide. Technical Report
.
The University of Alabama
,
Tuscaloosa, AL
,
USA
.
Zheng
C.
Bennett
G. D.
2002
Applied Contaminant Transport Modeling
,
2nd edn
.
John Wiley & Sons
,
New York
,
USA
.
Zheng
C.
Wang
P. P.
1999
MT3DMS: A modular three-dimensional multispecies model for simulation of advection, dispersion and chemical reactions of contaminants in groundwater systems; Documentation and User's Guide, Contract Report SERDP-99-1
.
US Army Engineer Research and Development Center
,
Vicksburg, MS
,
USA
.