Far-field transport modelling for a repository in the Boom Clay in the Netherlands

Abstract

A groundwater model was set up to study far-field transport for the potential of a radioactive waste repository the Boom Clay in the Netherlands. The existing national groundwater model, the Netherlands Hydrological Instrument, was extended in the vertical direction to include geological formation up to and beyond the Boom Clay. As the amount of hydrogeological data in the deeper subsurface is limited, simplifications in the model schematisation were necessary. Moreover, nationwide data about the tops and bottoms of many of the deeper geological formations and their members are lacking and required interpolation. Finally, values for hydrogeological parameters, such as porosity and hydraulic conductivity, are also lacking for the deeper formations. These values were estimated using relationships with depth and lithology. Moreover, no quantitative data about heterogeneity within the deeper geological formations or its members were available.

In the Dutch research programme on the geological disposal of radioactive waste (OPERA), the post-closure safety of a generic repository is assessed in either Boom Clay or rock salt. Disposal of Dutch radioactive waste is not foreseen in the next decades and a choice of host rock has not been made. In the early, conceptual phase of the radioactive waste disposal process in the Netherlands no potential repository locations were selected and a groundwater flow model for the entire Netherlands was build. As a starting point a geological disposal facility is assumed to be present at a depth of at least 500 m within a Boom Clay formation of 100 m in order to be able to make an assessment of post-closure safety with this geological formation in a disposal concept. With these assumptions, a general idea of potential flow patterns has been obtained and broken down into pathline trajectories. These trajectories were calculated to achieve input for the potential transport of radioactive isotopes (radionuclides) from this waste in the Netherlands after the closure of a disposal facility in Boom Clay.

The groundwater flow patterns in the deeper subsurface strongly resemble the larger scale flow patterns in the shallow subsurface, with flow from infiltration areas in the east and the south of the Netherlands towards to seepage areas of the polders in the west and the northern part of the country or towards the river valleys of the Rhine and IJssel. Groundwater flow velocities, however, are much lower in the deeper part of the model and consequently travel times are much larger. The conservative travel times from the pathlines range from a few 1000 years to more than 10,000,000 years depending on the location for the repository. Longer travel times are obtained for locations with a downward groundwater flow in the Boom Clay.

Because of the simplifications in the model schematisation and the uncertainty in the model parameters, the present results should only be considered as a first indication. Moreover, the model could not be validated due to a lack of validation data. However, the insight gained with the model may help to design a data collection strategy for dedicated model validation, such as measuring the hydraulic gradient over the Boom Clay to validate downward flow in the Boom Clay to obtain the necessary data for a post-closure safety assessment.

Introduction

The OPERA research programme studies the geological disposal of radioactive waste in the Netherlands. In this research program, the post-closure safety of a generic repository in the Boom Clay is investigated. A number of barriers are considered, starting from the engineered ones and the potential host rock, the Boom Clay itself. Knowledge concerning the transport through the formations surrounding the Boom Clay is essential to calculate how the radionuclides can enter the biosphere. This paper describes the hydraulic transport in the far-field, starting at the interface of the potential host rock, towards the biosphere, and translates the hydraulic transport into a simplified transport model that will be used for the radionuclide migration.

Disposal of Dutch radioactive waste is not planned in the coming decades. In the early, conceptual phase of the radioactive waste disposal process in the Netherlands, the geological disposal facility was assumed to be present at a depth of 500 m within a Boom Clay formation of 100 m in order to be able to make an assessment of post-closure safety (Verhoef et al., 2011). Because a generic repository was investigated, a groundwater model for the entire Netherlands was set up. This contrasts with many groundwater flow studies for disposal of nuclear waste, such as the Forsmark site in Sweden (Hartley & Joyce, 2013), a site in northeast Belgium (Gedeon et al., 2007), the Olkiluoto site in Finland (Blessent et al., 2011), the Sellafield site in the UK (Blum et al., 2005) and Yucca Mountain in the USA (Bodvarsson et al., 1999), where site-specific studies could be performed.

Moreover, the amount of data on which the deepest part of our model is based is limited. Much of the model data on the national scale has been obtained using interpolation, extrapolation and generic relations between certain properties, such as porosity and hydraulic conductivity, with depth and lithology.

