Numerical model of stratified flow – case study of the Neretva riverbed salination ( 2004 )

Two natural directions of saline sea water intrusion have been registered in the lower reaches of the Neretva River: the first direction of intrusion is through the Neretva riverbed, and the second one through deeper underground layers. The numerical model of stratified flow is described, and its application in the analysis of the Neretva River salination via the first direction of intrusion is presented. The model is calibrated based on measurements conducted in 2004. It was established that saline water appears in Metković at fresh water flow rates of less than 180 m3/s, while saline water is fully blocked out from the riverbed when water flow exceeds 500 m3/s.


Introduction
The area of the lower Neretva is a specific natural geographic area, surrounded by the karst rock massif.This is an area abounding in water and characterized by a mild Mediterranean climate, which is why an intensive agricultural production has been developed on meliorated land.However, the production has recently become endangered due to occasional or continuous salination of surface water and groundwater in the area downstream of Metković [1,2,3,4].This results in lower agricultural yields, and in destruction of fertile agricultural soil.
Previous studies have shown that there are two natural directions (routes) of salt intrusion in this area [4].The first route is via the Neretva riverbed, and the second one originates from deep layers where groundwater is heavily salinated.The movement of salt in alluvial aquifers depends on climatic, hydrological and hydrogeological conditions, but also on the use of land improvement and irrigation systems.The development and use of special numerical models, and appropriate "in situ" measurements, are needed to model the underground flow.The phenomenon of salination of watercourses and river mouths (estuaries and deltas) has been studied by many experts and scientists all over the world.This is not done solely for environmental reasons, but also because of economic considerations, as these areas are often important water resources, and hence of high significance for water supply, agricultural production, and industry.The interaction and mixing of salt sea water with the water flowing from the upstream part of the basin takes place at estuaries.The entire phenomenon is characterized by highly complex dynamic and stochastic processes.In order to comprehend and describe these processes, and apply appropriate methods to counter saltwater intrusion, modern methods must be used and, in this respect, a special use is made of appropriate mathematical (numerical) models that can fully describe surface flow of multi-phase fluid.The application of such models requires collection and processing of meteorological, hydrological and hydrogeological data, which characterize the situation in the watercourse and estuary.The study of salt and fresh water interaction gained in intensity about five decades ago.At that time, this activity mostly consisted of laboratory modelling [5].However, with the development of electronic computers, the use of numerical models gained in momentum.Owing to development of modern measuring equipment, extensive "in situ" measurements can now be conducted on rivers.Excellent results can be achieved by combining such measurements with numerical modelling [6,7].In Croatia, three biggest rivers that flow into the Adriatic, i.e.Krka [8], Zrmanja [9] and Neretva, are salinated in their lower reaches.The most comprehensive measurements were carried out at the Neretva River as, in its lower course, this river has water facilities forming a system that is of high economic significance for the Republic of Croatia.It was established that, in the period of low freshwater inflow, sea water enters into the riverbed through the Neretva river mouth, and forms a flow of fresh water above the sea water wedge (the so called "salt wedge").In the summer period, the salt wedge reaches the town of Gabela, which is situated 25 km away from the mouth.In the rainy period, when inflows from the basin are higher, the wedge is pushed downstream, until the point when it is completely removed from the riverbed.Salt wedge occurrences are found in sea bound rivers characterized by a relatively small tidal amplitude [10].According to Kurup [11], these are the seas with less than 2 m in amplitude.The corresponding river mouths are also called microtidal estuaries [12].The position of the salt wedge changes during the year, which is dominantly dependent on the change in fresh water inflow [13,14], but also on sea tides, bathymetry and wind activity.Maximum sea water wedge intrusion into rivers may vary from several kilometres to several hundred kilometres.The most significant examples of this type of river estuaries, often mentioned in literature, are the Ebro River (Spain) and the Rhône (France) in the Mediterranean [13,15], and the Mississippi River in the United States, where the maximum length of the salt wedge, amounting to as much as 190 km, was recorded [16].
In case of salt wedge, there is little mixing of salt and fresh water.In fact, this mixing is exclusively operated as a vertical one-way advection of saline water into the upper layer [17,18,19].There is however a very narrow transition zone between the two layers.Thus, there is a sharp interface between the layers, which is a halocline in the estuary (sudden change in salinity), and also a pycnocline (sudden change in density by depth), and often a thermocline (pronounced change in temperature).A number of numerical models are used to solve numerous water management tasks relating to salt (sea) water intrusion in coastal streams.Many software packages for modelling hydrodynamic and transport processes are available on the market.It should however be noted that some can not be used for solving stratified flows.Also, software packages are often very expensive and therefore inaccessible.Measurements were carried out in several verticals of the riverbed cross-section in the scope of initial studies of the Neretva riverbed salination level [4].The results show that salinity distribution is uniform over the entire width of cross section.For that reason, subsequent measurements were made in one vertical depth of flow, in the middle of the river, and one-dimensional model for flow simulation was selected.Basic features of a one-dimensional numerical model of salt and fresh water stratified flow, based on theoretical analyzes conducted by several authors over the past five decades [20,21,22,23], are presented.The application of this Numerical model of stratified flow -case study of the Neretva riverbed salination (2004) model in the actual Neretva riverbed, under the influence of time-dependent boundary conditions, is also depicted.The calibration of parameters and model verifications were made on the basis of measurements carried out in the Neretva River in 2004.The results show that salt wedge intrusion into the riverbed, under the influence of tides and fresh water inflow from the upstream part of the basin, can accurately be simulated by means of a one-dimensional model.Although this model analyzes hydrodynamic flow parameters only, it can be extended to also include biological and environmental components as well [14,24].
Several authors present the development and use of 2-D numerical models with hydrodynamic equations, averaged across the riverbed width [6,11].As geological, hydrological and hydrodynamic relations in the lower course of the Neretva are highly complex, it is very difficult to protect soil and surface watercourses from salination.In fact, this area is characterized by overlapping of different phenomena including water flow in karst, and the flow and mixing of fresh and salt water in an inhomogeneous alluvial environment, and in surface watercourses [25].The preservation of natural features, and biological and landscape diversity in the lower course of the Neretva, is not only significant for the Republic of Croatia, but has been recognised on an international level as well.Thus the Neretva lower course has been included in the Ramsar list of wetlands, and in the program of important ornithological areas in Europe.A special attention should be paid to these facts during definition and implementation of measures for protecting soil and water against salination.

