Pollutant Transport Via Sediments in Medium-Sized Meandering Rivers: The Example of Barium in Kupa River (Croatia)

The aim of our study is to expand the world-wide knowledge about pollutant transport in sediments of medium-sized meandering rivers and to explain barium distribution in Kupa River (Croatia), as it showed as an ideal “natural laboratory”. Sampling was performed twice at the same places with the time span of 15 years to see the changes in Ba concentrations after that period. Ba was analyzed using ICP-MS method on clay + silt sediment fraction (<63 µm). Modelling was performed in R environment (R version 3.5.1, 2018) and for solving the sediment transport equation of Ba in Kupa River a program based on �nite differences is constructed. After barite mine ceased operation in 1990, Ba is still entering Kupica and Kupa rivers, as karstic underground still contains large amounts of barite waste. Ba concentrations in Kupa River sediments are still increasing consequently along the whole river length, except in Kupica River and �rst downstream sample in Kupa River. It seems that in this part of river course erosion exceeded deposition processes. Also, there is an indication of changing the way of underground Ba transport. Erosion and deposition processes lead to formation of peak of Ba concentrations in the lower course of the Kupa River. Most important conclusion is that the mentioned peak was formed under one or two �ood episodes, when sediment eroded from the upstream locations and was brought and deposited in the lower course of the river. Our modelling predicts that Ba source will not cease completely during the next 20 years.


Introduction
A barium anomaly of a very high degree in sediments of Kupica and Kupa rivers was discovered in 2003 and results have been published afterwards (Frančišković-Bilinski 2005, 2006).Author of those studies found that Ba-contamination source originated from an abundant barite mine in Homer-Lokve area in Gorski Kotar region, Croatia.As a result of careless mine waste disposal in a channel leading to an abyss, the contamination penetrated through vulnerable karstic aquifers and appeared at Kupica River source.Biondić et al. (2006) performed a study in the Kupa River catchment.This study warned that the broader region of Kupa River basin, which is an area of deep karstic aquifers, is very vulnerable and that it should be protected as an important drinking water source.Downstream Kupica River source the barium anomaly in sediments reached more than 120 km along both (Kupica and Kupa) rivers.During 1990 barite mine in Homer was closed.But, this didn't stop in ow of Ba into Kupica and Kupa Rivers, as karstic underground still contained barite rich mine waste, which is slowly washing out of cavities in karstic underground.This was also con rmed with our data.Our study has a potential to be used in the future to study sediment transport processes in rivers since they act as an excellent "natural laboratory".In spite that sediment transport research, particularly in rivers with gravel bed lasts for decades, it is still a di cult issue with a lot of unknown facts (Ancey et al. 2014).
There are still plenty of open questions in this knowledge worldwide, also it is hard to nd appropriate river to follow sediment transport processes.Therefore, we repeated sampling on the same locations in 2018, exactly 15 years since the anomaly discovery, with the aim to see which processes happened during that time frame.
A uvio-geomorphological study of the Kupa River drainage basin was performed by Frančišković-Bilinski et al. (2012).Besides uvial geomorphology itself, this study investigated also ecological state of Kupa River.It was found that some parts of Kupa River, especially in its lower ow are in rather poor ecological state and under in ux of pollutants.For example, during the summer months some stretches of Kupa River are getting locally moribund and arrested within its course.Severe eutrophication with macroalgal blooms and generation of marsh vegetation are present in some places.All those natural characteristics and anthropogenic interventions in uence sediment transport in Kupa River and should be taken in account.
In the same quoted paper, Frančišković-Bilinski et al. (2012) calculated the sinuosity index for Kupa River, i.e thalweg length divided by the shortest distance from river source to mouth.It varies at different sectors of river with an overall sinuosity index of 2.20.Details and values of sinuosity index on some river stretches can be found in the quoted paper.According to its average sinuosity index, the Kupa River is classi ed into 'meandering' category of rivers (sinousity index=1.5-4).
Frančišković-Bilinski et al. ( 2014) investigated magnetic susceptibility caused by presence of magnetic spherules in the upper ow of sinking Dobra River, an important tributary of the Kupa River.They followed distribution of magnetic particles until the abyss in which Upper Dobra River sinks and have found that those particles were retained in the karstic underground, as they were not found in lower ow of Dobra River and thus didn't reach its con uence with Kupa.
Interdisciplinary study of Förstner (2004) showed that most important dynamic features of sedimentrelated processes in rivers are dramatic effects of stormwater ood events on particle transport, intense effects of sulphide oxidation during resuspension and biological accumulation.Potential release of toxic chemicals is also signi cant.Fine-grained sediments and suspended matter are crucial for pollutant transport, due to their large surface areas and high sorption capacities.Authors concluded that different underlying processes should be studied in rivers prior to modelling the overall effects of pollutants.
Suitable rivers for investigating sediment and contaminant transport processes are rare in the world and therefore there are not so many such studies.And because of lack of such "natural laboratories", some scientists perform laboratory experiments to test their models.
One of such studies was study of Yilmaz (2008), who performed experimental study of sediment transport in meandering channels measured in a laboratory ume and those measurements were evaluated and compared with previous model studies.They found that the values of the calculated sediment transport rates are generally overestimated.They used the Einstein's approach (Einstein 1942) to measure experimental sediment transport rates and sediment transport capacities.They have chosen Einstein's approach for bed load comparison because it has been partially tested.Park and Ahn (2019) performed study of primary ow patterns and of mixing, using experimental (construction of meandering channel) and numerical investigations.Obtained experimental results were used to develop a model.Authors nally concluded that same types of velocity characteristics and pollutant transport could be de ned with two-dimensional numerical models.
One of the rare studies performed on a real river was published by Nazir et al. (2016), who investigated sediment transport dynamic in a meandering uvial system of Chini River, Malaysia, with the aim to study the riverbed erosion and sedimentation processes.It was found that Chini River has sinuosity index of 1.98, what con rms instability of this river and that intense erosion processes are present on some parts of this river.Continuous erosion process in waterways has increased the riverbank instability due to the meandering factors, what is similar to situation in our study area -the Kupa River.
Ciszewski and Grygar (2016) found that metal remobilization can be more intensive in alluvia than in soils.Such processes are result of erosion of river banks and of prolonged oodplain inundation.They are also associated with reducing conditions, which are alternating with oxygen-driven processes happening during dry periods.Lauer and Parker (2008) performed a theoretical study on modeling framework, which was tested later on upper Clark Fork River, Montana, USA, which was exposed to mining activities.It was applied on sediment deposition, storage and evacuation in the oodplain of a meandering river and corresponding theory was developed.
The aim of the current work is to expand the world-wide knowledge about pollutant transport in sediments of medium-sized meandering rivers and to explain barium distribution in Kupa River (Croatia).