Methods and materials

As the potential pathline trajectories cover a large part of the country, a groundwater flow model of the entire country was used. Within the Netherlands a model of the entire country is available as the Netherlands Hydrological Instrument (NHI; De Lange et al., 2014), which is reported in detail at http://www.nhi.nu. This model was extended in the vertical direction to include the Boom Clay and its adjacent geological formations. Moreover, the groundwater model was made steady state to obtain average groundwater flow velocities, which were used in the computation of the pathline trajectories.

Model schematisation

The original NHI is a transient model using time-varied rainfall, evaporation and surface water data that combines saturated groundwater with unsaturated groundwater and with a surface water distribution model. For the present study, only long-term, time-averaged flow conditions in the saturated zone are of interest, therefore the NHI is simplified to a steady-state model for saturated groundwater only and groundwater recharge data were obtained by time-averaged fluxes from the original NHI.

The NHI consists of seven layers and in a large part of the country it does not reach as deep as the Boom Clay, but the lowest model layer belongs mostly to (parts of) the Maassluis, Oosterhout or Breda Formations in different parts of the country. The depth of the original NHI is shown in Fig. 1. The model was extended in the vertical direction with 11 additional model layers to the bottom of the Lower North Sea Group. The vertical model extension is several hundreds of metres in the larger part of the country and up to 1500 m in the Ruhr Valley Graben. Table 1 gives an overview of the various groups, formations and members starting from the Maassluis Formation. The table also shows if this layer is completely, party or not at all included in the original NHI and the source of a nationwide grid data, when available. It also indicates if and in which model layer the geological layer is included in the extended model. In some parts of the country, parts of or even the entire formations of Maassluis, Oosterhout and Breda are not included in the original NHI. The missing parts of these formations are included in model layers 8–10, respectively. Where these formations are completely included in the original NHI, the thickness of the corresponding model layers 8–10 is zero. A typical cross-section is shown in Fig. 2.

Fig. 1. Depths of the base of the original NHI model.

Table 1. Schematic overview of the geological formations and their members starting at the Maassluis Formation, their presence in the NHI, sources of grid data about the base and the layer number in the model schematisation.

+, completely included in the original NHI; ±, partly included in the original NHI; –, not at all included in the original NHI.

Fig. 2. Geological cross-section.

It was decided not to include several geological layers, such as the Ville, Inden and Tongeren Formations, as nationwide data because their tops and bottoms were missing.

The lowest model layer includes the entire Lower North Sea Group due to a lack of nationwide grid data of the formations and members within this group.

As the Lower North Sea Group contains some widespread thick clay deposits, such as the Ieper and Dongen Clay Members, it was expected that groundwater fluxes that go through the base of the Boom Clay would predominantly flow horizontally through the underlying Vessem Member and the formations of the Lower North Sea Group, where the permeability is higher. A small part of that flux through the Boom Clay may reach formations below the Lower North Sea Group, but the travel times of these pathline trajectories are expected to be higher and thus less critical. The omission of the geological formations below the Lower North Sea Group in our model is expected to have less influence compared to errors introduced by the present model schematisation and uncertainty in the parameter values of the formations that were incorporated in the model.