Site description
The Neretva River takes is source in Bosnia and Herzegovina to the southeast of Zelengora mountain, at 1095 m a.s.l.It is 215 km long, with only 22 km in the Republic of Croatia, at the lowest reaches of the river.At its mouth the river empties into the Adriatic Sea (Figure 1).The area of the Neretva River drainage basin is estimated at about 10,500 square kilometres, out of which about 280 square kilometres are situated in the Republic of Croatia [25].In its upper and middle course the Neretva is a typical mountain river running through narrow steeply-sloped valleys.To the downstream of Počitelj (km 35+800), the Neretva leave the canyon part of its course and runs with a meandering mainstream through a valley (Figure 2).This area has been formed over an extended geological time through deposition of eroded material originating from all parts of the drainage basin.A complex hydrographical network, with the Neretva delta and numerous smaller river Igor Ljubenkov, Mijo Vranješ branches, has developed in the alluvial area of some 19,000 hectares.The situation remained unchanged until the mid twentieth century when the soil amelioration campaign was initiated.Numerous flood protection dykes, dams and other water engineering facilities were built, which all greatly altered natural features at the delta.Presently, the river course to the downstream of Metković is fully channelled.On the sides of the main channel of this river, we can still see numerous traces of former delta -with many abandoned meanders and river branches.[25].Agricultural production is highly developed in many parts of the valley, while the remaining areas are occupied by wetlands and lakes.In the lower course of the river, the channel route used to vary considerably during the year, depending on hydrological conditions in the entire basin.The pronouncedly irregular natural hydrological regimen of the Neretva River (Figure 1) has been greatly regulated after construction of the hydropower plant system in the Neretva basin, with storage basins of Jablanica, Rama, Grabovica, Salakovac and Mostar, and with the Čapljina hydropower plant's balancing reservoirs in Popovo polje and in Svitava area.The greatest annual discharge rates in the lower course of the Neretva are registered in winter period, from November to April, while small discharge rates are typically registered from June to October.
The Neretva discharge rate in Metković (at KM 20+875) is primarily dependent on operation of the Mostar hydropower plant.The data from the water level recorder in Žitomislić (KM 47+000) are used for determining quantities of water coming from this hydropower plant.To the downstream of Žitomislić and all the way to Metković, the Neretva has three greater tributaries, namely Trebižat (right-side), Bregava and Krupa (left-side).Quantities of water these tributaries bring to the Neretva are measured via the existing gauging stations located near their mouths.In the rainy period, discharge rates of these tributaries may be significant: up to several tens of cubic meters per second.In dry season, discharge rates amount to several m 3 /s only.In the zone from Metković to the mouth, the discharge rate of the Neretva varies considerably and can not easily be measured, as this part of the watercourse is greatly influenced by the sea (backwater).The maximum discharge rate in Metković, measured before construction of hydropower plants on the Neretva River, amounts to 2180 m 3 /s.The mean perennial discharge rate of the Neretva at Metković cross-section is estimated at 355 m 3 /s [25].On the stretch from Metković to Opuzen, the Neretva River has only one right-side tributary called Norin.It is difficult to accurately define discharge rates of the Norin at the point where it joins the Neretva, as this rivulet is constantly influenced by the Neretva water levels.The sluice gate which controls discharge of water from the main canal of the Neretva into the Mala Neretva is of crucial significance for the distribution of water in Opuzen.That gate protects Mala Neretva against flooding and salt intrusion, and is therefore closed most of the time.In the dry season the gate prevents the Neretva salt wedge from penetrating into the Mala Neretva, and in the rainy season it does not allow the flood wave from the main canal of the Neretva river to penetrate into the Mala Neretva.The Mala Neretva riverbed is also protected at the downstream side by the gate, which prevents direct penetration of sea water into the riverbed.A large number of springs with very different qualitative and quantitative hydrological characteristics can be found in the wider area of the lower Neretva basin.However, in terms of modelling flow rate in the Neretva riverbed, it is important to note that in the dry part of the year the abundance of all springs is relatively small, and so the impact of these springs can be neglected.
Water levels at the Neretva estuary depend on flow rate of water coming from the upstream part of the basin, and on sea-level fluctuations the influence of which is felt even upstream from Metković (about 25 km from the mouth) at low water levels in the Neretva.Water level recorders installed in Metković and Opuzen are used for monitoring water levels in the Neretva River.Fluctuations of the Adriatic Sea levels at the Neretva mouth are caused by tides, but also by wind action (sea level rises along the coastal area).The tides are of mixed type, i.e. they are characterized by semi-diurnal and diurnal rhythm.The usual diurnal amplitude is 0.5 m, and it varies from 0.1 m a.s.l. to 0.6 m a.s.l.Data from the tide gauge installed near the Mala Neretva mouth are used for this area.