Study Area And Sampling
The Kupa River basin is presented in Figure 1.It is located in the western-central part of Croatia and its parts are also shared by two neighboring countries (Bosnia and Herzegovina and Slovenia).Kupa River is 294 km long.It is an important tributary of the Sava River, to which it in ows at Sisak (Croatia).The Sava River and its watershed occupy the south-western part of the Danube River watershed and in ows to the Danube River at Belgrade (Serbia).Kupa River and its whole watershed were intensively investigated during last two decades.Most important scienti c studies performed were those of geochemical assessment of river sediments of the whole watershed (Frančišković-Bilinski 2007) and of uvial geomorphology of the Kupa River (Frančišković-Bilinski et al. 2012).Most representative locations along both (Kupica and Kupa) rivers were selectd and sampling was performed there.Positions of those 18 sampling locations are presented in Figure 2 and details about sampling stations are given in Table 1.
Locations IŠ and 51 are situated on Kupica River, while all other locations are situated along Kupa River.From locations on Kupa River location 52 is situated upstream Kupica River in ow, while other locations are situated downstream its in ow.For sampling were chosen those locations where ne grained sediment accumulates along the river bank.On each sampling site a composite sample of about 1.5 kg was taken.Such sample was taken from at least three grab samples of active ne-grained surface sediment (0-5 cm deep), collected on different places within a surface of 5 m 2 , what decreased the effects of local variability.The same sampling procedure was performed in 2003 and 2018 and sampling was conducted on the same sampling points, except two points: due to border problems with neighbouring Slovenia sampling point 50 at Žaga (left bank) from the year 2003 was moved in 2018 to nearest accessible site from Croatian side -Čedanj, located on the right bank 5 km upstream and was marked as ČD.As Ba-anomaly arrived from Kupica River, in 2018 one additional point not sampled in 2003 (IŠ) was sampled on Kupica River near Iševnica, 1.8 km upstream sample 51.

