Numerical analysis of the transport of brine in the Odra River downstream of a mine's discharge

Underground mining and copper processing is one of the industries that use significant amounts of water [1]. A characteristic feature of most mines is the inflow of water from the surrounding rock mass, or from aquifers located above, from which the water flows through cracks that occur due to mining operations. To reduce the environmental impact, mine waters are currently used in processes related to the extraction and processing of solid minerals and are then stored together with mining wastes. In water-rich locations, the use of mine water may not balance the inflows, and any excess must be disposed in another way. The most commonly used method is the hydrotechnical method, which consists of discharging mine water into rivers [2, 3]. The design of discharge devices depends on the amount of water being discharged, the hydrological characteristics of the river, the geometrical parameters of the riverbed, as well as the variability of flows and the manner of their use, including the purposes of navigation and recreation. According to [4], the most common method of discharge construction involves the location of an outlet on one of the riverbanks. The shape and length of the contamination zone (until it is completely mixed) depend on the outlet's construction (Fig. 1a). Another way is to lay a steel or plastic pipeline on the river's bottom, perpendicular to the axis of the riverbed, with holes arranged, so that the outflow is as even as possible. However, such a solution is more prone to mechanical damage, including that resulting from the influence of sediment dragged along the river's bottom (Fig. 1b).

(a) Examples of mixing zone for river discharges: (A) a traditional shoreline discharge, (B) a submerged offshore single-port discharge, (C) a regular perpendicular diffuser design and (D) a diffuser design for navigating rivers or rivers with other cross-sectional limitations. B_MZ – the width of mixing zone, L_MZ – the length of the mixing zone downstream. B_MZ at distance L_MZ must be greater than B_MZ specified in the water legislation [4]. (b) Photo of a fragment of the discharge installation damaged by sediment (after disassembly).

Reducing the mixing zones of mine water in a river is essential for protecting the natural environment. The composition of mine water differs from that of surface waters. In the literature [5, 6], mine water is categorised into acid mine drainage (AMD), with a pH below 6, neutral mine drainage (pH above 6) and saline drainage (SD) with a pH above 6, and total dissolved solids (TDS) above 1 kg m⁻³. AMD is extremely hazardous to the environment due to its heavy metal content and the need for it to be treated to an environmentally safe level. The discharge of sewage, regardless of its type, to natural watercourses in the European Union is regulated in accordance with the Water Directive and must meet the requirements of emission limit values and environmental quality standards (EQS). The introduction of EQS requires the discharge operator to adjust the pollutant plume, so as to not exceed the values that are considered safe. In the zone where pollutants mix with surface water, EQS are exceeded, which in turn has a negative impact on the environment.

SD is characterised by a high content of TDS. Due to the geological structure of polish deposits, a large part of them, for example, hard coal [7] and copper ore, is susceptible to an inflow of salt mine water. The highest load of pollutants is related to the high content of chlorides and, to a lesser extent, sulphates. The concentration of chloride in rivers negatively affects their fauna and flora. The strength of this interaction depends on the species, temperature and duration of the brine exposure and its concentration. The most sensitive species of aquatic plants negatively react to the concentration of 0.1 kg m⁻³; however, fish species show resistance to salinity of several kg m⁻³ [8]. To limit the negative impact of increased salinity in the mixing zone, it is advisable to define the maximum range of the mixing zone [9, 10]. Discharge installations should be constructed in such a way as to ensure the fastest possible mixing [4].

The use of standard river quality assessment methods (along the river current) in the case of mine water discharges may lead to misinterpretation. High concentration zones may bypass the sampling point or vice versa – they may suggest much more contamination than there actually is. The correct assessment of pollutant concentrations downstream the discharge requires the performance of concentration tests in research measuring points located across the river. Additionally, the measurements must be related to the value of the water flow velocity [11].

The article describes the installation of KGHM Polska Miedź S.A. on the Odra River in Głogów, which consists of a perforated DN1420 pipeline laid on the bottom of the river and auxiliary equipment. The form of the discharge installation results from the fact that the Odra River is used as a waterway, and therefore, there is a need to ensure an appropriate depth. The aim of the article is to present 3D numerical modelling results of brine mixing process in the river downstream mine water discharge and their comparison to the measurement results. The modelling results and measurements will also determine whether the applied environmental impact monitoring system is designated correctly. This information may help relevant services to make decisions regarding the management of saline waters and the reduction of salinity in the Żelazny Most tailings storage facility (TSF). Thus, the undertaken decisions will be able to protect the aquatic environment against the negative effects of the introduced pollution more effectively. It is a comprehensive approach to determining the quality of mixing of the mine water stream in the Odra River.

Study site

2.1

Brine characteristics

