Assessment of numerical schemes for solving the advection – diffusion equation on unstructured grids : case study of the Guaíba River , Brazil

In this work, a first-order upwind and a high-order flux-limiter schemes for solving the advection–diffusion equation on unstructured grids were evaluated. The numerical schemes were implemented as a module of an unstructured two-dimensional depth-averaged circulation model for shallow lakes (IPH-UnTRIM2D), and they were applied to the Guaíba River in Brazil. Their performances were evaluated by comparing mass conservation balance errors for two scenarios of a passive tracer released into the Guaíba River. The circulation model showed good agreement with observed data collected at four water level stations along the Guaíba River, where correlation coefficients achieved values up to 0.93. In addition, volume conservation errors were lower than 1 % of the total volume of the Guaíba River. For all scenarios, the higher order flux-limiter scheme has been shown to be less diffusive than a first-order upwind scheme. Accumulated conservation mass balance errors calculated for the flux limiter reached 8 %, whereas for a first-order upwind scheme, they were close to 18 % over a 15-day period. Although both schemes have presented mass conservation errors, these errors are assumed negligible compared with kinetic processes such as erosion, sedimentation or decay rates.


Introduction
Ecology of lakes and ponds is closely linked to physical factors, especially to hydrodynamic variables such as velocity, turbulence, and diffusion/convection of suspended material (Brönmark and Hansson, 2005).These aquatic ecosystems present complex circulation patterns which vary over time and space depending on density, temperature, wind among others (Reynolds, 1984).Fragoso Jr. et al. (2008), for example, revealed a close relationship between spatial heterogeneity of phytoplankton and water circulation patterns -driven by gradients in viscosity and diffusivity, wind stress and bottom friction -in a shallow lake.In addition, Hodges and Dallimore (2013) showed that upward and downward fluxes of water in lakes and/or ponds are associated with vertical heat exchange throughout the water column.Thus, hydrodynamic processes affect transport of substances to, from and within lakes and ponds.The spatial variance of suspended particle concentrations may be great in shallow ecosystems once their distribution can be influenced by short-term physical factors (e.g.sediment resuspension) (Carrick et al., 1993).Thus, a better understanding of the water circulation of these ecosystems plays an important role in their water quality dynamics (Chapra, 2005).
Several studies have presented numerical models which are able to represent velocity field and water level by solving the shallow water equations.Some of these models use finite difference scheme in a horizontal plane formed by uniform rectangular grids (Casulli, 1990;Cheng et al., 1993).However, occasionally uniform grids are not flexible enough to represent complex geometries.Therefore, some numerical models use finite difference schemes on nonorthogonal curvilinear grids to allow greater flexibility (Ye and McCorquodale, 1997;Zhou, 1995).On the other hand, an disadvantage of using curvilinear grids lies in approximating oblique boundaries along the simulation domain which Published by Copernicus Publications on behalf of the European Geosciences Union & the American Geophysical Union.introduces errors to the numerical solution (Margolin and Shashkov, 1999).Moreover, additional terms derived from curvilinear transformation are also sources of errors arising from a numerical approach using curvilinear grids (Hirsch, 1988).
A detailed study on the effects of grid non-orthogonality carried by Sankaranarayanan and Spaulding (2003) reveals that truncation error terms due to first and second derivative terms are functions of grid angle and aspect ratio in a way that root-mean-square (rms) errors in surface elevation and velocities sharply increase as the grid resolution decreases for curvilinear grids with a grid angle less than 45 • (Thompson et al., 1985;Nielsen and Skovgaard, 2005).
In this sense, hydrodynamics models using orthogonal unstructured grids are considered an efficient tool to describe the dynamics of rivers, lakes and estuaries as they have high capacity of representing geometries and optimizing computational efforts with grid refinement in regions of interest (Casulli and Walters, 2000;Cheng and Casulli, 2001).
Once conservation of volume is assured by an efficient hydrodynamic solution, a numerical solution for the advectiondiffusion equation is needed to characterize the transport of a scalar variable such as salinity, temperature or any passive constituent that may represent sediment or biological communities.Analytical methodologies to solve the advectiondiffusion equation by mathematical substitutions and transformations (Mikhailov and Ozisik, 1984;Cotta, 2005) are widely applied to many experimental applications in heat and mass diffusion (Leij and van Genuchten, 2000;Guerrero et al., 2009).Nevertheless, a large number of parameters, coefficients, constants and functions are used in such methodologies, which may reduce their range of applications to systems under strictly controlled external conditions.
Several numerical models use a first-order upwind scheme to solve explicitly the advection-diffusion equation (Gross et al., 2002(Gross et al., , 1999;;Hetch et al., 1995).However, first-order upwind schemes are only effective in regions with low advection and low scalar concentration gradients (Le Veque, 1996).Greater accuracy may be achieved in first-order upwind schemes employing refined grids, which implies a higher computational effort.An alternative approach is to use high resolution schemes with flux-limiter functions.The highresolution scheme developed by Sweby (1984) has been successfully applied and improved over the past years (Casulli and Zanolli, 2005;Wang et al., 2008;Wood et al., 2008).
Therefore, this work presents the formulation and application of a two-dimensional depth-averaged circulation and scalar transport model for shallow aquatic ecosystems based on unstructured orthogonal grids.It has been developed to maximize flexibility in grid specification, ensuring a stable and robust numerical solution.A high-resolution scalar transport scheme using a flux-limiter function developed by Sweby (1984) was implemented and its results were compared with a first-order upwind scheme.The model was applied to the Guaíba River to test its scalar transport schemes in solving realistic problems.Their numerical solutions were compared with estimates and measurements of sedimentation and erosion rates, whereas their performances were evaluated using conservation of mass for the entire simulation domain.

Model description
IPH-UnTRIM2D is a finite-volume numerical representation of the two-dimensional depth-integrated continuity and momentum equations of water motion developed at the Instituto de Pesquisas Hidráulicas (IPH), and conceptually based on the unstructured version of the Tidal, Residual, and Intertidal Mudflat (TRIM) model proposed by Casulli and Walters (2000).It also represents integrated scalar transport (e.g.salinity, heat, suspended sediment) through the advection-diffusion equation.The governing equations and finite-volume numerical approximations are described here.

Hydrodynamic module
IPH-UnTRIM2D uses a finite-volume approach to solve the shallow water equations on unstructured grids.While it calculates values of free surface at the centres of each computational cell of the grid, values of velocity are estimated at the half length of their sides.As the shallow water equations are derived from depth-integrating the Navier-Stokes equations -under the assumption that horizontal flow components are much greater than vertical ones -the hydrodynamic module of IPH-UnTRIM2D expresses the principle of conservation of mass and momentum.
To integrate the Navier-Stokes equation over the water column, IPH-UnTRIM2D assumes two boundary conditions corresponding to effects of wind stresses and bottom friction.From the momentum equations for an incompressible fluid (Fox et al., 2009), the boundary conditions due to effects of wind stresses at free surface and bottom friction at bed are given by prescribing the term that represents the divergence of stress along the water column as: where A v denotes the vertical coefficient of kinematic eddy viscosity, ∂u ∂z and ∂v ∂z represent the gradient of the horizontal components of velocity along the water column, τ x and τ y are the horizontal components of wind shear stresses and γ indicates the bottom friction.
Another assumption used in the hydrodynamic module of IPH-UnTRIM2D is related to the term of pressure gradient force in the momentum equations.Theoretically, the pressure gradient force can be divided into barotropic and baroclinic components by using the Leibniz integration rule (Blaise et al., 2010).And, since the baroclinic components are driven by vertical gradients of density and/or temperature over the water column, IPH-UnTRIM2D includes only barotropic components of pressure in the calculations as part of its depth-integrated approach (i.e.well mixed water column).This assumption reduces the term of pressure gradient force in the momentum equations to the following form: ∇p (x, y, z, t) = g∇η (x, y, z, t) , (2) where p(x, y, z, t) is the pressure gradient force, g is the gravitational acceleration constant and η(x, y, z, t) means the free surface from an undisturbed reference level.
The hydrodynamic module of IPH-UnTRIM2D also incorporates the continuity equation.Similarly to the momentum equations, the continuity equation is integrated over the water column -from bed to free surface -and, along with the shallow water equations, is used by IPH-UnTRIM2D to describe free surface flows.Thus, the momentum and continuity equations for two-dimensional circulation solved by IPH-UnTRIM2D can be written as: where u(x, y, t) and v(x, y, t) are the horizontal velocity components in x and y directions in m s −1 , t is time in seconds, h(x, y) is water depth relative to an undisturbed reference level in m, f is the Coriolis parameter in s −1 assumed as constant and A h is the horizontal eddy viscosity coefficient in m 2 s −1 .
Wind shear stresses (τ x , τ y ) and bottom-friction (γ ) formulations are also calculated in this module.The wind shear stresses at free surface are proportional to the wind speed and they were calculated as follows: where C D is the wind stress coefficient, W x and W y the horizontal components of wind speed at the surface level in m s −1 , and W is the norm of the wind speed vector in m s −1 .The bottom-friction coefficient is calculated by: Here H (x, y, t) = h(h, y) + η(x, y, t) is total water depth and C v is the Chezy coefficient.Coriolis acceleration in the x and y directions used in the momentum equations are given by f v and f u, respectively, with f = 2 ω sin φ where ω means angular speed of rotation of the Earth about its axis in rad s −1 and φ is geographic latitude in degrees.As part of the assumptions for depth-averaged circulation models (e.g.AD-CIRC, Westerink et al., 1994;TRIM, Casulli, 1990;River2D, Heniche et al., 2000), both the effects of wind stresses and bottom friction are included at the same layer in the momentum equations.

Scalar transport module
The scalar transport module in IPH-UnTRIM2D calculates the concentration field of a tracer by means of the advectiondiffusion equation deduced from the law of conservation of mass for incompressible fluid, and is expressed as: where C is the mean concentration in the water column in mg L −1 , H is total depth in m and K h is horizontal eddy viscosity in m 2 s −1 .Considering a conservative tracer, sources/sinks and reaction terms are assumed to be negligible.

Unstructured grid
IPH-UnTRIM2D calculates free surface, velocity and tracer distribution at the centres and sides of computational cells on an unstructured orthogonal grid.An unstructured grid is a set of non-overlapping computational cells over a simulation domain.Each of these computational cells is composed of N p sides, N v vertices and only one centre, which is located neither at its centroid nor geometric centre but at the intersection point of the perpendicular bisectors of N p sides of the computational cell.Since the segment that joins the centres of two adjacent computational cells are orthogonal to each other, this set of non-overlapping computational cells is known as unstructured orthogonal grid (Fig. 1).MTOOL (Two-dimensional Mesh Tool) has been used as a grid generation package in order to generate the unstructured orthogonal grids used in this study.It is an interactive graphic program for two-dimensional finite element mesh generation developed by TeCGraf (1992).It allows one to combine coarser and finer triangular and/or quadrilateral computational cells on the same grid to represent transitions from broad and open regions and to regions of interest.
Figure 1 illustrates an unstructured orthogonal grid.For an unstructured orthogonal grid of N p computational cells, which i = , 2, 3 . . .N p , each ith computational cell has area denoted by P i and 3 or 4 sides represented by j = 1, 2, 3-4.And, for each side j , λ j indicates its length, whereas δ j 1116 F. F. Pereira et al.: Assessment of numerical schemes means the distance between the centres of two adjacent computational cells that share the same side j .

Numerical approximation
A semi-implicit finite-volume scheme was used to obtain an efficient numerical algorithm where its stability is independent of free-surface gravity waves, wind stress, vertical viscosity and bottom friction.In contrast to circulation models based on structured grids which numerically solve the momentum equations for each computational polygon on the Cartesian axes, IPH-UnTRIM2D differentiates the momentum equations for each side j of a computational polygon i with respect to its normal and tangential axes.Normal and tangential components of bottom friction and wind stresses are calculated by rotating their horizontal and vertical components.The bottom friction, wind stresses, Coriolis acceleration and horizontal viscosity are treated fully explicitly.Thus, the conservation of normal and tangential momentum for a side j takes the following forms: where u and v are the normal and tangential components of velocity at time step n + 1 and F is an explicit operator which combines all explicit terms with the backtracked velocity at time step n.The implicitness factor for temporal discretization (θ) may vary from 0 to 1, given that values of θ equal to 0 indicate a fully explicit scheme, whereas values of θ equal to 1 mean that an implicit scheme is being used.For semiimplicit schemes, the implicitness factor is in the range of 1/2 ≤ θ < 1 (Casulli, 1990;Casulli and Cheng, 1992;Zhang and Baptista, 2007).Stability numerical analysis for semiimplicit schemes carried out by Casulli and Cattani (1994) indicated highest accuracy and efficiency for an implicitness factor equal to 0.5, the same value of θ used in this study.Regarding the grid connectivity, i(j, 1) and i(j, 2) represent neighbouring polygons of the j th side, whereas v(j, 1) and v(j, 2) indicate its vertices.Water-surface elevations are computed for all polygons and vertices on the grid.At a generic vertex, the water-surface elevation (η v ) is estimated by area-weighted average using elevations at its surrounding polygons.On the other hand, a set of equations comes from the continuity equation in order to calculate the elevations at all polygons.Thus, the continuity equation can numerically be written for each polygon as: where N ij represents the number of sides of the polygon i and S (i,l) is a function related to flow direction through the lth side of the polygon i, and its value may be equal to 1 or −1.S (i,l) takes the following form: where every j th side of the grid presents neighbouring polygons of indexes [j (i, l), 1] and [j (i, l), 2].Substituting Eq. ( 10) into Eq.( 12), a linear system of N p equations with a symmetric and positive definite matrix can be composed and takes the following form: and being that where N j denotes the number of sides in computational cell i and N p is the total number of computational cells in the grid.
The numerical solution of this linear system can be efficiently determined by the conjugate gradient (Press et al., 1992).This numerical approach enables the hydrodynamic module included in IPH-UnTRIM2D to simulate unsteady flow under the effects of external flow disturbances caused by tidal forces and river inflows/outflows, for example.Such external forcings are interpreted by the hydrodynamic module as lateral boundaries, which are assumed strictly vertical.Theoretically, lateral boundary conditions are provided to the hydrodynamic module as described thereafter.
Lateral land boundaries are set at all sides along the landwater interface.For these sides of the grid, normal flow and drag on tangential flow are equal to 0. For open sea boundaries, boundary conditions are assumed at the centre of computational cells as water surface elevation.River boundaries are set at both sides and centres of computational cells by prescribing normal flow (i.e.discharge) and water elevation respectively.
Once the elevation field at time step n + 1 is known, normal and tangential velocities are also updated at time step n + 1 by using Eqs.( 10) and ( 11).Elevations and velocities computed in the hydrodynamic module are used to calculate the numerical solution of the advection-diffusion equation.
Currently, IPH-UnTRIM2D presents two numerical schemes for solving the advection-diffusion equation.In the first one, IPH-UnTRIM2D employs a first-order upwind scheme to numerically discretize the advection-diffusion equation over the simulation domain based on geometric features of the computational cells, their connectivity with each other as well as elevation and water velocity fields from the hydrodynamic module.The second one appears as an alternative to the first-order upwind scheme, and includes a fluxlimiter function that aggregates a higher-order term over regions of high-velocity gradients to reduce numerical diffusion in the solution of the advection-diffusion equation given by the first-order upwind scheme (Casulli and Zanolli, 2005).
The first-order upwind scheme solves temporal and spatial partial derivatives from the advection-diffusion equation using simple backward differences as shown in Eq. ( 17).
where m(i, j ) is the index of the polygon that shares the j th edge with ith polygon, are advective and diffusive fluxes, respectively.Since S function returns either 1 or −1, S + denotes only grid sizes of S function equals to 1, whereas S − indicates the negative ones.
Similarly to the upwind scheme, the higher-order fluxlimiter scheme also uses backward differences to numerically represent the advection-diffusion equation.However, the higher-order flux-limiter scheme incorporates a high-order term, which depends on a function ( ) that returns 0 (firstorder scheme), 1 (second-order scheme) or 2 (first-order but less diffusive scheme).As proposed by Sweby (1984), the higher-order flux-limiter scheme takes the following form: n j is a flux-limiter function which can be written as: where r is the ratio of consecutive gradients and φ is chosen when second-order accuracy for the flux limiter is sought.Thus, r and φ are calculated according to: Although the higher-order flux-limiter scheme uses a few more functions than the upwind scheme, both schemes are fully explicit.It means that in terms of computational efficiency, they are equally efficient.However, finer grids may impose severe numerical constraints on these explicit timestepping schemes (Chapra, 2005).Regarding computational costs, numerical analysis performed by Pereira (2010) on a Pentium ® 4 revealed that IPH-UnTRIM2D requires 5 min to complete a 24-h period simulation run for a grid composed of 293 computational cells of 100 × 100 km.Over a larger simulation domain of 7527 computational cells, IPH-UnTRIM2D needed a little over 40 min to complete the same 24-h simulation run.For both runs, IPH-UnTRIM2D presented the same CPU time when running with the upwind and flux-limiter schemes.
In the next section, IPH-UnTRIM2D is applied to the Guaíba River.Thereafter, its capability in predicting circulation and scalar transport patterns is discussed.).An auxiliary grid (dashed lines) composed of segments that join the centres of each computational cell is also shown to highlight geometric variables used in the section 2.9 Fig. 1.Illustration of an unstructured orthogonal grid (solid lines).An auxiliary grid (dashed lines) composed of segments that join the centres of each computational cell is also shown to highlight geometric variables used in Sect.2.9.