Sample preparation
After sampling the sediments were dried in air at room temperature and then sieved through a 63µm sieve (Fritsch, Germany) to obtain ne fraction containing clay and silt (<63µm), on which study was performed.Approximately 0.1 g of representative powdered sediment samples was dissolved with 2.5 ml of suprapur nitric acid and 7.5 ml of puriss hydrochloric acid for half an hour at 1000 W in Anton Paar Multiwave 3000 Oven according to ISO 11466 norm and later diluted to 50 ml with deionised water.

ICP-MS analysis of elements in sediments
The elements contents were detected by inductively coupled plasma-mass spectrometry (Agilent 8900 ICP-MS Triple Quad), with solution of 30 µg L −1 Ge, Tb, In and Re as internal standard according to HRN EN ISO 17294-1 and HRN EN ISO 17294-2 norms.All measurements were performed in ve replicas.The precision, evaluated directly as the relative standard deviation (RSD), was better than 10 % for all determined elements.The accuracy of the ICP-MS analytical procedure and method quality control was performed by the analysis of the elements of interest in standard reference material (RTC, Trace elements on fresh water sediment, catalog number: CNS392-050), which were analyzed at the beginning and after analyzing each series of samples.A generally good agreement within 15% was observed between our data and the certi ed values.

Mathematical modelling
Law discovery of decrease in elevation above the sea level versus distance from the spring of the river Kupa, curve tting and a transport law using modelling were performed in R environment (R version 3.5.

River elevation versus distance
Longitudinal water elevation pro le above the sea level of Kupa River from its source until in ow to Sava River is presented in Figure 3.This is an important feature for contaminant transport modelling in sediments.As one can see, a steeper slope prevails up to approximately 130 km downstream from the river source, the distance is found downstream the Ozalj hydropowerplant.Roughly, this point is approximately at the half of the total length of the river.From its source until its in ow to Sava River at Sisak, the height of Kupa River decreased for 217 m (from 313 m to 96 m) above the sea level.

Sinuosity index
Sinuosity index was calculated and discussed in details in the earlier research by Frančišković-Bilinski et al. (2012).This index is thalweg length divided by the shortest distance from river source to its mouth, but it can also be calculated for particular sections of a river.According to results from the mentioned paper, sinuosity index for the whole length of the Kupa River is 2.20.When it is divided into sections, in the uppermost part of the river (from source to Žaga) it is 1.96.In the middle part of the river index is 2.42, while in the downstream part of the river from Karlovac to the con uence with Sava at Sisak it is 2.10.The value of sinuosity index mathematically describes the type of the stream.With its average sinuosity index of 2.20, the Kupa River belongs to the category of "meandering rivers", which have sinuosity index between 1.5 and 4.0.
Barium concentrations in the Kupa River sediments in 2003 and 2018 In Table 2 are presented all available data of Ba concentrations in sediments of Kupa River listed downstream the Kupa River from its source under Risnjak Mountain in Gorski Kotar until its in ow to Sava River at the town of Sisak.For each point the exact distance from the river source in km is given.Ba-concentration ratios for years 2018 vs. 2003 In Table 3 data are presented for all measured concentrations of Ba in sediments, including those from Kupica River (IŠ and 51), which is the source of Ba-contamination of Kupa River.Calculated ratios of Baconcentrations in sediments for years 2018 vs. 2003 are also presented in the table, as they are illustrating how many times Ba-concentrations on some point increased or decreased.One could notice that Ba-concentration on sampling location 51 (Kupica River just before its con uence with Kupa) decreased almost 5 times during this 15-year period, so contamination input decreased a lot during that period.This is due to the fact that since 1990 there is no new active Ba-input into the karstic background, as Ba-mine in Homer ceased operation.Obviously, a process of rinsing of karstic underground which was enriched with Ba-ore waste continued but less and less Ba is reaching Kupica source.Consequently, Baconcentrations in Kupica River started to decrease and therefore total input of Ba-contamination from Kupica to Kupa River dropped signi cantly.The same process is observed on the rst downstream sampling location on Kupa River (50 and substitute location ČD).The most signi cant increase of Baconcentrations in sediments in the observed period of 15 years happened in the lower ow of Kupa River.On that stretch of the river in 2003 increased Ba-concentrations were not observed and they were approximately at the background level characteristic for this region (<100 ppm) while 15 years later Baanomaly reached until the end of Kupa River and its in ow to Sava River.

Mathematical modelling
Elevation versus distance from the spring of Kupa In order to describe a decrease of elevation, the graph of elevation above sea level versus distance from the spring is examined (Figure 3).In the rst attempt we tried to t an exponential function of the form: E = E o exp (-a x) where E o = 313 (m) is the elevation at the spring, a is a free parameter to be determined and x is the distance from the spring in (km).However, the best t shows the function to have insu ciently steep slope close to the spring and insu ciently gentle slope further from the source.This indicates that one should not take the distance as an independent variable but a function of x which grows slower.
Consequently, a law of the type: gives the best t for a = 0.0767±0.0015(with the RMSE = 11.15 on 30 df).The t to data is shown in Figure 5.
Further comparison of the function to the data in Figure 5 reveals that this law is completely satisfactory until 130 km when a sudden decrease in elevation occurs.Subsequently, elevation decrease of river Kupa has more gentle slope as it is laid on the Pannonian atland.

River slope
From the data used to obtain the Figure 5 one can compute the slope along the river (Figure 6) as the slope is the derivative of water elevation pro le.
From Figure 6 it follows that the negative slope starts as over six times greater in absolute value close to the Kupa source than in the region where the river meets the Pannonian Plain.

Barium in sediments
The concentration of Ba in Kupa River measured in 2003 y. may be tted by an exponential law: In the exponential function Ba 0 = 1060 (ppm) is the concentration value at the entrance of Kupica to the river Kupa, and x is the distance from the entrance.The best t is obtained for b = 0.0167 ± 0.0014 (with RMSE = 83 on 12 df).Let us also note that the law behind the exponential function is: d(Ba)/dx = -b Ba with the initial condition equal to Ba 0 .
In order to explain how did the Ba concentration arrive to the distribution shown in Figure 7 during 70 y since the barite mine in Homer started operating and Ba rich particles from the mine waste penetrated the vulnerable karstic underground and through it reached the Kupica River source, a sediment transport model is considered.
The transport model for the concentration of Ba (C) in Kupa is given by: with the boundary condition Ba| x=0 = Ba 0 and the initial condition Ba| t=0 = 0, for every x > 0.
The rst term in (3) accounts for diffusion of particles while the second term accounts for advection downstream from the entrance of Kupica into Kupa.Since we have no data from which to determine the diffusion coe cient and furthermore, since the turbulent diffusion coe cient is from10 −7 to 10 −5 cm 2 /s (Wörman 1998), this term may be safely neglected in considered longitudinal transport.
What remains is the advection of sediment downstream.There are at least two processes of importance for movement of sediment downstream.The rst is that during higher precipitation, the horizontal velocity of water increases, increasing vertical turbulent water movement and entrainment of sediment into the water column.This increases advection of sediment downstream.The second process is deposition of sediment at places further downstream where the slope and velocity of water are smaller.
The advection model is discretized with a step ∆x = 1 (km) and run for 70 y with a ∆t = 1 (y) given the fact that once in a year we have an exceptionally high river ow due to the snow melt, and a constant boundary concentration Ba 0 = 1060 (ppm) at the entrance of Kupica into Kupa where x = 0.The initial condition is set to Ba = 0 for all x > 0. Since no velocity of sediment was available, we used the discrete law which states: and adjusted c = 0.75 and d = 0.76 for the best t to data measured in 2003 y.
The resulting simulation is shown in Figure 8.
From Figure 8 it is seen that the initial concentration at t = 0 is located only at the entrance of Kupica into Kupa.Through the 70 years it has gradually moved downstream.Finally, it appears as an exponential distribution of Ba in sediments reaching over 200 km downstream Kupa River.
It is important to point out that an attempt to use a = b failed to t the data from 2003 y for any value of a since it produced too steep a slope of Ba concentration.This strengthens the discovery that indeed at any distance downstream from the entrance of Kupica into Kupa, 75% of the sediment containing Ba is added per km per year from upstream and 76 % of the sediment found at the location is transported downstream.
Veri cation using the data of 2018 y.
Using the same transport parameters a and b, the model was run until 2018 y.The comparison until 130 km downstream the entrance of Kupica into Kupa reveals a good agreement to the data (Figure 9).For the explanation of the "bump" in the distribution occurring from 2003 to 2018 y downstream from 130 km shown in Figure 9, one has to invoke at most two new factors to the model.The rst is that downstream beyond 130 km, the elevation slope is signi cantly smaller so the deposition potential is greater.The second factor is that in this time interval, one or more violent precipitation events may have happened which caused entrainment of heavier sediment into the water column at the stretch of the river prior to 130 km (where the water velocity is greater) and deposition downstream beyond 130 km where the river water ows practically along a at plain causing signi cantly slower movement of water and hence greater deposition potential of sediment from the water column.
In the model only the rst factor has been simulated as a gradual slightly greater deposition then loss past 130 km during 15 y, i.e. from 2003 to 2018 y.
The Figure 9 shows the result of two pieces of simulation.The last exponential function (the envelope) until 130 km represents the result that the distribution would have had if the same simulation, which is shown on Figure 8, is continued to 2018 y i.e. 15 year longer.
The concentration "bump" is obtained by simulation of the part beyond 130 km with c = 0.75 and b = 0.74.Another words, with the same import from upstream as in the period until 2003 y but with only slightly greater deposition, or equivalently, with smaller release of sediment downstream from 2003 until 2018 y.
However, as stated above, the same shape of the distribution beyond 130 km downstream the source may be explained by one or more violent precipitation events which have moved sediment from the upper part of Kupa over the sudden change of elevation slope into the at land.Indeed, the examination of historical data on violent precipitation events in the last 100 y. reveals the highest ood during 2015 y and the second highest ood during 2014 y.The oods in 2014 and 2015 alone could have moved the sediment downstream in agreement with the data in Figure 4.
Finally, high concentration of Ba in sediments measured in Kupa 200 m upstream from the entrance of Kupica (station 52) during 2003 and 2018 may be explained by diffusion of sediments upstream.As for the explanation of the two somewhat higher values above the background level of Ba closer to the source of Kupa, the likely mechanism may be that during an event of intensive precipitation in the Kupa watershed, Ba is transported into Kupa through the porous karstic bedrock from the mine deposit uphill.
The same mechanism was found to be responsible for occurrence of Pb in the beach sediments which were transported through porous karstic rock from uphill lead processing workshop into the Punat Bay of the Adriatic Sea (Legović et al. 1990).
From measurement and modelling results it may be predicted what will happen to Ba concentration in Kupa during the next 20 y.
First, the Ba source will not cease any time soon.It may decrease further but it will not completely vanish in the near future.Given the fact that barite mine was closed in 1990, we see Ba concentration upstream of Kupica entrance into Kupa is still increasing in 2018, 38 years later.This means that Ba entrance into Kupa through a porous karstic rock is still higher than erosion downstream.Likewise, Ba concentration downstream Kupica entrance into Kupa is still increasing testifying again that entrance is higher than erosion downstream.However, there is a sign that this process started decreasing as Ba concentration in sediment at the location 51 (Kupica River just before its in ow into Kupa River) was in the year 2018 about 4 times lower than in 2003 (see Table 3).Also, at the rst downstream location from the Kupica con uence (50) Ba concentration in sediment was almost 30% lower in the year 2018 than it was in 2003 (see Table 3).This is obviously a sign that parts of karstic underground leading towards Kupica River source are gradually being washed out.If this process persists, in the future we might expect further decreasing of Ba concentrations in the upper course of the Kupa River, while in the lower course this will not happen so soon.Finally, a peak of Ba concentration downstream from 130 km has a potential to increase further and in part be transported further downstream.Hence, by the end of next 20 y we do not expect that Ba concentration will decrease in the lower part of the Kupa River, especially downstream of Ozalj powerplant.To the contrary, we expect an increase throughout Kupa, except on part of its upper ow just downstream of Kupica River in ow.

Conclusions
Results obtained within the current paper led to the following conclusions: After barite mine in Homer ceased operation in 1990, Ba is still entering Kupica and Kupa rivers, as karstic underground still contains large amounts of barite waste.Our results have shown that between the years 2003 and 2018 Ba concentrations in sediments are still increasing consequently along the whole river length, with the exception of sampling point in Kupica River before its in ow to Kupa and rst sampling point in Kupa River downstream the con uence with Kupica.
The signi cant decrease of Ba content in Kupica River sediment and rst downstream sediment sample in Kupa River are signs that a reversal started due to the process in which parts of karstic underground leading towards Kupica River source are gradually being washed out of barite waste.Now in this part of river course erosion exceeded deposition processes.
An important nding is that between the years 2003 and 2018 Ba concentrations in sediments signi cantly increased at the sampling point upstream the Kupica River in ow.This could mean that present Ba-waste in the karstic underground found another way.Except owing to Kupica River source, now part of it enters Kupa directly at the stretch of river upstream Kupica in ow.But, obviously these quantities entering Kupa are now not so large as before.This is indicated by the fact that on the rst downstream location Ba concentrations started to increase i.e. that erosion now exceeds deposition.
About 150 km downstream from Kupica River in ow into Kupa there is a peak of Ba concentrations.
We assume that it was formed under one or two ood episodes, which happened between 2003 and 2018, namely 2014 and 2015.During such events of extremely high water discharge, sediment may have eroded from the upstream locations, brought downstream and nally deposited in the lower course of the river.
From measurement and modelling results it may be predicted what will happen to Ba concentration in Kupa River sediments during the next 20 years.We predict that Ba source will not cease completely during that period.Although, in most part of the river length, except the uppermost part downstream Kupica River in ow, Ba deposition will still be smaller than erosion.However, in the lower-most part of Kupa, concentrations will be increasing since there deposition is greater than erosion.
Further research of this interesting phenomenon is advised.Especially, it will be important to perform new sampling at the same locations in the future, e.g. in 5, 10 or 15 years after the last sampling campaign from the year 2018.First, this would enable computing an estimate of the rate of decrease of Ba in ow into Kupa.Second, it would enable answering te question of whether Ba will eventually enter Sava river or it will be accumulated, buried and compacti ed further into the present peak section of the Kupa.

Declarations
Funding: Research was partially funded by Ministry of Science and Education of Republic of Croatia (Grant No: 0411FI18).
Con icts of interest/Competing interests: We state that there is no con ict of interest or competing interests. Figures

Figure 1 Position
Figure 1

Figure 2 Position
Figure 2

Figure 5 The
Figure 5

Figure
Figure

Figure
Figure

Table 1
Details about sampling locations: locality, watercourse, to which watercourse it ows, geographical coordinates, elevation and sampling date

Table 2
Frančišković-Bilinski (2006)iments of Kupa River listed downstream the river for both years(2003 and  2018).The distance from the Kupa River source is given for each sampling point.Data for the year 2003 are taken from previous study ofFrančišković-Bilinski (2006), while data for the year 2018 are original data obtained for the current paper.n.a.= not available In Figure4are presented Ba concentrations in sediments (expressed in ppm) vs. distance from the source measured in the both studied years.As it could be seen from the gure, during those 15 years Baconcentrations increased signi cantly in the lowest course of the Kupa River and reached until its con uence with Sava River.Contrary to that, in the year 2003 in sediments of the last 100 km of the river course Ba-concentrations in sediments were not visible and they were approximately at characteristic background level for Ba in the studied region (<100 ppm).
* Data for the year 2003 are taken from previous study of Frančišković-Bilinski (2006), while data for the year 2018 are original data obtained for the current paper.n.a.= not available *

Table 3
Measured concentrations of all studied sediment samples, including those from Kupica River and their ratios for years 2018 vs. 2003