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

: The mining of underground deposits causes the inflow of water to workings and the necessity of pumping them to the surface. The mining plant of KGHM Polska Miedź S.A. extracts copper ore in plant branches with different hydrogeological conditions. The inflowing water into the workings is characterised by variable mineralisation, which depends on the location of the branch. In the south-western part of the deposit, a low-mineralised stream with a relatively high flow rate can be observed, while the outflow of highly saline waters occurs in the north-eastern branch. Despite the activities undertaken that aim at using the pumped-off mine waters industrially, it is necessary to deposit them into the Odra River. Reducing the environmental impact on the Odra River is one of KGHM’s goals, which is being implemented by stabilising its salt concentration at a safe level. The paper presents the results of a 3D simulation of brine plume propagation based on a numerical model of advection diffusion and turbulent flow. Bathymetric data from a section of the river approximately 500 m long and point data from an Odra water quality test were used to develop and validate the model. The paper discusses the types of factors that minimise the impact of brine discharge. The developed model will be used in the future to propose solutions that accelerate the mixing of mine waters with the waters of the Odra River.


Introduction
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).
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 -3 . 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 (a) (b) Figure 1: (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). 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 -3 ; however, fish species show resistance to salinity of several kg m -3 [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.

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

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.

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 3 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 3 s -1 , one dispersing nozzle works with a flow rate of 0.0625 m 3 s -1 .

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 3 s -1 and the average annual flow is 198 m 3 s -1 , 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 -3 . 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.   Figure 2: A diagram of the flow of the KGHM mine water [13].

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.

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.

Advection-diffusion of pollutants
The general form of the diffusion-advection equation is the position vector with Cartesian coordinates (the x axis is horizontal and consistent with the river flow direction and the z axis

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 timethe so-called Reynolds-averaged Navier-Stokes (RANS) methodor to appropriate filtering of the components of velocity fluctuationsthe 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]: where φ 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: where , p u 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: g ; g is the acceleration due to gravity) and ρ ′ ′ u 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 2 C ε . This model is given by the following system of equations: where is the tensor describing the production of turbulence, The RNG model contains an additional term in 2 C ε that is small in weakly strained turbulence and large in rapidly distorted flows [21].
where T ν is the eddy viscosity coefficient, ij S 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 k ε 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: 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.

Boundary conditions and the numerical program
The system of equations (1)-(3) must meet the following conditions: a. initial: - The concentration of dissolved salt in the river is negligibly low - The uniform velocity distribution at the inlet b. boundary: The normal component of the concentration gradient vanishes on the boundary - The discharge of the brine installation is constant - The concentration of the brine discharge is constant - The velocity vector of the brine seeping out of the river bottom is constant and normal where n is the normal vector to the boundary, d u is the velocity of the brine's outflow through the riverbed at the discharge site, d A ⊂ ∂Ω is the part of the bottom surface area through which the brine flows out and d q 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 0 t > 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. 0ˆo n F p S ⋅ = -T n n (11) ( ) 0 The 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 gridsthe 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.

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 2 (Fig. 6).
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, 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 0 c f = was assumed because the brine does not react chemically with the water and sediment in the river.

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.  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.
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.
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 3 s -1 was determined from the water gauge data. The brine discharge was 1.0 m 3 s -1 and its chloride concentration was 40 kg m -3 .
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.
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 3 s -1 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 3 s -1 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%.

Results of calculations
Calculations were performed for the following condition: six flow rates of the Odra River ranging from 41 to 318 m 3 s -1 , 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   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 -3 -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. 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 -3 .
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).

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. 1. 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. 2. 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. 3. 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. 4. 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 -3 , which is well below the 1 kg m -3 limit. 5. 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.