KGHM exploits a single copper seam of sedimentary origin in the south-western part of Poland. The deposit was made available from its shallowest side in the region of Lubin and Polkowice, and further exploitation of its deeper areas is mainly carried out towards the north. Water inflow is from the surrounding rock mass to the workings. In the shallowest regions of the deposit, the water inflow is characterised by the highest flow rate and a relative low mineralisation. In the areas located further to the north, the outflow of water from the rock mass is much lower, but this water contains high concentrations of salt. It comes from a rock salt seam located on average 80 m [12] above the pits used for copper extraction. Among the mine waters that are pumped out from the three mining plants of KGHM, the waters from Rudna mine have the highest salinity of approx. 150 kg m⁻³. The Lubin and Polkowice–Sieroszowice mines pumping stations, which are characterised by lower mineralisation, reduce the average salt concentration in the process water to the level of 38 kg m⁻³. The salinity of the water increases its density from 2% to 10%, and there is a need to take this parameter into account when selecting pumps and fittings, as well as during numerical analyses of the discharge installation. It is the only factor that can have a negative impact on the environment, and it prevents the direct unlimited discharge of mine water into water courses.

2.2

Mine water management in KGHM

Mine water that is pumped to the surface goes to equalising tanks at Ore Enrichment Plants (OEP), where it is used in enrichment processes. As a result of these processes, about 29 million Mg rock is produced, which is devoid of copper compounds and other useful metals, commonly known as mining waste. It is pumped in the form of a hydro-mixture to the Żelazny Most TSF, which is the only place where mining waste is stored for KGHM. On beaches and in the TSF reservoir, the solid parts sediment, which allows for the storage of clean water by pumping it back to the OEP for reuse (Fig. 2). Despite the repeated use of mine water in processing, its constant inflow makes it necessary to manage its excess. Due to the climatic conditions and the volume of the stream, the only technically justified method of naturalising the mine water is its discharge into natural watercourses.

A diagram of the flow of the KGHM mine water [13].

2.3

Brine discharge installation

The nearest river with a sufficiently large flow that can absorb excess water from mine deodorisation, without damaging the environment, is the Odra River. Pipelines that had been originally used to supply water (for production) from the Odra River were used to drain water from TSF. Along with the increase in the range of mining excavations, the mine water fully covered the demand for the technological water of KGHM. Currently, KGHM records approx. 25 million m³ p.a. of excess mine water. The previously existing pipelines were connected with a bottom discharge located approximately 420 m upstream the bridge along the route of national road No. 12 in Głogów. Until 2021, an installation in the form of a perforated pipe in the bottom of the river across its current was in operation; but due to the too shallow foundation on the river bottom, which caused its extensive mechanical damage, a decision was made to build a new installation several dozen metres upstream. The currently operated discharge consists of two reinforced concrete chambers on both banks of the river and a pipeline connecting them with a diameter of 1420 mm, which is laid at a depth of about 3 m beneath the bottom of the Odra River (Figs 3 and 4). The brine is delivered to the chamber on the left bank of the river. Such a deep location of pipe foundation is caused by the necessity to increase the protection against the possibility of the installation being exposed to river traffic, as the authorities are planning to increase the navigability of the Odra River. In the 50-m-long part, the pipeline is equipped with multi-nozzle ejectors to ensure an even outflow of technological water. With the design capacity of the entire system being 2 m³ s⁻¹, one dispersing nozzle works with a flow rate of 0.0625 m³ s⁻¹.

The discharge system during construction.

Cross section of the discharge installation, dimensions in mm [14].

2.4

Characteristics of the Odra River and the principles of brine discharge

The Odra River is the second largest river in Poland, the source of which is in the Czech Republic. It flows into the Baltic Sea. In Głogów, about 2.75 km downstream the discharge installation, there is a permanent water gauging station belonging to the Institute of Meteorology and Water Management. The data obtained from the gauging station shows that the lowest recorded flow rate is 41 m³ s⁻¹ and the average annual flow is 198 m³ s⁻¹, which is 100 times more than the maximum capacity of the discharge installation.

The quality tests of Odra water are periodically carried out upstream the discharge to determine the volume of brine that can be introduced into the Odra River through the discharge installation because the maximum concentration of the sum of chlorides and sulphates downstream the discharge cannot exceed 1 kg m⁻³. The salinity of the mine water is tested daily in a laboratory, while the natural salinity of the Odra River is tested three times a month upstream the discharge. Studies conducted over many years indicate the dependence between natural salinity and the flow rate of the Odra River (Fig. 5). A similar relationship was also found in the studies of other authors [15,16,17]. This correlation is used to determine the flow rate of brine that can be introduced into the river. Every day, KGHM specialists verify the maximum discharge flow rate on a given day on the basis of the table of permissible discharge flow rates, which depends on the water level in the Odra River and the salt concentration in the discharged waters. In order to control the effectiveness of the discharge installation, additional water quality tests are carried out 422 m downstream the discharge, where three water samples are taken from the left and right riverbanks as well as from the middle of the Odra River's current. Samples are collected from just below the water surface and then tested by an independent company.

