Hydrogeological modelling to support urban planning in harbour areas: a case study from Horsens, Denmark

Historically, industries were in harbour areas of cities for easy access to transportation of resources. Today, transforming former industrial areas into living spaces has become an attractive business. However, this transformation has often been challenged by high levels of soil contamination caused by the industrial use. Remediation measures are mandatory to ensure the public safety in the redeveloped areas. Detailed information about the contaminant type, distribution, and transport mechanisms is required to address the contamination issues. This paper suggests a workflow for investigations assisting decision making for construction work in redeveloped industrial areas. The workflow is applied to Horsens harbour (Denmark). In this area, renovation of the harbour walls introduces the risk of spreading of contamination to planned construction areas. The study demonstrates how detailed hydrogeological information about the site allows for scenario modelling of contaminant transport, guiding remediation efforts and aiding decision makers in developing the harbour area.


Introduction
Groundwater in urban areas is exposed to anthropogenic influences and is often contaminated (Vasin et al. 2016). Various sources of contamination can exist and challenge remediation efforts and future redevelopment due to changing land use (Schirmer et al. 2013). Therefore, the practical contaminated site investigations and measures to determine the transport of contaminants in urban areas are necessary and have been scientifically scrutinised for several years (Bauer et al. 2004;Greis et al. 2012;Trowsdale and Lerner 2007;Vasin et al. 2016). However, insufficient information about the hydrogeological setting and the numerical transport model have made the evaluation of contaminant concentration and distribution significantly uncertain (Hansen et al. 2019;Yang et al. 1999). Consequently, the effects of collection and remediation schemes developed seldom meet expectations (Graber et al. 2008;Khan et al. 2004).
Knowledge of the hydrogeological setting in urban areas is crucial because it forms the basis for understanding the fate and distribution of the contaminants. In recent years, geophysical methods have become a valuable, non-invasive source of continuous data sections, significantly improving the data coverage and accuracy of the investigated urban areas (Bockhorn et al. 2015;Pazzi et al. 2016;Prudhomme et al. 2019). Geophysical methods have been applied at contaminated sites to identify fill materials and, in some cases, map contamination plumes (Boudreault et al. 2010;Maurya et al. 2017;Vaudelet et al. 2011). However, the anthropogenic structures in urban areas represent abrupt changes which sometimes shortcut the geological units. Hence, the urban hydrogeological models often suffer from low fidelity which reduces the accuracy of numerical simulations.
Consequently, it is important to include those features in the urban models (Andersen et al. 2018).
Mapping contaminant plumes provides a better understanding of contaminated sites which enables sound risk assessment and targeted remediation schemes (Tao et al. 2019). The first step to delineate the geometry of contaminant plumes is usually point sampling. Recently, the mapping approach was developed from a simple interpolation to novel methods including stable isotopes tracing, geophysical mapping, geostatistics, and hypothetical numerical simulation (Baawain et al. 2018;Karges et al. 2018;McLean et al. 2019;Negrel et al. 2017;Pannecoucke et al. 2020;Rivett and Allen-King 2003;Zuo et al. 2018). Once the outline of contaminants has been revealed, the treatment technology takes over. Many in situ and ex situ approaches, including permeable reactive barriers, capture wells, and electrochemical transformation, have been tested and proven effective (Bertrand et al. 2016;Christ and Goltz 2004;Hyldegaard et al. 2020;Vaezihir et al. 2020). In some cases, the contaminated areas are too large or expensive to consider remediation on a shorter time scale. If construction work is planned in such an area, it is beneficial to model the resulting changes in the expected flow pattern. This may mitigate the risk of spreading contamination.
Recently, 3D numerical flow and solute transport software and code, including FEFLOW (Finite Element subsurface FLOW system), Visual MODFLOW, and CSTREAM, have attracted considerable interest due to their excellent functionality and flexibility (Bauer et al. 2004;Jeppesen et al. 2011;Trowsdale and Lerner 2007). FEFLOW has been successfully applied to simulate groundwater recharge and denitrification below a seasonal flooded restored riparian zone, as well as for constructing a freshwater management system and assessing seawater intrusion in coastal aquifers (Elad et al. 2017;Jensen et al. 2017;Kazmierczak et al. 2016;Nocchi and Salleolini, 2013;Priyanka et al. 2018;Rona et al. 2014;Sherif et al. 2012). In urban areas, FEFLOW also acts as a reference for groundwater resource management and underground space engineering (Du et al. 2016;Greis et al. 2012;Vazquez-Sune et al. 2005). In recent years, the solute transport simulation has played an increasingly important role in forming remediation schemes for groundwater (Bauer et al. 2004;Mulligan and Yong 2004;Vasin et al. 2016). This paper evaluates the resultant shape of a contamination plume simulated for different possible scenarios suggested in the context of renovation and construction works planned in the Horsens harbour area. The example demonstrates how the simulation may help decision making.
Horsens harbour (Fig. 1) was formerly an industrial area from 1860 to 1969 consisting of gasworks, an asphalt roofing factory and a tar production facility. It was found that the former gasworks site was characterised by multiple contaminant sources. The contaminants of concern had the potential to continuously feed the contaminant plumes (Birak and Miller 2009;Bockelmann et al. 2001;Mak et al. 2006). Today, the upper metres of soil are still heavily contaminated even though the area was remediated in the 1980s and 1990s. Nonetheless, converting the industrial port into living spaces has become an attractive business because people willingly pay higher rent to live near both the city centre and the sea.
This study suggests a workflow for modelling of contaminant transport in areas where construction works alter the hydrological situation. It demonstrates how different scenarios can visualise the variation in the resulting contamination plume and assist decision making in complex contaminated sites when detailed hydrogeological information is available. The study focuses on integrating the field and numerical methods such as geophysical mapping, 3D hydrogeological model construction, and solute transport simulation in an urban area. In the demonstration case, the purpose is to help plan the renovation at a harbour front to guide water to the sea while considering various interests. Andersen et al.'s work (2018) presents a detailed 3D geological model of the former gasworks at Gasvaerksgrunden in Horsens, Denmark. The previous work was based on joint interpretation of various geophysical methods, boreholes, water samples, and land use information. Building on these premises, the research objectives of this paper are: 1. To import the 3D hydrogeological voxel model to FEFLOW by assigning proper hydrogeological parameters for it, 2. To explore the transport of the contaminant plume according to the identified contamination source information, 3. To evaluate the changes in the contaminant plume for different scenarios related to the renovation of the sheet pile walls, 4. To summarise, providing a good foundation to the authorities for their current and future work schemes related to site construction and groundwater remediation in the contaminated harbour aquifer.