Measurements in the Neretva
Organized measurements of salinity in the Neretva riverbed have been conducted by experts from the Faculty of Civil Engineering and Architecture -Split since 1997.The purchase of a highly sophisticated Hydrolab miniature probe in 2002 has enabled high-quality measurements of the following waterrelated parameters: temperature, electrical conductivity, salinity, pH, and dissolved oxygen content.The speed is measured with the current meter.A detailed presentation of measurements and relevant results is given in [4].Studies have shown that the flow in the Neretva is stratified, i.e that the fresh water flows above the salt water wedge, which penetrates deep upstream along the riverbed.Due to influence of fresh water from the upstream part of the basin, and the influence of sea tides, the wedge sometimes moves upstream and sometimes downstream.The sea usually starts to penetrate into the Neretva riverbed in springtime, when freshwater inflow from the upper part of the basin is smaller.In summer months, the sea water wedge penetrates through the Neretva mouth, and continues upstream all the way to Gabela (about 25 km away from the mouth).In late summer or autumn, when the fresh water inflow becomes greater, the wedge is pushed downstream, until it is completely pushed out of the riverbed.
In the previous period, measurements were conducted on some ten sites along the riverbed, from the point of discharge into the Adriatic Sea to the bridge in Metković (on the stretch 21 km in length).Measurements are most frequent at the water level gauging site in Opuzen (km 11+880).Measurements from 2004 are used in this paper for calibration of model parameters, and for verification purposes.In that year, the measurements were the most detailed, and were conducted 8 times (days) in the period from April to October.The total of 30 measurements were made in 6 profiles, from S1 to S6, cf. Figure 2.These measurement points are located at 0+000 (mouth), 0+500, 1+000, 3+621.5 (Rogotin), 11+880 (Opuzen water level recorder), and 20+875.5 (Metković Bridge).Salinity levels measured on 7 July 2004 and 13 September 2004 at four typical cross sections along the Neretva River, i.e. at the estuary, and in Rogotin, Opuzen, and Metković, are presented in Figure 3.At the estuary, brackish water of very low salinity level (<5g/l) is present in the surface layer of water, which is 1 to 1.5 m in thickness.It can be seen at other upstream cross sections that two layers of constant salinity can clearly be differentiated in the riverbed.At the surface, the water is slightly brackish, while salt sea water (about 38 g/l) can be noticed in the nearbottom layer.The transition zone between the top fresh water layer and the bottom layer is very narrow (about 0.5 m).This is the main feature of river estuary of stratified type, and so the existence of a sharp interface between layers can be adopted for modeling purposes.The surface layer thickness gradually increases in the upstream direction, and amounts to some 4.5 m (7 July) and 2.5-3.0 m (13 September) in Metković.The gradient of boundary surface of layers (halocline) can be noticed in the longitudinal direction.
Velocities are the highest at water surface, and are gradually reducing along the depth [26].In the near bottom layer, velocities have positive and negative values, depending on boundary conditions (change in fresh water inflow, and influence of tides).