Study area
The Guaíba River is a large (surface area 436 km 2 ) and shallow (mean depth 6.0 m at full pool) water body that runs through Porto Alegre in the south of Brazil.It plays a fundamental role in transport, irrigation, drinking water supply and wastewater discharge for the region.Its length measures 50 km, whereas its width at certain points has a cross section of up to 15 km.Moreover, the Guaíba River comprises one of the most important freshwater systems in Rio Grande do Sul (Fig. 2).
According to the Guaíba directive plan, its major environmental problems are domestic sewage and industrial waste effluents from urban areas.In addition, rural areas present water pollution by pesticides.Therefore, the Guaíba River is often subjected to eutrophication processes due to the nutrient enrichment.It becomes more evident during the summer once high levels of temperature and sunlight provide optimal conditions for phytoplankton blooms.

Input data
The IPH-UnTRIM2D model requires multiple input data sets which are collected from different databases.All time series data should match the time period for which calibration and predictions are being made.Thus, input data sets were obtained for 1991 due to their consistency over the entire year.
Hourly wind speed and direction data were recorded by the Instituto Nacional de Meteorologia (INMET) at a meteorological station located in Porto Alegre.Daily water level and discharge data were provided by Agência Nacional de Águas, Brazil (ANA).The water level stations are situated at five sites along the Guaíba River: Ilha da Pintada, Ipanema, Cristal, Ponta dos Coatis and Ponta Grossa (Table 1).
The discharge stations at the Jacuí, Sinos, Caí and Gravataí rivers were used as upstream boundary conditions, whereas the water level station at Ponta dos Coatis was assumed as the downstream boundary condition.Since there are no discharge stations at the interface of the Jacuí and Guaíba rivers, continuous discharge values for the Jacuí River were estimated based on a regionalization method (Tucci et al., 1995;Samuel et al., 2011).The regionalization method involved transferring discharge data from subsidiaries of the Jacuí River to its outlet by using interpolation techniques that depend on geographic location.The discharge stations used as input to the flow regionalization method are located 20 km upstream from the interface of the Jacuí and Guaíba rivers, as shown in Fig. 2.
Once the critical flow conditions were reached in the main channel, a gradual variation of element size from the  shoreline to the main channel was used, where larger elements were placed along the main channel.MTOOL has therefore created a triangular grid composed of 4622 nodes, 7527 sides and 12 156 triangles (Fig. 3).
A bathymetric survey of the Guaíba River was carried out by the Diretoria de Hidrografia e Navegação do Ministério da Marinha in 1964 (Fig. 3).Therefore, average depths at each element were estimated by a linear interpolation.