Suggested workflow
This section suggests a workflow (Fig. 2) which has been developed and tested as part of this study. The field geophysical mapping, and lithological and anthropogenic data interpretation for the 3D geological model construction in GeoScene3D have been detailed by Andersen et al. (2018).
Based on their research, this paper focuses on demonstration of the last part of the workflow: the hydrological parametrisation of the 3D geological model and simulation of different scenarios in FEFLOW. The recommended workflow presented below has been developed based on the combined experience of the two studies. This section will go through the suggested workflow and clarify its advantages. Initially, all relevant data from the site are acquired. With regard to the suggested workflow, studying the fate of contaminants requires additional data. This includes geophysical data obtained from multiple high-density surveys and distribution of anthropogenic structures among others including pipes, basements, and information about former excavations and fill materials. Once the data coverage is reasonable, considering the expected geological complexity, a 3D detailed voxel model representing lithology is created.
For information on data collection and setup of the model for the area of Horsens harbour, please refer to Andersen et al. (2018). Once the hydrogeological model is parametrised and imported to FEFLOW, simulations are run and the model is calibrated to fit groundwater levels observed in the study area. When the hydrogeological model is calibrated, it is possible to build up the fluid flow and solute transport model. For the case of Horsens harbour, the sheet pile walls are adjusted in the simulations according to the suggested scenarios. Decisions can be made based on the observed differences for the chosen scenarios.
The suggested workflow has its strength in the following aspects. A traditional estimate based on groundwater samples and a simple ground model is unlikely to resolve the complex structure of urban areas originating from multiple construction works through time. The suggested high-density geophysical surveys can provide non-invasive investigations and improve the characterisation of the subsurface, which often shows a large lateral variation originating from previous construction works (Andersen et al. 2018). Based on the geophysical data, boreholes may be planned for calibration of the geophysical data. Anthropogenic features significantly alter the subsurface, potentially perturbing groundwater flow and contaminant transport especially in harbour areas where groundwater is often near the terrain. Consequently, this information should be incorporated in the model (Andersen et al. 2018). The detailed voxel model allows for this adding of anthropogenic features like pipes and basements, which cannot be represented in a meaningful way in a conventional layered model. Furthermore, in the suggested voxel model, it is numerically easier to represent pinching out layers and lenses. Scenario test and prediction based on the hydrogeological modelling and simulation can reveal changes in the geometry of the contaminant plumes derived from planned changes in an area. In this way, potential problems may be identified and mitigated prior to initiating construction works at a site. In this study, the suggested workflow is applied to evaluate different scenarios and their impact on contaminant transport in the research area. As a result, the optimal solution is selected and the results are presented to the decision makers in Horsens municipality for conclusion and application.