Mathematical formulation of stratified flow
Several authors [20,21,22,23,26] have presented mathematical models of stratified flow, using various assumptions.Therefore, only some basic equations, and assumptions adopted for simulation of one-dimensional non-stationary stratified flow at river estuary, are presented below.The following assumptions were adopted for derivation of equations: -homogenous, isotropic and incompressible fluid, -one-dimensional flow, -steady flow regime with slight changes in flow, -the riverbed is of complex cross-section, described with polygon, and so model members are integrated for the entire cross section, -stratified flow with sharp interface, without mixing, -fluid density in each layer is constant and unchanging, -mean flow velocity is taken for each layer in the riverbed cross section, -hydrostatic distribution of pressure is adopted along the depth, -water table, and boundary surface between two layers, are horizontal in the cross-section, -it is taken that he inclination of the riverbed bottom is small, and so it can be assumed that the angle cos is approximately one.
The upper layer (Layer A) is located between the water table and the interface and represents fresh water, and the lower layer (Layer B) is located between the interface and the riverbed bottom and represents sea water (Figure 4).-continuity equation for Layer A -continuity equation for Layer B -dynamical equation for Layer A τ τ ρ ρ -dynamical equation for Layer B Index A denotes the upper layer, and index B denotes the lower layer.Variables A and O are areas of flow and wetted perimeter.B S is the interface width, g is gravity, ρ is density, and h is the upper layer thickness.Tangential stress has two components -riverbed shear (τ 0 ) and interfacial shear (τ S ).
As equations ( 1), ( 2), ( 3) and ( 4) are hyperbolic and do not have an analytical solution, the Finite Elements Method (FEM) was applied to solve the problem.This is one of approximate numerical methods that are often used for solving engineering tasks [27].
The discretization of continuum on finite elements, approximation of solution by elements, and selection of interpolation functions, were conducted in a way most suited to the nature of the problem.Two-node linear elements were used, with shape functions in form of first-degree Lagrange polynomials (linear functions) [28].A strong formulation of mathematical model is applied, using Galerkin methods.The equations of continuity and dynamic equations for each element were integrated in space [x 1 , x 2 ] and time [t 1 , t 2 ].According to numerical integration rules, they were reduced to nonlinear algebraic equations.Therefore, four algebraic equations are defined on one element.There are 4xM such equations for the entire system with the total of M elements.The above equations, obtained by applying the mass and momentum conservation law, are set on an isolated part of the riverbed, i.e. on a finite element, and are therefore called Element Equations.
For the system in which we have M elements and N nodes, 4M discharges and 2N water levels are calculated, which gives the total of 4M+2N unknowns.Additional equations that close the system are obtained by setting the equation of continuity in nodes (the sum of flows allocated to a node should be zero).For each node "i" it is assumed that the sum of all Layer A GRAĐEVINAR 64 (2012) 2, 101-113 flows, which enter and exit the node, should be zero.The same applies to Layer B. These equations are called nodal equations.
The following is valid for node "i" and Layer A: where m is the total number of elements in node "i", and i AV Q is the "external" discharge of freshwater in node "i".Parameter p = 1 for the upstream element node, p = 2 for the downstream element node.The similar is valid for Layer B: where i BV Q is the "external" discharge of salt water in node "i".
Nodal equations in algebraic form are obtained through application of numerical integration in time.The system of 4M+2N nonlinear equations (4M element equations and 2N nodal equations) is obtained through element and nodal equations ( 5) and ( 6): x i is the vector of requested solution.As the previous system ( 7) is nonlinear, the Newton -Raphson iterative method is applied for its resolution, and so the problem is reduced to the issue of solving a linear system of 4M + 2N equations [27,28].An unambiguous solution of unknown functions h and Q is obtained in finite-elements nodes for given initial and boundary conditions.
For the initial state of the system, it is necessary to specify the water level value for Layer A (h A ) and Layer B (h B ) in all nodes.In addition, four flows are set for each element, and this on the upstream element node for Layer A (Q Aj1 ) and Layer (Q Bj1 ) and on the downstream node for Layer A (Q Aj2 ) and Layer B (Q Bj2 ).Boundary conditions are entered in form of hydrographs (Q-t) or water level graphs (h-t) in boundary nodes of the system (initial and final node of the topological scheme).The hydrograph is most often prescribed at the downstream boundary of the system, while discharge is prescribed at the upstream boundary.to the bridge in Metković.These cross sections were spaced at about 250 m intervals.In addition, ten additional fictitious cross sections, situated in the sea right in front of the river mouth, were taken into account, in order to annul possible errors in estimating the downstream boundary condition.
The selected cross sections also represent finite element nodes.The topological scheme is formed by numbering the nodes (profiles) and elements (Figure 5).

