Modelling the suspended matter distribution in an estuarine system : Application to the Odiel river in southwest Spain

ELSEVIER Ecological Modelling 87 (1996) 169-179 E(OLO61(gL mODELLIn6 Modelling the suspended matter distribution in an estuarine system Application to the Odiel river in southwest Spain R. Perififiez b,., J.M. Abril b, M. Garcla-Le6n a a Opto. Ffsica Atdmica, Molecular y Nuclear, Universidad de Sevilla, Apdo. 1065, 41080-SeviUa, Spain b Dpto. Ffsica Aplicada, E.U. Ingenierfa Tdcnica Agr{cola, Universidad de Sevilla, Ctra. Utrera kin. 1, 41014-SeviUa, Spain Received 28 March 1994; accepted 25 January 1995 Abstract A numerical model which solves the advective-diffusive dispersion equation for suspended matter and includes the deposition and resuspension is presented. The model requires the simultaneous solution of the hydrodynamic equations under tidal dynamics and atmospheric forcing, using time steps of a few seconds. The model has been applied to the Odiel river (southwest Spain). The hydrodynamic module has been widely validated for neap and medium tides, whereas the dispersion was calibrated against the dissolved 226Ra dispersion pattern. This 226Ra is discharged to the Odiel river from a phosphate fertilizer factory. The model was able to reproduce the observed behaviour of the suspended matter in the estuary. The sedimentation rates have shown that a net, although slow, sedimentation is being produced. Sensitivity tests were inconclusive with respect to parameters describing settling and resuspension, as internal processes within the estuary are overridden by the high influx and efflux of particulate material from the sea. Keywords: Estuary ecosystems; Hydrodynamics; Sediments; Tides 1. Introduction An estuary is a semi-enclosed coastal body of water which has a free connection with the open sea. Estuaries form the transition from river to sea, so they are influenced by conditions in the river as well as in the coastal sea. The suspended matter is supplied from various sources: atmo- * Corresponding author. sphere, river inflow, erosion of the estuary bottom, shore erosion, coastal sea, waste disposal, runoff etc. Depending on the local situation one or more sources will dominate the supply. In large rivers, with a high suspended load, the river supply dominates. In tidally mixed estuaries the river inflow will dominate in the inner part and the supply from the coastal sea will dominate in the outer part of the estuary. The knowledge of suspended matter dynamics is particularly important to study the dispersion of non-conservative substances in aquatic envi0304-3800/96/$15.00 © 1996 Elsevier Science B.V. All rights reserved SSDI 0304-3800(95)00026-7 170 R. Peridfiez et al. / Ecological Modelling 87 (1996) 169-179 ronments (Abril and Garcfa-Ledn, 1993a,b), since its distribution and interaction with sediments will influence the behaviour of such substances. Recently two of the authors developed a model to describe the basic aspects of suspended matter dynamics over large time scales (Abril and Garcfa-Ledn, 1994). The model was formulated in terms of residual water circulation and mean settling and resuspension velocities. This produces a mean (annually averaged) suspended matter and sedimentation rate distributions. Here we are going deeper in the short scale (both spatial and temporal) aspects. Thus, we present a mathematical model which solves the advective-diffusive dispersion equation with the resuspension and deposition terms, both depending on the instantaneous water state. The model requires the simultaneous solution of the hydrodynamic equations under tidal dynamics and atmospheric forcing, using time steps of a few seconds. The formulation of the deposition and resuspension processes is more complex than the one used in models which work with residual circulations. Our description allows studying the influence of tidal oscillations in suspended matter concentrations and in sedimentation rates. This information is essential in modelling the dispersion of non conservative substances in a non equilibrium situation: the transfer of such substances among the dissolved and solid phases depends on the instantaneous suspended matter concentration and instantaneous sedimentation rate. The model has been applied to the Odiel river (southwest Spain), which is an estuarine system affected by tidal dynamics (M 2 is the main component). The Odiel river is surrounded by a marsh area and an industrial complex, in which two fertilizer processing plants are located (see Fig. 1). These factories release part of their wastes directly to the Odiel river and so important concentrations of U-, Th-(Martfnez-Aguirre et al., 1994) and Ra-isotopes (Perififiez and GarclaLe6n, 1993; Perififiez et al., 1994a) have been measured in its waters and suspended matter. A 2D, with high temporal resolution, model of the suspended matter dynamics is necessary to study the 226Ra dispersion in the Odiel river, ~2 e /V ~ HUELVA ! FER~I~ FZER COMPLEX 11kin I SEA Fig. 1. Map of the Odiel river (southwest of Spain). The part of the river which is covered by the grid is also shown. since ionic exchanges must be taken into account. The model must be 2D because the 226Ra source is very local and, on the other hand, there is a point source of suspended matter in the Odiel river from a mining factory. This produces a bidimensional structure in suspended matter concentrations and sedimentation rates. The high temporal resolution is needed to study the influence of tidal oscillations in the results. The mathematical model is presented in see section 2 and its application to the Odiel river is shown in see section 3. Finally, results are presented and discussed in see section 4, where several sensitivity tests are also shown. These tests are used to study the model response to changes in the parameters involved in the equations. 2. The model The estuarine system in which the suspended matter dynamic is to be studied is represented by a grid containing a certain number of compartments or grid-cells. Each compartment has coordinates (x,y) which define its position in the grid and contains a certain concentration m(x,y) of suspended matter in ppm (parts per million). We assume a homogeneous distribution of m in the water column, that is, vertical discretization is not IIIIII IIIIII IlIIIl IIIIII IIlIIlI IIIIIII IIIIIII IIIIIII IIiliII IIiIIII IiliIII IIiIiII IIIiIII IIIIIII IIIIIII IIIIiII IIIIIII mmmmmmm mmmmmmm mmmmmmm mmmmmmm mmmmmmm mmmmmmm mmmmmmmm mmmnmmmm mmmmmnmmm mmmmmmmmm mmmmummmm mmmmmmmmm mmmmmmmm~ nmmmmmm~m mmmmmmmmmm mmmmmmmmmm mmmmmmmmm IIImIIIIII mimmmiiiii iiiiiiiii IIIIIIIII R. Peridfiez et al. / Ecological Modelling 87 (1996) 169-179 b ¢ 171 ,go 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1000 0 1000 0 1000 Fig. 2. (a) Bathymetry of the Odiel river; depths are shown in m. (b) The grid used in the model. (c) Distribution of small particles in the top sediment layer of the Odiel river, pmf(X,y)/pw is given in %. considered. This assumption is realistic since the estuary is very shallow (see Fig. 2a) and the stream flows are very low, ranging from 4 to 50 m 3 s1 in usual conditions (Perififiez et al., 1994b). Only in the case of very heavy rainfalls stream flows could reach higher values. Nevertheless, the mixing of salt and fresh water takes place far upstream of the studied area. A fast dispersion of the fresh water into a much larger volume of salt water occurs (Borrego and Pend6n, 1988), as usual in tidal mixed estuaries of small rivers (Eisma, 1993). By convention, particulate matter in suspension is defined as the material that is retained on a 0.4-to 0.5-1xm pore size filter. Smaller material is considered to be dissolved (although it may be colloidal or particulate). The upper size limit of particulate matter in suspension is not fixed; heavy particles sink rapidly to the bottom but very large, low-density structures can remain in suspension for a long time. Nevertheless, in order to simulate the suspended matter dynamics, we have to accept a simplification which is usual in this kind of studies. Two size fractions of particulate matter are defined: one with a diameter < 62.5/.~m and one with a diameter > 62.5 txm. We will consider that the first class can remain as suspended matter; the second class will sink rapidly so it will not be included in this study since its horizontal movement is negligible (Belderson, 1964; Guburt et al., 1987). On the other hand, particles with 4> > 62.5 /xm will not be resuspended from sediments. The water velocity at which particles begin to be resuspended is called the critical resuspension velocity. Its value depends on the particle diameter and the roughness length factor (see section 2.1). Typical values for the roughness length range from 0.02 cm for muds to 0.3 cm for gravels (Pugh, 1987). The critical resuspension velocity increases with the particle diameter (larger particles need more energy to be resuspended) and decreases with the roughness length factor: if the roughness is small the sediment is more cohesive and then it is more difficult to resuspend. Roughness lengths in the Odiel estuary range from some 0.1 to 0.3 cm, thus, water velocities of about 0.25 m s -~ are needed to resuspend sediments with ~b ~ 100 ~m (Pugh, 1987). Water velocities in the Odiel river are larger than this value only during a short time, when the water level is increasing or decreasing in the southern part of the river, but as we move 172 R. Peridhez et al. / Ecological Modelling 87 (1996) 169-179 upstream water velocities decrease. Thus, it seems that to be a good approximation to take the limit diameter as 62.5/xm, since larger particles will be resuspended only in the south of the river during very short times. Several processes contribute to changes in m(x,y): the deposition of particles on the estuary bottom and the resuspension of particles from it because of the drag provided by the turbulent motion of the water. They provoke a vertical movement of m(x,y) inside the compartment. The exchange of matter with the adjoining compartments contributes to a horizontal movement of the suspended matter. As a boundary condition, suspended matter input from waste disposal and/or runoff can take place in the compartment. 2.1. Vertical transport The suspended matter in a water column of height h falls down with a mean settling velocity v s. So the sedimentation process contributes to a variation in time of m which, in the case of a homogeneous distribution of m in depth, is: 0m v~ ( Iv[ t 0~- - ~m 1 - ~-~¢~ ] (1) if [v[ < yes, where v is the depth-averaged horizontal water velocity and yes is a critical settling velocity. Thus, deposition of suspended matter occurs only when the water velocity is below a critical value, as observed in nature. If Iv[ > Vcs, there is no sedimentation due to the turbulence: particles remain in suspension and, on the other hand, flocs are broken. A value of 0.18 m s -l has been measured for Vcs (Eisma, 1993). Below this velocity, deposition of all particles in suspension Occurs. The settling velocity depends on the concentration of suspended matter: clouds of particles settle faster than single particles because a large number of particles in the cloud settle in the wake of others. It has been found (Mentha, 1989) that the settling velocity can be written as: v s = 1.74 × 10 -9 m 1"6 (2) if v s is measured in m s-1 and m in ppm. This equation can be used when the concentration of suspended matter is not higher than ~ 1000 ppm because at such concentrations the falling particle is hindered by the other particles in suspension. Indeed, at concentrations of 105 ppm the settling is negligible (Eisma, 1993). Anyway, the actual concentrations remain below ~ 100 ppm. The resuspension effect, produced by the shear stress, will contribute to the variation of m over time with an amount of (in ppm) (Abril and Garcia-Ledn, 1994): ~m vr Pm f -- 10 6 (3) ~t h Pw where u r is the resuspension velocity (m s-1), Pm is the dry matter density of the sediment, f is the weight fraction of small particulates (~b < 62.5 /zm) in the top layer of the sediment in compartment (x,y) and Pw is the water density. Laboratory measurements have shown that as the current speed is gradually increased from zero, there is a speed at which the sediment begins to move. This is called the critical resuspension velocity vet. Its value depends on the roughness length factor, z 0, which depends on f. For instance, z 0 = 0.02 cm for muds and z 0 = 0.3 cm for gravels (Pugh, 1987). From data shown in Pugh (1987) we have constructed an analytical function which relates Vcr with z 0 for a particle diameter of about 100 /zm (since there is no information for a diameter closer to 62.5 tzm). Thus, Ucr : 0.441 e - 1'117z°°231 (4) if Vcr is measured in m s-1 and z 0 in cm. Taking into account the critical resuspension velocity, the resuspension term has been rewritten as: Om VrPmf lo6( IUI _1 ) fi~w ~Ucr (5) This formulation is valid if Iv[ > Ucr. If the current speed is not higher than the critical resuspension velocity, the resuspension term will be zero. R. Perid~ez et al./ Ecological Modelling 87 (1996) 169-179 173 As shown in (see section 4), there are no important changes in the model results if v r is slightly varied, so we have used the value found in the current literature: u r = 0.03 m year1 (Abril and Garcia-Le6n, 1994). 2.2. Horizontal transport The suspended matter can remain in the water column for a long time, so it will participate in the water movements and will behave as a conservative substance. The exchange of suspended matter with the neighbouring compartments is expressed by the advective-diffusive dispersion equation (Peri~ifiez et al., 1994b): 3m ~m 3m l[~x( ~-X) + U + V- = -- hK u Ot Ox ~y h +3-7 where u and v are the components of the water velocity in the direction of the x-and y-axis respectively and K, and K,, are the diffusion coefficients in the corresponding directions of the x-and y-axis, which can be written as (Perififiez et al., 1994b): K. =/3,1ulf~u z + v z (7) K,~ =/321vl~u2 + v 2 where /31 and /32 (dimension [T]) have to be calibrated for each specific site. To solve the advective-diffusive dispersion equation, the water velocities and the height of the water column must be known for each compartment and for each time step. They are obtained by solving simultaneously the hydrodynamic equations (conservation of mass and momentum), which include the Coriolis term, bed friction, response to wind stress and response to changes in atmospheric pressure. Spatial gradients in atmospheric pressure have not been considered because of the small dimensions of our estuarine site (about 4 km length). A complete description of the hydrodynamic equations can be seen in Perififiez et al. (1994b). The equation which describes the dynamic of suspended matter m in compartment (x,y) is obtained by summation of all the above-mentioned terms: advective, diffusive, deposition, resuspension and a source term, which takes into account external sources of suspended matter. The solution of that equation provides information on the sedimentation processes that take place in the area under study. The sedimentation rate w is obtained as the net balance between sedimentation and resuspension (or erosion). Expressing w in kg m -2 s -~, we have: ( Ivl t [Ivl ) w=GmPw 1--- 10 -6Vcs ] VrPmf[~cr 1 (8) 3. Application of the model to the Odiel river 3.1. Geological structure of the Odiel river The Odiel river mouth forms a sedimentary environment. It is an estuarine system affected by tidal dynamics. The outer part of the river was invaded by the sea approximately 5000 years ago, when a big bay was formed over neoquaternary detritus (Borrego and Pend6n, 1988). Sediments are introduced in the estuary from several sources: River inflow." The Odiel river introduces, in the inner part of the estuary, thick-grained sands. This material contains important amounts of iron oxides, which originate from erosion of Palaeozoic rocks. Muds are also introduced during torrential rain episodes. Sea supply: In the outer part of the estuary. Siliceous medium and thick-grained sands are observed, as well as carbonated material (principally shells). Muds: Sea origin, from the flocculation of dissolved particles. This is the greatest amount of material introduced into the estuary. A map of the part of the Odiel river under study is shown in Fig. 2, where the grid used in our model is also represented. It includes 456 compartments, each of which can be water (white) or land (black). Compartments are described by 174 R. Perid~ez et al. / Ecological Modelling 87 (1996) 169-179 the spatial coordinates (x,y) and depths were introduced as input data from marine charts. The dimensions of the compartments are Ax = Ay = 100 m. The available information for f(x,y) (Universidad de Sevilla, 1991) is given in Fig. 2c, where pmf(X,y)/pw is represented in %. 3.2. Hydrodynamics of the Odiel river The hydrodynamic equations have been calibrated for the Odiel river for neap and medium tides, as well as the advective-diffusive dispersion equation, which has been used to study the 226Ra dispersion in the Odiel river treating it, as a first approach, as a conservative substance (Perififiez et al., 1994b,c,d). In these references the values obtained by calibration can be found for parameters which appear in the hydrodynamic and dispersion equations. These include the bed friction coefficient and /3 t and /32 in the diffusion terms. Boundary conditions are also described in detail. Elevations (m) -1 500 1000 1500 2000 2500 Velocities (m/s) 0.5 05 I i 500 1000 1500 2000 2500 time (min) Fig. 3. Time evolution of the water elevations and velocities along the y-axis in a compartment near the southern border (--) and a compartment in the middle of the grid (---). The differences in elevations are a few mm, so elevations for both compartments are represented by the same curve. mBq/I 600 Lines: model results • Points: experimental data 00 20O 0 i ~- 10 20 30 Position in the grid -- Low water • Low water -- High water A High water Fig. 4. Model results for 226Ra dispersion. Points are experimental data (Peri~fiez and Garcla-Le6n, 1993) and lines are the computed 226Ra concentrations for high (---) and low (--) water. Values in mBq/l. The x-axis corresponds to the localization in the grid. For details on the method used to solve these equations we also refer to the references. A centred finite differences scheme was applied. The time step was At = 6 s, so the CourantFriederich-Lewy criterion (Prandle, 1984) is met and the numerical dispersion is negligible (Perififiez et al., 1994b). In Fig. 3 the course of the water elevations and velocities along several tidal cycles for medium tides are shown. Water velocities are similar to those experimentally obtained. The calculated maximum water velocities are 0.45 and 0.61 m s -1 when the water level is increasing and decreasing respectively, while the measured ones were 0.48 and 0.66 m s-l. In the case of neap tides, the maximum calculated velocities are 0.38 and 0.29 m s -1, while the measured maximum velocities were 0.40 and 0.28 m s -1 when the water level is increasing and decreasing respectively. In Fig. 4 results of the application of the dispersion module to 226Ra dispersion are compared to field data. As can be seen in Fig. 4 the agreement is rather good; the real peak is more intense than the calculated one because that sample was collected just in the effluent of one of the fertilizer plants, and we are assuming instantaneous homogenization inside each compartment. R. Peridfiez et al. / Ecological Modelling 87 (1996) 169-179 175 Elevations (cm) ppm 350 80 300 250 200 150 100 50 0 60 40 20 9 10 11 12 1 2 3 4 5 6 7 8 Time (hours) Elevations (cm) ppm 350 80 300 250 60 200 40 150 1 O0 20 5O 0 0 8 9 10 11 12 1 2 3 4 .5 6 7 Time (hours) ~on -F S. matter~ Fig. 5. Suspended matter concentrations (ppm) in the northern (a) and southern (b) borders of the grid over a tidal cycle. Elevations are given with respect to an arbitrary reference. 3.3. Boundary conditions and sources of m(x,y) following a method similar to that of Durrieu de Madron et al. (1992). As said before, suspended matter concentrations in the southern border are used as boundary conditions, while those of the northern border will be used to calibrate the model. On the left shore of the Odiel river there is an industrial complex. One of the factories located here discharges important amounts of material into the river, so the source term S must be considered. The compartment in which the source term is located is (7,9). In order to evaluate the magnitude of this term, we must take into account that the =~Ra and U-isotopes activity concentrations measured in sediments collected at the point of discharge are about a factor 10 smaller than the activity concentrations measured in sediments of the surroundings (MartinezAguirre et al., 1994). Thus, the sedimentation rate when the discharges of material are performed must be about 10 times larger than when there are no discharges. This is due to the fact that the material discharged is basically Fe, and consequently there is a reduction in the U and Ra content (Respaldiza et al., 1993). These considerations have allowed us to find the source term in a calibration exercise, which is 0.78 kg/s. The boundary conditions for m(x,y) in the northern border of the grid are the ones which were used when the advective-diffusive equation was calibrated. They consists of taking the suspended matter concentration in the last row equal to that of the previous row: m(x,38) = m(x,37) Along the southern border, suspended matter concentrations were specified for each time step. These concentrations were experimentally measured in the Odiel river estuary at about 0.5 m below the water surface. To do this, a sampling campaign was performed: water samples were collected every 20 min during a complete tidal cycle in both the southern and northern borders of our grid. The results are shown in Fig. 5. Suspended matter concentrations were measured using a nephelometer, previously calibrated to convert ntu (nephelometric turbidity units) in ppm 4. Results and discussion 4.1. Suspended matter distribution and sedimentation rates To study the suspended matter distribution, a simulation over several tidal cycles was performed in a medium tide situation. The course of the suspended matter concentration in three compartments of the grid can be seen in Fig. 6a. As we move towards the northern border, the oscillations in suspended matter concentrations, which are due to the tidal oscillations, decrease. Indeed, the suspended matter concentration in the compartment near the northern border is almost constant. This behaviour has been observed in the Odiel river: Fig. 5 shows that concentrations at the northern border are quite constant, the mean value is 28 + 4 176 R. Peridfiez et al. / Ecological Modelling 87 (1996) 169-179 ppm. The computed mean value is 24.4 ppm, so the agreement is rather good. This behaviour of the suspended matter concentrations along the 0 o 0 ....................................................................... i . . . L i . i , i , 2500 3000 1500 2000 TIME (rnin) a O O o @4 O 8 b o O o ro cq O v-- o o 0 I000 0 I000 Fig. 7. Suspended matter concentration (ppm) maps when water level is increasing (a) and decreasing (b). ' i i i b tO q o 1500 2000 2500 COO TIME (rain) Fig. 6. (a) Time evolution of the suspended matter concentrations in three compartments of the grid: near the southern border (--), in the middle of the grid (---) and near the northern border (- • -). (b) Time evolution of the sedimentation rate (g cm -2 year -l) in a point located in the middle of the grid (--) and in a point near the northern border (---). The moments at which high and low water occur are also shown. Odiel river is due to the fact that the water velocities decrease as we move towards the northern boundary of the grid (see Perififiez et al., 1994b and Fig. 3), so the resuspension becomes negligible and the deposition must be almost constant. Thus, oscillations in suspended matter concentrations vanish. On the other hand, it can be seen in Fig. 6a that there are two peaks in the compartment near the southern border for each tidal cycle: when approaching high and low water levels, the velocities are higher and there is an important resuspension of matter from the river bed and concentrations reach the maximum values. But during high and low water only deposition occurs and the concentrations decrease, reaching minimum values. Suspended matter concentration maps when water level is increasing and decreasing can be seen in Fig. 7a and b respectively. The decrease in concentrations as we move towards the north is again clear, as well as the effect of the source of suspended matter, which produces a local concentration increase. The development of the sedimentation rate in R. Perid~ez et al,/ Ecological Modelling 87 (1996) 169-179 177 o 0 C,,I All these results correspond to a medium tide situation. Extreme conditions, such as torrential rains and equinoctial spring tides, can alter the results. These extreme conditions are presently under investigation. In the case of neap tides, the time evolution of suspended matter concentrations and sedimenta0 ,r-- C:I O (D 150 100 0. 50 ii. 0 lO 20 30 C~ 0 1000 Fig. 8. Sedimentation rates (g cm -2 year -1) averaged over several tidal cycles. two compartments of the grid is shown in Fig. 6b. The sedimentation rate ranges from some 0.02 to 0.07 g cm -2 year-~. It remains positive, so there is a net sedimentation in the area although the process is slow, as confirmed by the low values of the sedimentation rates. As can be seen in Fig. 6b there are two peaks for each tidal cycle in the compartment located in the middle of the grid: one during low water and one during high water. In these situations the water velocities are very low, being the deposition the dominant process. In the compartment near the northern border the pattern is the same but the peaks are lower. This is again related to the fact that, due to the low water velocities, the deposition process is almost constant. The sedimentation rates have been averaged over several tidal cycles all along the Odiel river. The results are shown in Fig. 8, where a map of the Odiel river with the sedimentation rates expressed in g cm -z year-~ is shown. The highest value of the sedimentation rate, 0.3 g cm -2 year-~, corresponds to the point of discharge. In the rest of the river the sedimentation rates are smaller, showing a net deposition of material all over the river bed. 150 100 50 0 __ , ~ ' - , , , b 0 10 20 30 v) i r i C .fi J , f~ i, q ,' \ , ~ ! , '~ ~ "~ ! , ~. ! ~,~, ., !- ; • i ~"~;J t C:~ , ~ L • , J , i 1 1500 2000 2500 3000 TIME (mln) Fig. 9. Suspended matter concentration for low and high water (a and b) and time evolution of the sedimentation rate (c) when the source term is doubled (- • -) and when there is not a source term (---). The continuous lines (--) are the model results.