Hydrodynamic module calibration
The numerical approximation of the momentum equations provided parameters related to fluid viscosity and diffusivity, wind stress, bottom friction and implicitness for the temporal discretization.Many of these parameters have already been tested and calibrated for rivers and open-channels (i.e.Chezy's roughness and wind drag coefficients) (French, 1986).Therefore, Chezy's coefficient (C z ) was adopted to be constant and equals to 44.7 m 1/2 s −1 over the entire domain.A calibration of wind drag coefficient (C D ) for hydrodynamic models was carried out Escalante Estrada et al. ( 2010).According to their numerical experiments, the wind

F. F. Pereira et al.: Assessment of numerical schemes
drag coefficient was assumed to be 0.016.Regarding the eddy horizontal viscosity (A h ) and diffusivity (K h ) , Fragoso Jr. et al. (2008) has successfully used 5 and 10 m 2 s −1 , respectively, for modeling spatial heterogeneity of phytoplankton in a shallow lake at South Brazil.For practical applications, the implicitness factor (θ) is recommended to be in the range 0.5 ≤ θ ≤ 1 (Cheng et al., 1993;Zhang and Baptista, 2007).However, numerical analysis performed by Casulli and Cattani (1994) has shown that the θ method was stable and presented the highest efficiency and accuracy for a value of θ equals to 0.50.
For calibration, the IPH-UnTRIM2D model used a data network composed of a precipitation station, five water level stations and a meteorological station.A 150 day-period from 1 January to 1 June in 1991 was chosen in order to represent a wide range of hydrologic events (i.e.daily floods and droughts).
The hydrodynamic module was calibrated by reproducing water level observations at the Ilha da Pintada, Ipanema, Cristal and Ponta Grossa stations (Fig. 2).Afterwards, the conservation of volume was also tested where the balance between incoming and outgoing water fluxes must be equal to the volume due to water level variations in the Guaíba River.

Assessment of advection-diffusion numerical schemes
The water quality module was tested based on the conservation of mass for both numerical schemes to solve the advection-diffusion equation.The efficiency of the higherorder flux-limiter scheme was compared with a first-order upwind scheme by their conservation of mass for the Guaíba River.In order to perform these comparisons, two scenarios were considered for development of numerical analysis.The mass balance computations for the higher-order flux-limiter and first-order upwind schemes were compared with each other at every time step.The first scenario reproduces a deliberate release of a tracer at a steady rate of 5 mg L −1 into the Guaíba River over 15 days.The tracer is released at all interfaces of the Guaíba River and its subsidiaries.As an initial condition, the Guaíba River is assumed to be spatially homogeneous and well-mixed, with a tracer concentration equal to 1 mg L −1 over its entire domain.
In the second scenario, the same initial conditions of tracer concentration are adopted over the Guaíba River.However, the release of a tracer remains at a constant rate of 5 mg L −1 only over the first 10 simulation hours.After the first 10 simulation hours, the tracer concentration released at the interfaces of the Guaíba River and its subsidiaries is instantaneously reduced to 1 mg L −1 .
Errors due to numerical diffusion were assessed by monitoring the tracer concentration of a computational cell located in the Delta of the Jacuí River.

Hydrodynamic module
In general, the hydrodynamic module of the IPH-UnTRIM2D model successfully reproduced water level observations in response to tidal and runoff daily inflows.Comparisons between observed and calculated elevations showed that IPH-UnTRIM2D model captured the trends of peaks and valleys throughout 150-day period (Fig. 4).
Of all the water level stations, Cristal and Ponta Grossa presented better agreement with the observed water levels once they were influenced by the downstream boundary condition.To measure the agreement between observed and calculated elevations, the Pearson correlation coefficient was calculated for each water level station, as given by Eq. ( 22).
where X and Y denotes the time series of observed and simulated water level of length size n, respectively.X and Y are the series averages.Table 2 shows the Pearson correlation coefficient computed at each gauging station.The accuracy of the numerical approximation given by the hydrodynamic module was compared with the continuity equation for the water budget of the Guaíba River.Since the hydrodynamic module calculates free surface variations at every time step, the incremental volume at time step t + 1 is estimated as the differences between free surface variations between t and t + 1. Theoretically, the incremental volume must be equal to the differences between incoming and outgoing water fluxes into the Guaíba River.
Results obtained from the simulations indicated that percentage differences between the incremental volume and incoming/outgoing water fluxes in the Guaíba River were lower than 1 % of the total volume of the system.As this error in volume conservation is lower than uncertainties associated with data collection procedures (Moradkhani and Hsu, 2005), it suggests that global (and also local) volume conservation was achieved.
According to the velocity fields (Fig. 5), the hydrodynamic behaviour of the Guaíba River was characterized by higher current velocities in the main channel (ca.0.1 m s −1 ) and    lower values along its shorelines (ca.0.005 m s −1 ).On the other hand, at the Delta of the Jacuí River, velocities and free surface elevations varied depending on the narrowing and widening effects of its set of stream channels.Although the hydrodynamic solution reproduced wetting and drying zones, these processes were not considered during the simulations due to the lack of topographic data along the shorelines of the Guaíba River and the Delta of the Jacuí River.

Numerical schemes for the advection-diffusion transport equation
Once the hydrodynamic solution was calculated, its velocities and elevations field were employed to solve the advection-diffusion equation using the higher-order fluxlimiter scheme.In the first scenario, accumulated mass balance errors for both first-order upwind and higher-order fluxlimiter schemes were computed at every time step of the simulation (Fig. 6).
Both schemes have showed that their accumulated mass balance errors were lower than 20 % over a period of 15 days after tracer was released.However, the higher-order fluxlimiter scheme used by IPH-UnTRIM2D led to a less diffusive approach compared with a first-order upwind scheme.At the end of the simulation, accumulated mass balance errors reached 18 % for a first-order upwind scheme, whereas the higher-order flux-limiter scheme yielded errors below 8 %.
In terms of tracer concentration, 1-day accumulated error in mass balance generated by the first-order upwind scheme ranged from −0.12 to 0.18 mg L −1 over the entire ecosystem.When derived from the higher-order flux-limiter scheme, these values fluctuated within a range of −0.11 and 0.13 mg L −1 .Since suspended sediment behaves as a passive tracer, the magnitude of these errors in mass balance is compared with estimates of daily sedimentation rates performed by Stevenson et al. (1985) and Ogston et al. (2004); Ogston and Field (2010), as well as measurements of sediment resuspension in estuaries (Bokuniewicz et al., 1991;Hill et al., 2003).
Sediment fluxes measured 0.3 m above the bottom in Long Island Sound (Bokuniewicz et al., 1991), the Mersey estuary and Dover Straits (Hill et al., 2003) (i.e.estuaries located along the shore, like the Guaíba River) varied from 0.00063 to 0.00100 mg cm −2 s −1 .Under the same sediment resuspension or sedimentation rates, the total particle gain or loss due to resuspension or sedimentation over a 1-day period in the Guaíba River may achieve up to 144 mg L −1 .Estimates of sediment resuspension and sedimentation in bays and estuaries (Stevenson et al., 1985;Ogston et al., 2004;Ogston and Field, 2010) showed that sediment resuspension and  sedimentation rates can fluctuate from 50 to 80 mg cm −2 d −1 throughout the year that leads to minimum and maximum daily deposition rates of 83 and 133 mg L −1 , respectively.Therefore, both measured and estimated particle gain or loss due to sediment resuspension or sedimentation were by far larger than error in mass balance generated by the numerical schemes.
Differences between the diffusion of higher-order fluxlimiter scheme and first-order upwind scheme were evaluated by plotting the evolution of tracer concentration fields over 26, 130 and 312 h (Fig. 7).Although both numerical solutions are similar after the first 130 h of simulation, the higher-order flux-limiter scheme presented faster spreading of tracer concentration over pelagic zones (open water).After 312 h, tracer concentrations reached Patos Lagoon, which indicates that the Guaíba River has a residence time approximately equal to 13 days, in accordance with (Rosauro, 1982;Silveira, 1986).
In the second scenario, after 15 days, accumulated mass balance errors do not exceed 3 % for both schemes (Fig. 8).Moreover, a slight and constant difference between their mass balance errors was lower than 0.5 % for whole simulation, where the higher-order flux-limiter scheme proved to be more conservative than a first-order upwind scheme.As the higher-order flux-limiter scheme may employ either first-or second-order approaches, it leads to a less diffusive solution compared with the solution calculated using a first-order upwind scheme.By analysing the pulse tracer shapes, the higher-order flux-limiter scheme showed slight spreading of tracer and peak tracer concentrations of close to 5 mg L −1 , whereas a first-order upwind scheme presented a higher diffusing capacity that smoothes gradients of tracer concentration (Fig. 9).

Conclusions
This paper presented comparisons between two different numerical schemes for solving advection-diffusion equations on unstructured grids using IPH-UnTRIM2D.As shown in previous applications of unstructured grid hydrodynamic models (Casulli and Walters, 2000;Zhang and Baptista, 2007), IPH-UnTRIM2D was willing to represent the time series of water level observations and the movement of suspended material throughout the Guaíba River.Its application to the Guaíba River showed that unstructured grids presented higher flexibility in representing the shape of the Guaíba River than previous studies on uniform grids (Rosauro, 1982;Silveira, 1986) once unstructured computational cells varied depending on bathymetry, geometry and shoreline-fitting boundary of the Guaíba River.
The efficiency of the higher-order flux-limiter scheme was tested by comparing its solution with a first-order upwind scheme solution for two scenarios in the Guaíba River.According to comparisons, the higher-order flux-limiter scheme was more mass-conservative than a first-order upwind scheme once accumulated mass balance errors achieved 18 % after a 15-day period for a first-order upwind scheme while flux-limiter errors were below 8 % in accordance with with numerical analysis performed by Casulli and Zanolli   (2005) in a curved channel under controlled boundary conditions.Moreover, accumulated mass balance errors showed that, independently of the numerical scheme employed, mass conservation was proportional to the total amount of mass released in the system.Although both schemes have presented mass conservation errors, these errors are assumed negligible  compared with losses due to erosion, sedimentation or decay of a substance in estuaries (Stevenson et al., 1985).
Ongoing studies have been developed to include the higher-order flux-limiter scheme on a three-dimensional complex dynamic model for aquatic ecosystems (Fragoso Jr. et al., 2008;Fragoso Jr et al., 2009) in order to create an efficient and capable tool for performing analysis of long-term ecosystem dynamics.

Fig. 1 :
Fig.1: Illustration of an unstructured orthogonal grid (solid lines).An auxiliary grid (dashed lines) composed of segments that join the centres of each computational cell is also shown to highlight geometric variables used in the section 2.9

FFig. 2 :
Fig. 2: Location of the River Guaíba and municipalities on its banks.It also shows the composed of a group of 5 water level stations (circles) and 2 discharge stations (triangles).

Fig. 2 .
Fig. 2. Location of the Guaíba River and municipalities on its banks.It also shows the composition of a group of 5 water level stations (circles) and 2 discharge stations (triangles).

Fig. 3 :
Fig. 3: Spatial discretization of the river Guaíba using an orthogonal triangular grid (a) and its bathymetry (b).

Fig. 3 .
Fig. 3. Spatial discretization of the Guaíba River using an orthogonal triangular grid (a) and its bathymetry (b).
Fig. 4: Calculated and observed elevation of water surface at four control points along the River Guaíba: (a) Ilha da Pintada; (b) Ipanema; (c) Cristal and (d) Ponta Grossa.

Fig. 6 :
Fig. 6: Accumulated mass balance errors for the flux limiter and upwind schemes over 15 days of simulation.

Fig. 6 .
Fig. 6.Accumulated mass balance errors for the flux-limiter and upwind schemes over 15 days of simulation.

Fig. 7 :
Fig. 7: Concentration field of a tracer plume calculated using the flux limiter and upwind schemes for 26, 130 and 312 hours in the River Guaíba where red computational cells denote values of concentration higher than 4.5 mgL −1 while blue cells mean values of concentration lower than 2.0 mgL −1 .

Fig. 7 .Fig. 8 :
Fig. 7. Concentration field of a tracer plume calculated using the flux-limiter and upwind schemes for 26, 130 and 312 h in the Guaíba River, where red computational cells denote values of concentration higher than 4.5 mg L −1 , while blue cells mean values of concentration lower than 2.0 mg L −1 .F. F. Pereira et al.: Assessment of Numerical Schemes

Fig. 8 .
Fig. 8. Comparisons between accumulated mass balance errors considering a pulse input of tracer concentration into the Guaíba River during 10 h for the flux-limiter and upwind schemes.

Fig. 9 :
Fig. 9: Pulse of a passive tracer at the Delta of the River Jacuí.

Fig. 9 .
Fig. 9. Pulse of a passive tracer at the Delta of the Jacuí River.

Table 1 .
Streamflow and water level gauging stations of the Guaíba River and its tributaries used in this study.

Table 2 .
Comparisons between observed and simulated water elevation data at the gauging stations (R-squared).