Boundary conditions
The longitudinal cross section of the bed, with marked boundaries of the system, is shown in Figure 6.The downstream boundary of the model is placed in the sea, 2.5 km to the downstream of the river mouth.The upstream boundary is located in Metković at the site of the old demolished bridge.The remaining parts of the bridge, which have not been removed, create a threshold at the bottom of the bed, which prevents stronger upstream penetration of sea water [2].In addition to this artificial barrier, the riverbed is elevated a few kilometres upstream from Metković, and this natural barrier also prevents further intrusion of seawater.

downstream water levels (h A and h B )
Although there is a water level recorder in Metković where water levels are continuously measured, this profile is subjected to sea storm surges, which is why the flow can not be accurately determined.Therefore, the discharge is taken from the first available upstream station, namely from the water level recorder in Žitomislić (km 47 + 000).During the summer, the fresh water discharge in Metković is primarily dependent on the operation of the Mostar hydropower plant, because lateral flows into the Neretva are very small.Thus, the upstream boundary condition is Q A RUB = Q ŽITOMISLIĆ .Therefore, the boundary discharge is the function of time which is given in form of a polynomial in the most upstream node (node No. 93).Hourly values of discharge are shown in Figure 7a.For the lower salt water layer the boundary discharge is constantly zero Q B RUB = 0, which in physical terms means that an impermeable barrier can conditionally be adopted for the lower layer at the upstream boundary, or that there are no upstream or downstream fluctuations.At the river mouth, freshwater generally discharges into the sea in surface layer that spreads like a fan.This is a phenomenon of multi-phase fluid flow, which will not be analysed here, and so the thickness of the surface layer of fresh water at the downstream boundary needs to be estimated (h A RUBh B RUB ).Due to above mentioned considerations, cross sections (profiles) situated in the sea were taken into account, in order to reduce the influence of the downstream boundary condition assessment error on the solution of the problem, i.e. on definition of water levels and discharges in the riverbed.Tide gauge measurements (hourly values) taken at the Mala Neretva estuary h A RUB = h MAREOGRAF (Figure 7.b) were accepted at downstream boundary of the system as Layer A water level.Boundary water level of saline sea water (h B RUB ) was estimated based on a set of measurements as a function of the Layer A water level h A RUB and boundary fresh water flow Q A RUB in such a way that the fresh water layer thickness at the downstream boundary h A RUB -h B RUB is variable, and ranges from 0,6 m at small discharges to 2,0 m at greater discharges (Figure 7b).