The relationship between the flow of the Odra River and its natural salinity.

Problem statement and formulation

Once brine is discharged into a river, it is transported and transformed by physical processes. During the brine transport, its concentration decreases due to the phenomenon of turbulent mixing and dissolution [18]. Therefore, the concentration of brine along the length of the stream will be locally variable and dependent on the geometry of the riverbed and discharge, the brine outflow velocity and the field of turbulent flow velocity in the riverbed. Physically, the initial concentration as well as the processes of advection and diffusion are responsible for the phenomenon of changes in the concentration over time and space. Their mathematical description is the equation of advection and diffusion of a scalar quantity, that is, the concentration or density of the pollutant [19]. To determine the components of the velocity field of the turbulent flow, which carries the pollutant particles, the methods of solving the system of Navier–Stokes equations are used. The description of the assumptions made and the solution of the problem of transporting pollutants are presented in the following sections of the article.

3.1

Mathematical formulation

3.1.1

Advection-diffusion of pollutants

The general form of the diffusion–advection equation of the pollutants in the area Ω bounded by the bank ∂Ω = S_B ∪ S_F, which consists of the surface of the bottom of the river and riverbanks S_B and the free water surface S_F, can be represented as follows: (1) $\frac{\partial}{\partial t} C (t, x) + (u (x) \nabla) C (t, x) - (\nabla \cdot D (x) \cdot \nabla) C (t, x) + f_{c} (x) C (t, x) = q (t, x)$ {\partial \over {\partial t}}C\left( {t,{\bf{x}}} \right) + \left( {{\bf{u}}\left( {\bf{x}} \right)\nabla } \right)C\left( {t,{\bf{x}}} \right) - \left( {\nabla \cdot {\bf{D}}\left( {\bf{x}} \right) \cdot \nabla } \right)C\left( {t,{\bf{x}}} \right) + {f_c}\left( {\bf{x}} \right)C\left( {t,{\bf{x}}} \right) = q\left( {t,{\bf{x}}} \right) where C(t,x) is the concentration of pollutants determined for t > 0 and x ∈ Ω, x = [x, y, z] is the position vector with Cartesian coordinates (the x axis is horizontal and consistent with the river flow direction and the z axis is vertical), D(x) = [D_x, D_y, D_z] is the turbulent diffusion tensor, u(x) is the vector of the turbulent velocity of the water flow, f_c (x) is a function that describes the velocity of the chemical process that changes the concentration of the pollutants and q(t,x) is the intensity of the source of the pollutants. The solution of the advection-diffusion equation must be combined with the solution of the problem of momentum conservation for the turbulent flow, which is described by the Reynolds equation, the mass conservation equation and the corresponding initial and boundary conditions.

3.1.2

The Reynolds-averaged Navier-Stokes–Renormalization Group method

The 3D modelling of the flow in an open channel consists of determining the distribution of velocity and pressure, as well as the position of the free surface. This movement is described by a system of equations, which includes the mass conservation equation and the momentum conservation equation that are supplemented with appropriate initial and boundary conditions. The mathematical description of turbulent flow usually comes down to averaging the system of equations with respect to time – the so-called Reynolds-averaged Navier–Stokes (RANS) method – or to appropriate filtering of the components of velocity fluctuations – the so-called large eddy simulation (LES) method. In this paper, the RANS method was chosen due to its lower computational requirements.