The model layers directly adjacent to the Boom Clay, such as the Someren, Veldhoven Clay, Voort, Steensel and Vessem Members, were analysed in more detail as it was assumed that they would have more influence on the model results. While analysing the Boom Clay characteristics, Vis and Verweij (2014) also interpreted the underlying and overlying Formation Members, and their results are shown in Fig. 3. The left-hand figure shows that the Steensel member is only present in the orange areas. The above-lying Voort Member is present in the plain yellow area and possibly also in the area where the Steensel Member is present. Together with some specific information in the DINOloket (see https://www.dinoloket.nl/tertiary-and-quaternary and linked web pages), for example that the Voort member rapidly thins in the northern direction as well as limited borehole data, fractions of the Someren, Veldhoven Clay, Voort and Steensel Members on a number of locations were obtained. By linear interpolation using these fractions and using the entire thickness between the top of the Boom Clay and the bottom of the Breda Formation, the thicknesses of all the members between the Breda Formation and the Boom Clay were estimated. The thicknesses of these formation members range between 0 and 216 m (Someren), 0 and 288 m (Veldhoven Clay), 0 and 310 m (Voort) and 0 and 30 m (Steensel). Finally, the Vessem Member was classified in thicknesses of 25 or 50 m based on the yellow and orange regions in the right-hand part of Fig. 3.

Fig. 3. Information on the geological layers overlying (left) and underlying (right) the Boom Clay (figures from Vis & Verweij, 2014).

It is known that there is heterogeneity within the formations and its members, but again nationwide data about it is lacking. In the present model, these formations were modelled as homogeneous layers but with different values for the vertical and hydraulic conductivity in which the up-scaling effect of vertical layering is incorporated. Preferential flow through more permeable deposits within a formation or its member is likely.

Hydraulic conductivity and porosity

For the geological layers that are not included in the database REGIS, no direct grid data of hydraulic conductivity is available. For these formations, a methodology was applied that is used by TNO for modelling oil and gas reservoirs. It makes use of tables which describe (1) the lithology per group, formation or member, (2) a relation between lithology, depth and porosity, (3) a relation between porosity, lithology and horizontal and vertical permeability, and (4) the relation between the viscosity and depths using increasing temperature to relate the permeability to the hydraulic conductivity. The methodology is shown in the supplementary material for this paper (available online at http://dx.doi.org/10.1017/njg.2016.13). The hydraulic conductivity in both the horizontal and vertical directions of the Boom Clay is shown in Fig. 4. Hydraulic conductivity measurements for the Boom Clay show comparable values for samples that have been categorised as non-mud samples. Mud samples have hydraulic conductivity values two orders of magnitude lower. These values have also been used in scenario calculations but these results are not shown in this paper.

Fig. 4. Horizontal (left) and vertical (right) hydraulic conductivity of the Boom Clay (in m/day).

Boundary conditions

Detailed information about the hydraulic heads near the borders in these deeper layers is lacking, therefore constant head boundary conditions were applied in the extended model layers in Belgium and Germany. The leakage factor, which is the square root of the product of the transmissivity of an aquifer with the vertical resistance of the layers above, for the deeper layers is a few kilometres at most. The effect of fixed boundary head values on the modelled hydraulic heads will therefore become negligible within a distance of 10 km from the border. The fixed hydraulic heads at the deeper layers are set equal to the lowest NHI model layer, as these values are expected to be the best estimate for the heads in the deeper layers considering the limited amount of information. At the North Sea, no flow boundary conditions were imposed for the extended model layers, which is equal to the boundary condition of the lowest model layer in the NHI. For the deeper model layers, this choice for a boundary condition neglects inward and outward fluxes at the lateral boundary but forces all pathlines to go to the surface and the calculated travel times are equal to the arrival times at the biosphere.

Transport

The conservative transport was calculated using pathlines. Pathlines were started at the top and the bottom of the Boom Clay on 4618 locations where Boom Clay is expected to be found at a depth of at least 500 m with a thickness of at least 100 m from currently available knowledge. For each location, the particle with the shortest travel time is considered to be the most critical pathline. Using the most critical pathline starting at all locations, a distribution of travel times was obtained. The 10, 50 and 90 percentiles of this distribution are used to obtain three pathlines that are the bases of a one-dimensional radionuclide transport model that is not presented in this paper.

Neglected processes

Besides the assumptions needed due to the lack of nationwide data, there are also a few processes that were not taken into account. These processes include (1) the effect of density variations on groundwater flow due to variations in the salt distribution and temperature, (2) the effect of faults on the groundwater flow and (3) the processes of dispersion and diffusion.

The effect of density variations due to the salt distribution and temperature variation are neglected due to a lack of data on the salt distribution and the computational demands if the transient behaviour of the salt distribution were taken into account: this would require a complete salt transport model for a simulation period of more than 100,000 years. The effect of faults is also neglected due to a lack of data (the hydrological impact of faults quantitatively) and computational demands as it would require more model layers to create direct flow connections between different geological units.

Changes in the physical system are also likely to take place during the time scale of 100,000 years. Geological scenarios are planned to be incorporated as well as some changes in the human interference with the physical system. These scenarios have not yet been defined.

Analyses and results

The groundwater flow and pathline were calculated using the software iMOD (Vermeulen et al., 2013), which is based on the model codes MODFLOW (McDonald & Harbaugh, 1998) and MODPATH (Pollock, 1994). Because of all the assumptions on the model parameters and model simplification, the results are only considered to be a first estimate of the pathline trajectories and their travel times.

Hydraulic heads

Fig. 5 shows the hydraulic head distribution in model layer 16 (Vessem Member). This is the first model layer underneath the Boom Clay. This figure gives information about the flow pattern in the deeper part of the model, for instance at the Veluwe area higher hydraulic heads are present in comparison with its surroundings. It can be concluded that flow through the Boom Clay is large enough to result in a clear rise in hydraulic head levels at the Vessem Formation. The large-scale pattern of infiltration and seepage areas that is observed in the NHI model is also present in the deeper model layers. An important part of the flow resistance for the deeper groundwater flow system is due to the low horizontal permeability in the aquifers below the Boom Clay and its long flow distance.

Fig. 5. Hydraulic head distribution in the Vessem Member below the Boom Clay (in m).

Pathlines

Pathlines were obtained by starting particles at the top and bottom of the Boom Clay on locations where Boom Clay is expected to be found at a depth of at least 500 m with a thickness of at least 100 m from currently available knowledge. The flow pattern for these pathlines is from infiltration areas, such as in the south and east of the Netherlands, towards the polder areas in the western and northern part of the country and in some areas towards the main rivers of the Rhine and IJssel.

For a single pathline, the cross-sectional view is shown in Fig. 6. This cross-section follows the curved trajectory of the pathline. The colours of the pathline indicate the conservative travel times starting from the interface with the Boom Clay. This pathline starts at the bottom of the Boom Clay and goes into model layer of the Lower North Sea Group. There, it travels a significant horizontal distance before it flows almost vertically upwards through the Boom Clay and travels a significant horizontal distance through the model layers of the Veldhoven Clay, and Breda and Oosterhout Formations. The travel distance and the residence time in the part of the original NHI model are relatively small in comparison with the deeper model layers. The details of this pathline are also shown in Table 2 in the columns for the slow pathline. The phenomenon of the pathline staying most of its travel distance and residence time within the model layers below the original NHI model is typical for all pathlines.

Fig. 6. Cross-sectional view for a single pathline.

Table 2. Residence times and travel distance per geological layer for three selected pathlines.

Travel times

For each pair of pathlines starting at the top and bottom of the Boom Clay on the same location, the pathline with the shortest travel time was considered the critical pathline for that location. No pathline left the model domain at the lateral boundaries. The travel time distribution of all critical pathlines is shown in Fig. 7. The distribution shows a large variation, with travel times from a few thousand years to more than 10 million years and a majority of travel times over 100,000 years. It is clear that with the present data the range in travel times for a generic – that is not site-specific – geological disposal facility is large. In general, the longer travel times coincide with areas where the flux through the Boom Clay is downward and the most crucial particle is started at the bottom of the Boom Clay. At some location downstream, this particle has to break through the Boom Clay and both the residence time there and in the Vessem and/or Lower North Sea Group model layers result in relatively long travel times.

Fig. 7. Travel time distributions of the critical pathlines.

Finally, for three selected pathlines the residence time and the travel distance in each model layer is summarised in Table 2.

Discussion and conclusions

Discussion

The present NHI is still insufficient to provide input for groundwater flow in cases of geological disposal of radioactive waste in the Netherlands. A first attempt was made to extend this instrument for geological disposal of the waste. The extension presented in this paper is based on numerous assumptions, simplifications and indirectly derived parameter estimates. The main model simplifications that are expected to have a significant effect on the model results are the disregard for density effects, faults and the heterogeneity of hydrological parameters within the geological layers and within the model layers that consisted of more geological formations or members due to the lack of nationwide data. Our expectation would be that adding the effect of present density differences and adapting them dynamically would decrease the deeper water fluxes in cases of deep saline water. Any displacement of this more saline water towards the surface would counteract the driving force of the deeper groundwater flow. When adding the effect of salt dissolution from diapirs, our expectation would be that it can, in different locations, both counteract and add to the present hydraulic forcing in the deeper groundwater system, therefore temporally part of the water from a deeper origin would reach the surface more quickly and some additional mixing will take place. With respect to faults, our expectation is that incorporation of faults is likely to decrease the calculated travel times but as long as these faults do not have high vertical hydraulic conductivities and provide a connection over a larger vertical distance, the impact will not be too severe as in the present model a significant part of the flow resistance is due to the horizontal hydraulic conductivity in combination with the large horizontal scale of the deeper flow systems. With respect to the heterogeneity, a distinction should be made between the heterogeneity at large and relatively small scales. Heterogeneity at the large scale is due to different long-term geological conditions, such as the depth of the sea, which resulted in deposition of finer sediment in the present northern part of the Netherlands in comparison with the nearer shore parts of the present southern Netherlands. Our expectation is that the characterisation method that was used overestimates hydraulic conductivity in the northern part of the country and therefore calculated travel times there are being underestimated. Heterogeneity at a more local scale will be present, as a single Formation Member often consists of sequential layers of different lithologies with variable horizontal extensions. Our expectation is that for such heterogeneities the travel times can be significantly smaller if the layers with higher hydraulic conductivities are connected over long distances and thus create shortcuts in both horizontal and vertical directions. More information about this heterogeneity and the effect of faults is needed to improve the model.

The present model was not validated due to the lack of good validation data for transport under the assumed steady-state groundwater flow conditions. A limited amount of isotope data with groundwater ages of several tens of thousands of years are available (Glasbergen, 1987). However, the groundwater flow system was considerably different during large parts of this period due to meteorological (sea level, climate) and anthropogenic (decreasing artificial surface water levels in the polders areas due to land subsidence and groundwater abstractions) influences. Using the isotopic groundwater age data for model validation would require a model that incorporates these changes in the hydrologic system.

Nevertheless the model gives important insight. It shows that the groundwater flow pattern in the deeper part of the model domain resembles flow patterns from the shallow domain, i.e. flow from infiltration areas in the south and eastern parts of the country towards seepage areas, such as the polders in the western and northern parts of the country and the rover valleys of the Rhine and IJssel. In contrast, the groundwater flow velocities are much lower and travel times are much higher. For the majority of other repository locations where Boom Clay is expected to be present at a depth of at least 500 m with a thickness of at least 100 m for currently available knowledge, travel times would exceed 100,000 years. Longer travel times coincide in general with locations with a downward groundwater flow in the Boom Clay.

Finally, the effect of geological scenarios should be taken into account. These scenarios are presently being constructed within the OPERA and their effect on the hydrogeology will be incorporated once they become available. In reality, there will be a gradual change from the present situation to a new geological scenario on the time scale of the calculated travel times. Using travel times based on steady-state flow paths does not exactly mimic reality. In our opinion the contribution of the safety function of the far-field can also be evaluated when steady-state travel time calculations for the hydrogeological scenarios are added. Uncertainty about the locations of geological features such as rivers, tunnel valleys, ice coverage extensions, etc., should be incorporated.

Conclusions

It has been possible to construct a groundwater flow model for the entire Netherlands that is extended from the NHI vertically to include the Boom Clay and some geological formations below in order to contribute to the post-closure safety of a geological disposal facility. The extension is based on information that was readily available. Data on heterogeneity with the geological formations and the quantification of the effect of faults are expected to improve the model. The effect of geological scenarios still needs to be taken into account in order to use the model results in a post-closure safety assessment. Applying an uncertainty or sensitivity analysis to the model could provide input for a data collection strategy.

The flow pattern in the deeper part of the model resembles the large-scale flow pattern in the shallow groundwater. It clearly shows groundwater flow from the infiltration areas in the east and south of the country towards the seepage areas of the polders in the western and northern parts of the country and the river valleys of the Rijn and IJssel. The flow velocities in the deeper part of the model are much lower, however, mainly due to the low horizontal hydraulic conductivity in the deeper aquifers.

Using the present model, conservative travel times from the interface of the Boom Clay towards the biosphere can be large. Starting at locations where Boom Clay is expected to be found at a depth of at least 500 m with a thickness of at least 100 m, travel times range from a few thousand years to more than 10,000,000 years. The majority of pathlines starting at these locations have a conservative travel time of more than 100,000 years.