Description of simulation
The dynamics of the salt wedge was simulated during the period from 10th June 2004 (10:00) to 31 st Oct 2004 (00:00).
In order to choose the best combination of model parameters (integration time interval Δt, interpolation parameter, etc.), both in physical and numerical terms, as well as their optimum values which best approximate measurements, the model was run several times [13].The time interval of Δt = 3 min, and the time integration parameter of θ = 0,75, were adopted.Coefficients of friction were obtained by calibration with measurements conducted in 2004, and the values of c f = 0,003 (coefficient of friction at the wall of the bed), and c fs = 0,0015 (coefficient of friction on the surface), were adopted [21].It should be noted that actual coefficients of friction on the surface are spatially variable along the boundary surface in the longitudinal direction, and are also time variable due to unsteadiness of speed and depth [22].
Constant values of the coefficient of friction, representing average values of real coefficients, were adopted.Available measurements and measurements used for calibration of model parameters are shown in Table 1.Depths from measured profiles of salinity S = 20‰ (g/l) were adopted for interface position, i.e. water level of Layer (B hB ) Fresh water surface layer density of ρ A = 1000 kg/m 3 and the salt water layer (Layer B) density of ρ B = 1025 kg/m 3 were adopted.
Water levels h A and h B were defined for the initial stage (10th June 2004, 10:00 a.m.) in all nodes, and they correspond to the equilibrium steady state of interface layers.Thus the flow is the same in all elements, and amounts to Q A = 247,52 m 3 /s, at the upper layer, and to Q B = 0,0 m 3 /s (salt wedge stopped).Salinity measurements conducted in Opuzen at that initial stage show that the actual situation is different (h B Opuzen = -4,77 -meters above sea level).The interface is slightly higher when compared to the assumed initial state, and flow gradients exist in both layers in longitudinal direction.Such an unsteady solution cannot be assumed, because it is the result of joint activity of all climatic and hydraulic factors over the past 2-3 days.Hence the initial steady state is assumed.Such a condition is desirable in numerical terms because water levels have a slight gradient in longitudinal direction, and discharge rates are constant.In this way, during simulation of this problem, the error in estimation of initial condition disappears in a very short time (2-3 days), i.e. after 2-3 days the initial condition no longer influences solution functions -water levels and discharges which, at the end of this period, become balanced in accordance with previously defined boundary conditions.

Presentation and analysis of results
The calculated water table (h A ) and the interface (h B ), and the corresponding measurements for the three dates (7 July, 20 August, and 5 October), are shown in Figure 8.The water level calculated for Layer A (h A ) is only slightly different from water levels recorded in Opuzen and Metković, and amounts to just a few centimetres.Average water level errors (measurement deviations) for Layer B differ in individual cross sections and amount to 46 cm at the estuary, 23 cm in Rogotin, 34 cm in Opuzen, and 55 cm in Metković.It can be seen from the figures that the model interface averages measurement results.Therefore, the solution obtained can be regarded as a good-quality assessment of the real halocline, and it can be assumed that the friction coefficients adopted are real.
In the course of June, fresh water inflows were highly variable: an average value was 200 m3/s, and peak values were in excess of 350 m 3 /s.For that reason, the wedge tip was mostly situated between Opuzen and Metković (Figure 9.a).The ensuing period was characterized by practically constant fresh water inflow of about 70 m3/s, and the situation remained unchanged until the end of September.Brief and relatively small increases in flow rate, registered in that period, did not have a considerable influence on wedge displacement.The wedge penetrated deeply inland, to the upstream of Metković (Figure 9.b).The fresh water flew as surface layer of variable thickness, which amounted to about 1,0 m at the estuary, about 2,5 m in Opuzen, and about 3,5 m in Metković.
In October, due to increased rainfall in the basin, the inflow of fresh water increased in the lower course of the Neretva, and so in the second half of October the discharge rates reached values of over 450 m 3 /s (Figure 7.).As a consequence, the salt wedge tip was moved up to Rogotin (Figure 9.c).