Turbulence is characterised by fluctuations in velocity and in other physical quantities that occur in the flow. Reynolds formulation involves the decomposition of any physical quantity ϕ into an average value and its fluctuation [20]: (2) $ϕ = \bar{ϕ} + ϕ^{'}$ \phi = \bar \phi + \phi ' where $\bar{ϕ}$ {\bar \phi } is the time-averaged value of a scalar quantity and ϕ′ is the fluctuation of that quantity.

After applying the averaging procedure to all quantities that appear in the Navier–Stokes equations, mass conservation equation and the boundary and initial conditions, a system of equations known as the Reynolds equations is obtained: (3) $ρ \frac{\partial \bar{u}}{\partial t} + ρ (\bar{u} \nabla) \bar{u} = - \nabla \bar{p} + μ \nabla^{2} \bar{u} - ρ \nabla \cdot (\bar{u^{'} u^{'}}) + f$ \rho {{\partial {\bf{\bar u}}} \over {\partial t}} + \rho \left( {{\bf{\bar u}}\nabla } \right){\bf{\bar u}} = - \nabla \bar p + \mu {\nabla ^2}{\bf{\bar u}} - \rho \nabla \cdot \left( {\overline {{\bf{u}}'{\bf{u}}'} } \right) + {\bf{f}} (4) $\nabla \cdot \bar{u} = 0$ \nabla \cdot {\bf{\bar u}} = 0 where $\bar{u}, \bar{p}$ {\bf{\bar u}},\bar p are the averaged values of velocity and pressure, μ is the dynamic viscosity coefficient, f is the vector of mass forces (in this paper, only the force of gravity is taken into account: f = ρg; g is the acceleration due to gravity) and $ρ \bar{u^{'} u^{'}}$ \rho \overline {{\bf{u}}'{\bf{u}}'} is the Reynolds stress tensor. This is a symmetric tensor and its values are the components of turbulent stress. There are 10 variables in equation (3) – three components of the average velocity and pressure and six terms of the Reynolds stress tensor. Therefore, apart from the equations of conservation of momentum (equation 3) and mass (equation 4), additional equations are needed to solve equation (3) unequivocally. In this study, a turbulence model was used, which involves the use of two additional equations, the conservation of kinetic energy of turbulence k and the dissipation rate of kinetic energy ɛ. The RNG (Renormalization Group) model was derived from the classical model k − ɛ that describes the flow of an incompressible fluid with isotropic turbulence and a relatively large Reynolds number, but with a modified coefficient C_ɛ2. This model is given by the following system of equations: (5) $\frac{\partial k}{\partial t} + \bar{u} \cdot \nabla k = ℘ - ε + \frac{\partial}{\partial x_{i}} (\frac{ν_{T}}{σ_{k}} \frac{\partial k}{\partial x_{i}})$ {{\partial k} \over {\partial t}} + {\bf{\bar u}} \cdot \nabla k = \wp - \varepsilon + {\partial \over {\partial {x_i}}}\left( {{{{\nu _T}} \over {{\sigma _k}}}{{\partial k} \over {\partial {x_i}}}} \right) (6) $\frac{\partial ε}{\partial t} + \bar{u} \cdot \nabla ε = \frac{\partial}{\partial x_{i}} (\frac{ν_{T}}{σ_{ε}} \frac{\partial ε}{\partial x_{i}}) + C_{ε 1} \frac{ε}{k} ℘ - C_{ε 2} \frac{ε^{2}}{k}$ {{\partial \varepsilon } \over {\partial t}} + {\bf{\bar u}} \cdot \nabla \varepsilon = {\partial \over {\partial {x_i}}}\left( {{{{\nu _T}} \over {{\sigma _\varepsilon }}}{{\partial \varepsilon } \over {\partial {x_i}}}} \right) + {C_{\varepsilon 1}}{\varepsilon \over k}\wp - {C_{\varepsilon 2}}{{{\varepsilon ^2}} \over k} where $℘ = ν_{t} (\frac{\partial {\bar{u}}_{i}}{\partial x_{j}} + \frac{\partial {\bar{u}}_{j}}{\partial x_{i}}) \frac{\partial {\bar{u}}_{i}}{\partial x_{j}}$ \wp = {\nu _t}\left( {{{\partial {{\bar u}_i}} \over {\partial {x_j}}} + {{\partial {{\bar u}_j}} \over {\partial {x_i}}}} \right){{\partial {{\bar u}_i}} \over {\partial {x_j}}} is the tensor describing the production of turbulence, (7) $ν_{T} = C_{μ} \frac{k^{2}}{ε}$ {\nu _T} = {C_\mu }{{{k^2}} \over \varepsilon } (8) $C_{ε 2} \equiv C_{ε 2}^{*} + \frac{C_{μ} η^{3} (1 - η / η_{0})}{1 + β_{1} η^{3}}$ {C_{\varepsilon 2}} \equiv C_{\varepsilon 2}^* + {{{C_\mu }{\eta ^3}\left( {1 - \eta /{\eta _0}} \right)} \over {1 + {\beta _1}{\eta ^3}}} The RNG model contains an additional term in C_ɛ2 that is small in weakly strained turbulence and large in rapidly distorted flows [21]. (9) $η = \frac{Sk}{ε}, S = \sqrt{2 S_{ij} S_{ij}}, S_{ij} = \frac{1}{2} (\frac{\partial {\bar{u}}_{i}}{\partial x_{j}} + \frac{\partial {\bar{u}}_{j}}{\partial x_{i}})$ \eta = {{Sk} \over \varepsilon },\;\;\;S = \sqrt {2{S_{ij}}{S_{ij}}} ,\;\;\;{S_{ij}} = {1 \over 2}\left( {{{\partial {{\bar u}_i}} \over {\partial {x_j}}} + {{\partial {{\bar u}_j}} \over {\partial {x_i}}}} \right) where ν_T is the eddy viscosity coefficient, S_ij is the rate of the deformation tensor and S is the average strain velocity in the flow. The non-dimensional parameter η is the ratio of the turbulence time scale κ/ɛ to the time scale of the mean S and has a large value in regions of rapid distortion. The remaining dimensionless coefficients were derived analytically and experimentally verified for various cases of turbulent flows [21]. Their values are as follows: (10) $\begin{matrix} C_{μ} = 0.085, C_{ε 1} = 1.42, C_{ε 2}^{*} = 1.68, σ_{k} = \\ = σ_{ε} = 0.72, β_{1} = 0.012, η_{0} = 4.38 \end{matrix}$ \matrix{ {{C_\mu } = 0.085,{C_{\varepsilon 1}} = 1.42,C_{\varepsilon 2}^* = 1.68,{\sigma _k} = } \cr { = {\sigma _\varepsilon } = 0.72,{\beta _1} = 0.012,{\eta _0} = 4.38} \cr } The RNG method enables the solution of the Navier–Stokes equations on relatively thicker grids for the flow with a high Reynolds number value. Its advantage is better mapping of the length of the recirculation zones in flows with separation from obstacles, and thus better mapping of the mixing of pollutants in the turbulent flow stream.

3.1.3

Boundary conditions and the numerical program

The system of equations (1)–(3) must meet the following conditions:

initial: –

C (0, x) 0= The concentration of dissolved salt in the river is negligibly low

–

u (0, x) = u_xi The uniform velocity distribution at the inlet

boundary: –

(n·∇) C (t, x), = 0 on ∂Ω The normal component of the concentration gradient vanishes on the boundary

–

$\bar{u} = 0$ {\bf{\bar u}} = 0 on S_B ‘No-slip’ condition

–

q (t > 0, x), = q_d on A_d The discharge of the brine installation is constant

–

C (t > 0, x) = C_d on A_d The concentration of the brine discharge is constant

–

$u \cdot \hat{n} = u_{d} = q_{d} / A_{d}$ {\bf{u}} \cdot {\bf{\hat n}} = {{\bf{u}}_d} = {q_d}/{A_d} The velocity vector of the brine seeping out of the river bottom is constant and normal

where

\hat{n}

{{\bf{\hat n}}}

is the normal vector to the boundary, u_d is the velocity of the brine's outflow through the riverbed at the discharge site, A_d ⊂ ∂Ω is the part of the bottom surface area through which the brine flows out and q_d is the intensity of the source of the brine.

The existence of a free surface implies that the position of this boundary part ∂Ω is not generally known for t > 0 and must be defined as part of the solution. Therefore, additional equations should be formulated for this surface. The following dynamic (equation 11) and kinematic (equation 12) equations can be used here. (11) $T \cdot \hat{n} = - p_{0} \hat{n} on S_{F}$ {\bf{T}} \cdot {\bf{\hat n}} = - {p_0}{\bf{\hat n}}\;\;\;{\rm{on}}\;{S_F} (12) $\frac{\partial F}{\partial t} + \nabla \cdot (u F) = 0$ {{\partial F} \over {\partial t}} + \nabla \cdot \left( {{\bf{u}}F} \right) = 0 The dynamic equation expresses the equilibrium of normal forces on the free surface, where $T = ν_{T} (\nabla \bar{u} + \nabla {\bar{u}}^{T}) - p I$ {\bf{T}} = {\nu _T}\left( {\nabla {\bf{\bar u}} + \nabla {{{\bf{\bar u}}}^T}} \right) - p{\bf{I}} is the stress tensor, p₀ is the atmospheric pressure and $\hat{n}$ {{\bf{\hat n}}} is the normal vector to the free surface. Surface tension has been ignored in this equation. The kinematic equation is expressed in the volume-of-fluid (VOF), where the position of the fluid is expressed as a function of the volume of the fluid F(x, t)

This equation is associated with the method of numerical solution of the flow problem, which is the finite-volume method. The function takes the value of 1 in the cells occupied by the fluid (full cell) and the value of 0 everywhere outside the fluid region (empty cell). A cell, which is contained in the stationary edge, is blocked for the calculations. If the function in a cell has a value in the interval (0,1) and if this cell is adjacent to at least one cell with a value of 0, it means that a free boundary runs through this cell.

In the proposed numerical method of analysing the propagation of pollutants in the turbulent flow stream, the diffusion of the pollutants’ concentration is proportional to the momentum diffusion by means of the coefficient α ranging from 0 to 1.

The FLOW-3D program was used to solve the problem of brine transport in the riverbed, which is considered as the best tool for simulating free surface flows. The program uses the numerical method of finite-volume VOF, which was transformed into the current best algorithm for tracing free surface flow – TruVOF (used in the FLOW-3D program). A characteristic feature of this software is also the free definition of numerical grids – the program modifies their edges to accurately map the geometry of the area. Its other advantage is the ability to use several grids (including nested grids) with different grid sizes to increase the accuracy of calculations in areas with larger gradients of velocity and pressure fields. The calculations were made using the Reynolds turbulent flow equation with the RNG model and the transport equation of a scalar quantity, that is, brine concentration.

3.2

The computational grid and data

The numerical model covers a river section consisting of an upstream part 78 m long and a downstream part 422 m long from the discharge. The discharge area measures 50 m × 8.5 m = 425 m² (Fig. 6).

Three-dimensional model of the Odra riverbed, with the brine outflow zone marked.

The calculations were made with the assumption of keeping the maximum discharge expenditure that is permitted by law – in order to use the entire available absorption capacity of the river (determination of the absorption capacity is described in Section 2.4). Depending on the water level, the absorptivity of the river increases and decreases, and therefore, calculations were made for various flows of the Odra (41, 51, 77.2, 113, 118, 169 and 318 m³ s⁻¹). The lowest flow is the same as the historical lowest level of the Odra River. Other flow rates were determined by increasing the river depth by steps of 0.5 m. The absorptive capacity of the river for the two largest flows (169 and 318 m³ s⁻¹) is so large that the KGHM can discharge a maximum brine flow rate of 2 m³ s⁻¹.

The numerical grid of the model was created by dividing the calculation area into 750 elements in the x-direction, 150 elements in the y-direction and 14–22 elements in the z-direction, depending on the water depth in the riverbed. The results of the brine concentration calculations are presented in the form of figures of longitudinal profiles and cross-sectional profiles of the riverbed regarding the flow rate.

The momentum diffusion multiplier used in the numerical calculations was determined using the factor α = 1 (see Section 3.1.3). This is due to the fact that the viscosity of brine is of the same order as that of water [22]. The value of the function. f_c = 0 was assumed because the brine does not react chemically with the water and sediment in the river.

Results and discussion

4.1

Calibration and validation of the model

To verify the convergence of the solution, a change in grid density due to the grid size was applied. The following values of the dimensions of the volumetric element were assumed: (i) Δx = Δy = 1.0 m and Δz = 0.5 m, (ii) Δx = Δy = 0.667 m and Δz = 0.3 m and (iii) Δx = Δy = 0.5 m and Δz = 0.25 m. The calculations were performed for Q = 169 m³ s⁻¹, q_d = 2.0 m³ s⁻¹ and C_d = 150 kg m⁻³. The results of the maximum concentration value and the component values of the velocity vector u_x in the cross section – in the middle and at the end of the calculation section of the river – are presented in Table 1. A comparison of these values shows that the difference between the maximum concentration and the velocity for grids (ii) and (iii) is 1.5% and 0.7%, respectively. It was, therefore, decided to carry out further analyses for the computational grid (ii), and not for case (iii), as the computation time for case (iii) is very long. A comparison of the values for cases (i), (ii) and (iii) is shown in Fig. 7.

Table 1

The results of the maximum concentration value and the component values of velocity vector u_x.

Distance	Case	(i)	(ii)	(iii)	Error (ii) - (i), %	Error (iii) - (ii), %
X = 250 m	C_s (kg m⁻³)	0.991	0.999	0.997	0.80	−0.14
X = 250 m	u_x (m s⁻¹)	0.999	1.000	1.000	0.12	−0.13
X = 500 m	C_s (kg m⁻³)	2.022	2.140	2.173	5.52	1.52
X = 500 m	u_x (m s⁻¹)	1.442	1.448	1.437	0.43	−0.71

Image of the brine concentration in the cross section at the end of the model: (a) grid 1 × 1 × 0.5 – 500,000 elements, (b) grid 0.67 × 0.67 × 0.31 – 1,912,500 elements and (c) grid 0.5 × 0.5 × 0.25 – 4,400,000 elements.

The validation of the model was carried out on the basis of salinity measurements which were conducted in the cross section located approximately at the inlet to the calculation area, that is, about 78 m upstream the discharge, and in the cross section approximately at the end of the model area, that is, approximately 422 m downstream the discharge. Measurements were taken from a boat using a portable multiparameter meter Hanna HI9829. To verify the performed salinity measurements, the collected samples were tested in an accredited laboratory. The measured salt concentration values are presented in Table 2.

Table 2

Measured values of the brine concentration at selected depths in the Odra River upstream and downstream the discharge.

Location		Salinity of the Odra River at selected depths

		kg m⁻³	kg m⁻³	kg m⁻³	kg m⁻³

		m	m	m	m
80 m upstream the discharge		0.188	0.188	0.181	0.185
80 m upstream the discharge		2.2	1.5	0.6	0.0
420 m downstream the discharge	Left riverbank	0.614	0.594	0.594	0.565
	Left riverbank	2.5	2.0	0.6	0.0
	Middle of the river's current	0.549	0.566	0.565	0.565
	Middle of the river's current	2.5	1.5	0.5	0.0
	Right riverbank	0.266	0.246	0.263	0.257
	Right riverbank	1.5	1.0	0.5	0.0

In addition, the exact ordinates of the water level at the ends of the model section were measured from the riverbank using a GPS geodetic device. The water level elevation at the inlet was 70.80 m asl and at the outlet was 70.70 m asl, with the stage difference being 0.10 m. This value was necessary to create a numerical model in FLOW-3D. The flow rate of 118 m³ s⁻¹ was determined from the water gauge data. The brine discharge was 1.0 m³ s⁻¹ and its chloride concentration was 40 kg m⁻³.

Brine propagation simulation was performed and the results are presented in Table 3. The salinity measurements were reduced by the natural salinity of the Odra River upstream the discharge and these results are also presented in Table 3.

Table 3

Brine propagation simulation results (C_s) compared with the measurements of chloride content (C_m) at the end of the model section.

Left riverbank				Middle of the river's current				Right riverbank

h	C_s	C_m	Δ	H	C_s	C_m	Δ	h	C_s	C_m	Δ

m	kg m⁻³	kg m⁻³	kg m⁻³	m	kg m⁻³	kg m⁻³	kg m⁻³	m	kg m⁻³	kg m⁻³	kg m⁻³
0.0	0.400	0.380	0.020	0.0	0.287	0.380	−0.093	0.0	0.218	0.072	0.146
0.5	0.401	0.413	−0.012	0.5	0.287	0.384	−0.097	0.5	0.219	0.082	0.137
1.5	0.407	0.406	0.001	1.5	0.290	0.379	−0.089	1.0	0.220	0.059	0.161
2.5	0.410	0.426	−0.016	2.5	0.292	0.361	−0.069	1.5	0.224	0.078	0.146

The error observed on the right bank may be related to the failure to meet the assumptions of uniform distribution of water velocity in the initial cross section in the simulation model. Another factor influencing the model error may be the brine discharge value q_d = 1.0 m³ s⁻¹ on the measurement day due to the temporary condition of KGHM's water supply network. The installation was designed for an even outflow of brine with a capacity of 2.0 m³ s⁻¹ and for smaller flows, the assumed uniformity may not be met and the outflow velocity from nozzles located close to the right bank of the river may be lower.

As can be seen from the table, the differences are not significant and justify the conclusion that the constructed numerical model reflects the phenomenon of transport and mixing of brine in the flow stream in the Odra River with a satisfactory accuracy. The difference in the form of percentage error for the entire cross section showed a value of 6.89%.

4.2

Results of calculations

Calculations were performed for the following condition: six flow rates of the Odra River ranging from 41 to 318 m³ s⁻¹, in which the lowest value corresponded to the low flow that occurred on the Odra River in 2015 and the highest value corresponded to a flow rate 50% higher than the average annual flow rate. Higher flows were not analysed due to the negligible influence of brine on water quality. The boundary conditions for the modelling were: the water ordinate at the inlet H_i and the average velocity u_xi determined on the basis of the Odra flow rate Q_i and the cross-sectional area A_i at the inlet, the ordinate of the water level at the outlet H_o and the brine outflow velocity u_zd through the bottom in the area of discharge. It was assumed that the velocity vectors are directed towards the z-axis and the value of the mean flow rate of the brine was determined from the quotient of the intensity of the discharged brine q_d and the area of its outflow through the bottom with the value of 425 m². The value of q_d is adopted on the basis of the table of admissible discharges used by KGHM, which depends on the water level of the Odra River and the salt concentration in the discharged water. The last data that was taken into consideration was the concentration of the discharged brine C_d, which was the same for all analysed brine cases and was equal to 40 kg m⁻³ – average salt concentration in water discharged at the turn of 2020 and 2021. The values of the input data and the results of the mean and maximum concentrations in the outlet cross section are presented in Table 4.

Table 4

Input data and results.

No.	Q_i	A_i	u_xi	q_d	u_zd	C_d	C_avg	C_max

	m³ s⁻¹	m²	m s⁻¹	m³ s⁻¹	m s⁻¹	kg m⁻³	kg m⁻³	kg m⁻³
1	41	63.58	0.645	0.282	0.0007	40	0.2729	0.341
2	51	75.25	0.678	0.425	0.0010	40	0.3306	0.341
3	77.2	106.11	0.728	0.767	0.0018	40	0.3933	0.497
4	113	139.07	0.813	1.467	0.0035	40	0.5125	0.632
5	169	173.75	0.973	2.000	0.0047	40	0.4678	0.560
6	318	246.77	1.289	2.000	0.0047	40	0.2500	0.340

Figure 8 shows the spatial distribution of brine for the two extreme flow rates. Figure 9 shows a comparison of the brine distribution in the cross section at the end of the model for all flows. This comparison shows that the brine discharge strategy used by KGHM's specialists is appropriate because the obtained average concentration of chlorides did not exceed the permissible value equal to 1 kg m⁻³.

Spatial distribution of chloride for low flow 41 m³ s⁻¹ at elevation (a) 68.45, (c) 69.08, (e) 69.38 m asl and for high flow 318 m³ s⁻¹ at elevation (b) 68.45, (d) 69.08, (f) 72.13 m asl.

Distribution of chloride in the cross section at the end of the model for flow values: (a) Q_i = 41 m³ s⁻¹, (b) Q_i = 51 m³ s⁻¹, (c) Q_i = 77.2 m³ s⁻¹, (d) Q_i = 113 m³ s⁻¹, (e) Q_i = 169 m³ s⁻¹ and (f) Q_i = 318 m³ s⁻¹.

The efficiency of the mixing of brine in the turbulent flow depends on the river flow rate and the geometry of the bed. The studied section of the river is slightly curved in plan, and the main flow of the stream is along the left concave bank. Therefore, the concentration in this area is greater than at the opposite bank. The flow rate also affects the spatial density distribution of the brine. The higher the flow, the more complete the mixing of brine and water. This was presented graphically for the outlet cross section in the form of diagrams of the product of concentration and the fraction share of the area occupied by this concentration (Fig. 10).

Distribution of the concentration of brine in the outlet cross section for (a) Q_i = 41 m³ s⁻¹ and (b) Q_i=318 m³ s⁻¹.

Conclusions

The article discusses the results of calculations of the spatial distribution of pollutants in the river downstream of the discharge from the copper mine. In the case of discharges of salt mine waters into rivers, an important element of environmental protection is the determination of the distance at which the pollutants mix with river waters to a level that is not harmful to fish fauna and benthos. The study of pollutant propagation in the river was performed using numerical modelling of the advection–diffusion equation, whereas the turbulent flow velocity field in the river was modelled using Reynolds equations and the RNG turbulence model.

The greatest calculation error of the brine concentration in the outlet cross section was at the right bank of the river, which may be due to a different velocity distribution in the inlet section of the model (uniform distribution) than it was in the river during the measurements. Another reason for this discrepancy may be the uneven flow velocity from the nozzles along the length of the discharge line during measurements. The brine moves with the main stream, and if it is concentrated in only one area of the riverbed, the brine flow will also be concentrated within this stream and the mixing distance will be significantly longer.

The flow velocity in the river is also important for the efficiency of mixing, not only because of the proportionally higher absorption capacity for higher flow rates, but also because of the greater level of turbulence in the current.

Field tests as well as calculation results indicate a discrepancy between the degree of mixing of mine waters in the vertical and horizontal axes. In verticals, the variation in salinity is small; however, the results differ significantly between the right and left banks of the river. This may be due to the location of the section in question on the bend of the river and the related shift of the main stream towards the left bank.

The performed numerical analyses show that the brine mixing in the Odra River is effective for both small and large flow rates. In the worst case, the salt concentration is reduced from 40 to 0.632 kg m⁻³, which is well below the 1 kg m⁻³ limit.

Due to the low variability in the vertical axis, the currently used environmental monitoring system from KGHM can be used to measure the efficiency of the discharge installation. The interpretation of the measurement results should, however, take into account the differentiation of concentrations between the banks of the rivers.

eISSN:: 2083-831X
Sprache:: Englisch

Zeitrahmen der Veröffentlichung:: 4 Hefte pro Jahr
Fachgebiete der Zeitschrift:: Geowissenschaften, andere, Materialwissenschaft, Verbundwerkstoffe, Poröse Materialien, Physik, Mechanik und Fluiddynamik

Zeitschrift RSS Feed

Numerical analysis of the transport of brine in the Odra River downstream of a mine's discharge

Article Category: Original Study

Online veröffentlicht: 22. Dez. 2021

Seitenbereich: 366 - 379

Eingereicht: 23. Okt. 2021

Akzeptiert: 19. Nov. 2021

DOI: https://doi.org/10.2478/sgem-2021-0036

Schlüsselwörter
Water quality, mining, wastewater management, brine, dilution, numerical modelling, CFD

© 2021 Szymon Zieliński et al., published by Sciendo

This work is licensed under the Creative Commons Attribution 4.0 International License.

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6

Figure 7

Figure 8

Figure 9

Figure 10

Numerical analysis of the transport of brine in the Odra River downstream of a mine's discharge

Article Category: Original Study

Online veröffentlicht: 22. Dez. 2021

Seitenbereich: 366 - 379

Eingereicht: 23. Okt. 2021

Akzeptiert: 19. Nov. 2021

DOI: https://doi.org/10.2478/sgem-2021-0036

SchlüsselwörterWater quality, mining, wastewater management, brine, dilution, numerical modelling, CFD

© 2021 Szymon Zieliński et al., published by Sciendo

This work is licensed under the Creative Commons Attribution 4.0 International License.

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6

Figure 7

Figure 8

Figure 9

Figure 10

Schlüsselwörter
Water quality, mining, wastewater management, brine, dilution, numerical modelling, CFD