Site description
The study is conducted on the Horsens Gasworks site, Denmark (Fig. 1). The site is designated as one of the nine national test sites for the development and demonstration of new techniques related to the investigation and prevention of soil contamination (https:// www. danis hsoil. org/ tests ites/ index. php? lang= uk). Various urban features including roads, houses, parking lots, and green areas constitute this area. An estimated 45,000 m 3 of soil is still heavily contaminated in spite of the Fig. 2 Map presenting the suggested workflow remediation efforts in the 1980s and 1990s. Today, 5 up to 7 m of the upper soil is contaminated with volatile aromatics, phenols, polycyclic aromatic hydrocarbons (PAH), and inorganic compounds. The degradation rate of the contaminants is expected to be extremely slow due to the anaerobic conditions which hinder the reproduction of aerobic microorganisms that might degrade the contaminants. A contamination plume has been detected along the western part of the area that extends towards the harbour. It is believed to originate from the gasworks vicinity field.
The land and sea areas are originally separated by two sheet pile walls (named Wall 1 and Wall 2, respectively; see Fig. 1). Wall 2 is a concrete structure and considered impermeable, while Wall 1 is an old wooden structure which is permeable. Due to the ongoing degradation, Wall 1 is planned to be reinforced by adding an impermeable steel structure in front of the existing wooden wall. This construction operation is to avoid the sliding of soils behind the existing wall. This modification is expected to change the hydrological conditions and, consequently, influence the transport pattern of the contaminants. Simultaneously, several construction sites including commercial and residential projects have been scheduled in the Horsens harbour area. An impermeable Wall 1 is expected to widen the contamination plume and increase the concentration of contaminants in the planned construction sites. This should be avoided for several reasons: Increased concentrations of pollutants in groundwater will increase construction costs if lowering of the groundwater table is needed. Additional costs would be associated with handling of a larger volume of contaminated soil and if volatile pollutants migrate from soil to indoor air then installation of membranes and increased ventilation will be needed. Thus, different solutions need to be tested.

Geology
The geology at Horsens harbour is complex, with sediments ranging from glacial to post-glacial and anthropogenic sediments. The glacial deposits include moraine clay, thin layers of meltwater deposits overlain with lacustrine clays and silt. The glacial sediments were deposited during the Pleistocene (Houmark-Nielsen 2004). The postglacial sediments include several thin layers of sand, clay, and organic clay (Gyttja) deposited during the Holocene when the area experienced several regressions/transgressions (Houmark-Nielsen 2004). The glacial and postglacial deposits are overlain with fillings that display abrupt lateral variations and thicknesses up to 5 m. The fillings around anthropogenic structures, for e.g. pipelines, sewers, and foundations, are predominantly sand and gravel.

Borehole data
The lithological information from 98 boreholes at the field site was collected from the national borehole database, Jupiter, hosted by The Geological Survey of Denmark and Greenland (GEUS), various geotechnical investigations at the field site, and the Central Region of Denmark's database (Fig. 3). The borehole depth varies from 0.5 to 28 m below terrain (m.b.t.) with most boreholes being drilled to around 7 m.b.t. In addition to the lithological information for several boreholes, the databases provide information about the groundwater level and chemical analyses. The quality of the borehole information varies due to differing borehole ages, drilling purposes, and drilling methods.

Groundwater level
Geophysical surveys and borehole information suggested the presence of two sandy aquifers in the study area. The upper aquifer is mainly distributed in the central and northern parts of the model, whereas the lower aquifer is thinner but generally continuously spread over the study area. The groundwater levels in the aquifers were derived from the available borehole information (section "Borehole data"). Data was filtered to obtain the authentic water level representing the continuous groundwater flow field in the aquifers rather than stagnant water lenses, discontinuous or deeper water levels. Hence, water level data from boreholes were selected only if the screen, or the bottom, in case of open boreholes, was in the upper or lower sandy aquifer. Additionally, outliers, where the hydraulic head listed was above the terrain or far deeper than most of the measured values, were discarded. Finally, 52 qualified boreholes were selected based on the above conditions (Fig. 5). In some boreholes, series of measures were available at different times. In such cases, the median value was adopted. Still, variation in water levels sustained, most likely due to seasonal changes, tides, and rainfall differences since measurements were not conducted at the same time. Ideally, a synchronous series of measurements should be undertaken. This was not possible here since the boreholes were temporary and often originated from geotechnical investigations before construction works. The upper aquifer has a hydrostatic pressure around 3.383 to 3.708 m above sea level. The hydraulic head in the lower aquifer varies from 0.8 m below sea level to 2.3 m above sea level. The groundwater elevation generally decreases from north to south.

Analyses
Analyses of water samples from 29 boreholes, soil samples from 12 boreholes, and pore air from 5 boreholes are available from the study area (Fig. 6). The samples were analysed for various components among other volatile aromatics. For e.g. benzene, phenols, PAH, naphthalene, and inorganic components such as cyanide, ammonia, and sulphate. Furthermore, NVOC (non-volatile organic carbon), electrical conductivity, pH, redox potentials, and temperature. This article aims to investigate the transport of contaminants from the gas works stations. Phenols were selected as an indicator of spreading contamination based on their high water solubility and slow degradation (Broholm et al. 1999). Information from the analyses was applied for locating the source of the contamination and verification of the simulations.

Geophysical data
It was concluded that continuous and area covering data was essential to resolve the geology in detail for modelling of contaminant transport. Therefore, geophysical investigations were conducted on a 200 m × 330 m test site to construct a high-resolution geological model (Andersen et al. 2018). The extensive geophysical datasets detailed in Andersen et al. (2018) included DualEM-421, GPR, and DC/TDIP (Fig. 3). During the investigation, two ground-penetrating radar (GPR) profiles and ten direct current (DC) resistivity and time domain-induced polarisation (TDIP) profiles were measured. Moreover, a total of 5500 m of profiles was measured in two surveys with the DualEM-421.

Anthropogenic structures
The buried anthropogenic features expected to affect the flow of groundwater in the area include sewer systems, water pipes, cables, and houses with basements. Information on these were provided as georeferenced features by Horsens municipality. Cables are generally expected to be located no deeper than 0.8 m.b.t., while sewers and water pipes are placed from terrain to approximately 1.25 m.b.t. The depth of basements is assumed to be 1.8 m (Fig. 4).

Voxel model
The information from the geophysical mapping, boreholes, and location of anthropogenic features were imported to an interpretation software-GeoScene3D (www.I-GIS. dk)which conducted an integrated interpretation, as described in Andersen et al. (2018). As an alternative to the traditional surface modelling approach, the geological property variations and anthropogenic structures were represented in a regular 3D grid. Material parameters are assigned to each 3D cell (voxel) in the grid, allowing for heterogeneity, abrupt changes, and properties that vary within lithologies (Jessell 2001). In this study, the georeferenced information about anthropogenic structures has been included in the geological 3D voxel model. For the gasworks model, voxels can delineate the geological structures and characterise anthropogenic features within the uncertainty of the mapping methods. Structures smaller than the voxels can possibly exist, but they will be below the resolution of the applied mapping methods. In the following optimisation process, the voxel model was reviewed and adjusted by geologists to achieve the best match between the model presentation and geological interpretation.
The voxel model contains information about a 5-7-m layer below the crust. Because of topography and local information about greater depth, the model extends from 3.75 m above sea level to 11.25 m below sea level with voxels measuring 1 m × 1 m × 0.5 m (x, y, z). The topsoil fillings are dominated by sand and minor interactions of clay and silt ranging between 1 and 5 m in thickness. Under the fillings, a 1-3-m-thick layer of medium to coarse-grained sand is observed. Below the sand, a series of thin discontinuous 0.5-2-m-thick clay and silt deposits are observed. A 2-10-m-thick unit alternating from medium-grained sand to silt and clay layers (Fig. 4) lies under the clay and silt layers. Below this, boreholes indicate the presence of primarily moraine clay.

FEFLOW modelling
The high-resolution hydrogeological voxel model forms the basis for estimating the groundwater flow field and the contamination transport from the gasworks site to the harbour basin. A 3D steady-state FEFLOW (Finite Element Modeling of FLOW, mass and heat transport in porous and fractured media) model simulating confined groundwater flow and contaminant transport was utilised.

Model discretisation
The model domain covers 553,819.5 to 554,019.5 m Easting and 6,190,599.5 to 6,190,929.5 m Northing, as per European Terrestrial Reference System (ETRS) 89/Universal Transverse Mercator (UTM) zone 32 N in the horizontal plane. The elevation of the top and bottom of the model is 3.75 m above and 11.25 m below the mean sea level, respectively. The model is a cuboid spanning 200 m (UTM-X) × 330 m (UTM-Y) horizontally and 15 m in vertically. Voxel dimensions are uniformly set to 1 × 1 × 0.5 m (x, y, z). Thus, the model has a total of 1,980,000 elements (i.e. 66,000 elements in each layer) and 2,062,461 nodes. To simplify the calculation, part of the elements representing seawater, impermeable anthropogenic structures, and elements without geological information were set as computationally inactive. (1)

Model inputs
The initial and calibrated model parameters have been listed in Table 1.

Inflow on the top
Land use information was provided as georeferenced features by the Horsens municipality. Infiltration areas for numerical simulations were defined based on this information. The precipitation infiltration was the main inflow on top of the aquifer. Buildings and impermeable land cover such as roads were assigned zero infiltration. The recharge was computed by multiplication of the recharge coefficient with the annual precipitation data. The mean annual precipitation and soil infiltration coefficient were estimated to be 880-950 mm/year and 23.5%, respectively, by Severinsen et al. (1996) and Stisen et al. (2012). A uniform precipitation of 900 mm/year was applied for the study area which corresponds with an annual infiltration of 212 mm/year.

Boundary conditions
As shown in Fig. 5, on the interface between Wall 1 and seawater, a constant head boundary condition (BC) was set to 0 m over the sand unit. The northern boundary of the aquifers was assigned by the Kriging interpolation of observed groundwater levels for each aquifer (Bartier and Keller 1996;Yao et al. 2014). Due to the impenetrability of the structures located outside the southern study area, the southern boundary was set as a no-flow condition over the entire depth. The no-flow boundary was also defined along the eastern and western boundaries since the overall flow direction of the groundwater was north to south (Andersen et al. 2018). The bottom of the model was also considered as a no-flow boundary. At nodes where the model was not specifically defined, the default impermeable boundary condition was maintained.

Material properties
Voxel cells from the hydrogeological model were assigned with hydraulic conductivities on import to FEFLOW. The chosen hydraulic conductivity was based on lithological definitions, field observations, empirical values, and results of previous studies (Fredericia 1990;Henriksen et al. 2003;Heron et al. 1998;Sidle et al. 1998). The conductivity value of each unit was found to be spatially varying and anisotropic in the vertical direction, i.e. K xx = K yy ≠ K zz . The vertical hydraulic conductivity was assumed to be 10% of the horizontal conductivity in the x or y plane (Priyanka et al. 2018). The porosity was considered uniform and set to 0.3 for all units.

Flow model calibration
The flow model was calibrated by comparing observed and calculated hydraulic heads in the lower sand aquifer in boreholes using varying horizontal and vertical hydraulic conductivities. The accuracy of the calibrated models was estimated with the root-mean-square error (RMSE): Table 1 Input parameters for numerical simulation a The hydraulic head boundaries in the upper layers of sand were acquired from a regional model and interpolated by the linear method, while the boundaries in lower layers of sand were acquired from observation wells and interpolated by the Kriging method b Vertical hydraulic conductivity (K zz ) = 10%·Horizontal hydraulic conductivity (K xx , K yy ) where, H obs is the observed hydraulic head, H sim is the simulated hydraulic head, and n is the number of data points.

Contaminants transport
The initial mass concentration of the model was set at 0 to simulate the impact of the contaminants, which originated from the identified contamination source, on the whole model area. The mass boundary conditions on the northern nodes were assigned with a zero boundary since the groundwater flow direction was expected to be roughly southwards.
The Scheidegger-Bear model was used for the mechanical dispersion in which the transverse dispersion was assumed to be isotropic (Raats 1973). Both longitudinal and transverse dispersion coefficients depended on the spatial scale of the model domain and were uniformly set to 7.5 m and 0.75 m respectively (Gelhar et al. 1992). The contamination source corresponds to the hot spot area at the gasworks site, as shown by Andersen et al. (2018) (Fig. 4). In the vertical direction, boreholes indicated the presence of sandy sediments below the groundwater table. Moreover, resistivity and chargeability profiles from 10 DC/TDIP profile lines (Fig. 3) suggested that sandy sediments located in the saturated zone were heavily contaminated (Andersen et al. 2018). Building on these premises, the source of the contaminants was simulated with the elements from layer 5 to layer 7 that were predominantly sand (Fig. 4). Based on the groundwater analyses, a source area of phenols was outlined and assigned with a concentration of 100,000 µg/L. The geometry of the spreading of contaminants rather than the different compounds is of interest in this article; thus, results are shown for phenols only. The remaining components could follow a similar path, but with different properties that could be estimated by a chemist. The contaminant transport was assumed to be conservative.

Results and discussion
The hydrological model Figure 5 illustrates the boundary conditions in the model, the simulated hydraulic head distribution in the lower aquifer, and 52 observed hydraulic heads in the study area. The simulated contours show how the hydraulic head is highest near the northeast and gradually decreases from north to south, as it approaches 0 at Wall 1. Some discrepancy is observed between manual and simulated values as well as between closely lying manual measures. Observed values are generally lower than the ones simulated in the harbour area and higher than the simulated values at the north of the model. Manual measures have not been conducted simultaneously; some are years apart. It was concluded that measurements made at the same time appear to show similar levels, decreasing from north to south, towards the harbour front in accordance with the simulations. Boreholes with over five measures have values varying by more than a metre. With asynchronous measures, variations may originate from seasonal change, human activities, tides, and recent precipitation. Nevertheless, measuring boreholes with high-density distribution offers possibilities to acquire relatively authentic data representing the continuous groundwater flow field.  Figure 6 shows the modelled original distribution of the contaminant plume in the sandy aquifer. The values represent the maximum modelled concentration of contaminants in the aquifers presented below that point in the model. The original distribution of contaminants was simulated from known pollution source (line in red) and calibrated with the observed contour with the concentration of 10,000 µg/L (line in black). The dashed line in blue represents the Danish criteria of phenols in groundwater, thus approximately outlining the contamination plume. It can be noted that there are three green dots which represent no phenols detected are located in the edge of the plume This probably results from the conservatism of the calculation method. This figure represents the original scenario and will be maintained in the following figures for comparison to highlight the changes in the flow direction and extent of contaminants when the harbour Wall 1 is renovated.
The figure reveals a southward migration of contaminantsaway from the source. This is in accordance with the general north to south flow of groundwater. The plume widens and the concentrations decrease southwards. The areas located south of the source, including the gasworks area, streets, and the western part of the site 1, are the main contaminated areas while site 2 is at a lower risk of contamination. The permeable Wall 1, representing the shortest hydraulic distance to the sea, implies a discharge at the western part of the wooden harbour wall. Consequently, installing the outflow control facilities on the wall is recommended. This will drain the contaminated groundwater to prevent the expected spreading if the contaminated water cannot escape and is forced onto a longer path. Furthermore, it will have contaminants flowing out at a known position, ready for remediation work in the future.
The figure further compares the modelled plume and observations from groundwater samples (circles), soil samples (triangles), and pore air samples (crosses). The red symbols show locations where phenols were detected. Green symbols represent locations of analyses with no phenols detected. With phenol analyses carried out at different times, on different sample types and possibly using different methods, the figure shows whether phenols were detected rather than the detected phenol concentration. Phenols do not evaporate into the pore air easily. That might be the reason why some green crosses representing pore air samples were located where groundwater samples indicate presence of phenols. Furthermore, pore air analyses were collected between the depths of 0.3 and 0.8 m below the terrain, often in fillings or organic soils located above the aquifers from which water samples had been collected. Results from soil samples showed an unclear pattern. These samples were also collected at shallow levels, generally less than 1 to 2 m depth, typically in layers consisting of gyttja or fillings. The harbour area was a former industrial area with piles of various materials lying on the ground available to the industry needing it. Thus, it is expected that the findings in shallow soil samples originated from the local contamination related to former storage spaces.

Scenario and prediction
The outflow control facilities comprise a perforated drainage pipe supplied with two outlet pipes. The drainage pipe is 10 m long and installed parallel to the wall. The two 2-m-long outflow pipes connect the drainage pipe to the sea; see Fig. 7b-d.
Geophysical surveys and borehole information suggest the presence of two sand layers in the study area. An upper, 1-3-m-thick sand layer to the north thins out and becomes discontinuous to the south, pinching out near the Walls (Fig. 4). The lower sand layer gradually increases in thickness from north to south. Although the lower sand generally appears thinner than the upper sand layer, it seems to be continuous over the area and can provide ideal conditions for the groundwater and contaminants. The analysis of groundwater samples and electrical resistivity data revealed that both the sand layers were contaminated (Andersen et al. 2018). The upper sand layer failed to provide a continuous migration path for contaminants due to the pinch out before reaching the wall, where the facilities were supposed to be located. Hence, to avoid disturbing the original flow pattern, the outflow control facilities were installed in the lower sand layer.
Subsequently, three scenarios have been designed to determine the best location for the outflow control facilities (Fig. 7b-d). In the first two scenarios, the facilities have been placed on the western part of Wall 1 with the concentrated contaminant discharge into the sea, as shown in Fig. 7a. Conversely, the third scenario is placed on the Wall 2, to test whether installing the facilities on this wall can reduce the plume width. In these scenarios, the outflow control facilities are arranged at intervals of around 10 m. The outflow control facilities are implemented in the numerical model by assigning the nodes on the drainage pipe as a fixed hydraulic head (i.e. 0 m) since the drainage pipe is connected to the sea through the outlet pipes. Figure 7 (colour-filled contours) shows the distribution of the contaminants' plumes in the three scenarios with different locations of the outflow control facilities and compares them with the original case (grey area enclosed by black dashed lines). Compared to the original case, the covered areas of the contamination plumes in site 1 for all scenarios (Fig. 7b-d) were effectively reduced. Despite the largest reduced areas having been observed in scenario 3 (Fig. 7d), newly contaminated areas are predicted to derive from this case. This acts in contravention of the in situ remediation regulations. Thus, scenario 3 was discarded. In scenarios 1 and 2 (see Fig. 7b, c), the outflow control facilities were placed on the western part of Wall 1, where the concentrated contaminants originally discharged into the sea. Both scenarios appear to fulfil the requirements.
Scenario 1 may be the preferred solution because of the reduced contaminated area in site 1. Moreover, concentrations near the facilities in scenario 1 were higher than those in scenario 2. Therefore, more contaminants were discharged Fig. 7 Results of original situation and scenarios. (a) Original situation, Wall 1 was permeable and Wall 2 was impermeable. In (b) scenario 1 and (c) scenario 2, the outflow control facilities were both installed on Wall 1 but a different location. (d) Scenario 3, the outflow control facilities were installed on Wall 2. The colour-filled contour represents the contaminant plume distribution in different scenarios, while the binary filled contour showed in Fig. 7a and partly covered in Fig. 7b-d shows the original distribution of the contaminant plumes. Coloured arrows indicate the displacement of contours by the facilities installed in scenario 1 than those in scenario 2 during the same period. It can also be noted that in scenario 1, displacement of contaminated areas was mainly located in site 1. Although different scenarios were tested to minimise the impact of contaminants, it was inevitable that site 1 and the streets will still be located in the path of the plume.
The main purpose of installing the outflow control facilities is avoiding an additional spread of contaminants compared to the original situation while renovating Wall 1 and making it impermeable. Under this condition, the recommended solution would be installing facilities west of Wall 1, near Wall 2, as shown in Fig. 7b. In addition to fulfilling this requirement, the most important thought behind the planned solution is the installation of these facilities provides an accurate position for future collection and remediation work compared to the widespread discharge areas in the original case.
Phenols were chosen as the solute for transporting simulations due to its relatively higher mobility compared with other components among identified volatile aromatics. Moreover, the mass transport was conservatively set. There are several types of contaminants characterising different chemical reactions and physical processes in the soil and groundwater. However, acquiring the accurate composition and related parameters as well as systematic concentration testing work are absent. Therefore, based on the results of this study, subsequent research should not only focus on the collection and treatment of contaminated groundwater at outflow control facilities but also establish more accurate solute transport-reaction and degradation models, building on the relevant data and parameters.

Conclusions
This study suggests a workflow for planning the transformation of a polluted industrial area into a living space. The workflow includes data assessment and collection, and construction of a 3D hydrogeological model integrating anthropogenic features and information on infiltration. It also includes simulations and scenario testing designed to answer questions relevant for the decision about the transformation process. The main point of the workflow is that detailed information is needed for reliable results. Thus, geophysical mapping of the subsurface is recommended to tie together information between boreholes. In urban areas, abrupt changes from construction works and anthropogenic features call for voxel models rather than the traditional layer models in order to better represent these features. Finally, scenario calculations allow for testing of the influence of different construction types on the contaminant transport. The workflow has been successfully applied at Horsens harbour, Denmark. Simulations have been verified by comparing the models and the available data from the site. Scenario testing designed for decision making was undertaken. The results of numerical simulations suggest the optimal location for installing outflow control facilities guiding a flow of contaminants. By knowing the location of the outflow of the contaminants, they are ready for treatment and remediation work. The installation scheme effectively mitigates the impact of contaminants on two planned construction sites. This study demonstrates the way modelling of contaminant transport based on detailed hydrogeological information for different scenarios can assist decision making about complex urban contaminated sites.