Consclusion
Numerous experts, scientists, and organizations operating all over the world, deal with various problems relating to coastal waterway salination.Physical parameters are measured on rivers, and appropriate numerical models are developed, so that these processes can be better understood and adequately described.
Intensive measurements of salinity have been conducted over the past few years in Croatia at the lower course of the Neretva River [4].It was established that in this river the flow of water is stratified at any given moment.Measurements made so far represent an appropriate basis for development of numerical models.The stratified flow (salt wedge) is typical for rivers emptying into seas and having a relatively small amplitude of tides (micro-tidal seas).
Processes such as watercourse salination and circulation of freshwater, sea water, and brackish water at river mouths, rank among highly complex dynamic and stochastic processes.By introducing certain assumptions, an one-dimensional model of stratified flow, without mixing between layers ("sharp interface"), has been developed.The finite element method Igor Ljubenkov, Mijo Vranješ is presented.In this area, the intrusion of salt sea water is most intensive in summer period and, on the other side, the need for clean fresh water is highly pronounced because of agricultural production which is highly developed in this region.The model calibration was made based on measurements conducted in 2004 [4].The interface of layers obtained during this study deviates only slightly from corresponding measurements in longitudinal cross section of the riverbed, which is why this solution can be considered as a good-quality estimation of actual haloclines.It can be concluded from results of this study that the salt wedge reaches the town of Metković in case discharges amount to about 180 m 3 /s, the town of Opuzen if discharges are approximately 280 m 3 /s, and Rogotin if discharges amount to about 450 m 3 /s.It is estimated that, in case of discharges in excess of 500 m 3 /s, the sea water would not intrude into the riverbed, as the sea is in such cases fully pushed out of the river.These values approximately correspond to stationary conditions, and should be taken as information only.It should be noted that the salination process is a pronouncedly non-stationary phenomenon.In addition to quantities of freshwater coming into the system, it is also significant to take into account the duration of such flows as, due to inertia of water, the system is not able to react instantaneously.For that reason, the length of riverbed salination can vary considerably for a given flow, depending on whether the sea water is in the phase of entering the river, or in the phase when it is being pushed downstream.
If we take into account the mentioned typical values of flow, and annual fresh water discharges from the upper part of the basin, it can be concluded that the wedge reaches the town of Metković and/or goes more upstream in several intervals during the year, and that the total annual duration of such intervals is about four months.The riverbed is not salinated only in brief intervals when freshwater discharges are high, and the total annual duration of such high discharges amounts to approximately one month.
In case of minimum flow of 70 m 3 /s, fresh water flows in a surface layer of about 1 m in thickness at the estuary, while its thickness varies from 2 to 3 metres in the upstream part from Rogotin to Metković.
For further research and better modelling of salination, it would also be possible to focus on other possibilities such as the mixing between layers, more accurate establishment of the downstream boundary condition, determination of the model for freshwater spreading in the coastal basin, etc.It is first of all necessary to continue with the current monitoring as this will certainly enable us to gain more knowledge about these highly complex physical processes.
Regardless of the above mentioned additional possibilities, the results shown in this paper provide a good quality simulation of salt wedge displacement at the river estuary, and this model can be considered appropriate for practical use in development of various water management solutions.
To protect the Neretva riverbed against salination, it would be possible to realize a hydrulic barrier that would stop the salt wedge penetration.

Figure 4 .
Figure 4. Schematic presentation of stratified flow (Layer A -fresh water, Layer B -salt water)For the previous assumption, the stratified flow is fully described with four dependent variables -discharge in Layer A (Q A ), discharge in Layer B (Q B ), water level in Layer A (h A ), and water level in Layer B (h B ).These variables are functions of space x and time t.Four equations are needed to define four dependent variables Q A , Q B , h A and h B : two continuity equations and two dynamic equations (derived from the principles of conservation of mass and momentum).
The boundary condition function (water level graph or hydrograph) is discretized by polygon in time domain.The arrangement of discrete points depends on the rate of change of the boundary condition function.5.Example -mouth of the Neretva river 5.1.Topological scheme 82 cross sections were used to describe the river Neretva riverbed from the point where it empties into the Adriatic Sea GRAĐEVINAR 64 (2012) 2, 101-113 Numerical model of stratified flow -case study of the Neretva riverbed salination (2004)

Figure 6 .
Figure 6.Longitudinal section with model boundaries

Figure 7 .
Figure 7. Boundary conditions: a) freshwater inflow (Q A ), b) downstream water levels (h A and h B )

Figure 8 .Figure 9 .
Figure 8.Comparison of model solutions for water level h A and interface h B with measurements for: a) 7 th July, b) 20 th August, c) 5 th October 2004.(measurements profiles S1-S6)