An Analysis of Total Phosphorus Dispersion in Lake Used As a Municipal Water Supply

In Belém city is located the potable water supply system of its metropolitan area, which includes, in addition to this city, four more municipalities. In this water supply complex is the Água Preta lake, which serves as a reservoir for the water pumped from the Guamá river. Due to the great importance of this lake for this system, several works have been devoted to its study, from the monitoring of the quality of its waters to its hydrodynamic modeling. This paper presents the results obtained by computer simulation of the phosphorus dispersion within this reservoir by the numerical solution of two-dimensional equation of advection-diffusion-reaction by the method θ/SUPG. Comparing these results with data concentration of total phosphorus collected from November 2008 to October 2009 and from satellite photos show that the biggest polluters of the water of this lake are the domestic sewage dumps from the population living in its vicinity. The results obtained indicate the need for more information for more precise quantitative analysis. However, they show that the phosphorus brought by the Guamá river water is consumed in an area adjacent to the canal that carries this water into the lake. Phosphorus deposits in the lake bottom should be monitored to verify their behavior, thus preventing the quality of water maintained therein.

this assertion is the Metropolitan Region of Belém (MRB) in the state of Pará, in northern Brazil, as can be seen in Figure 1.
The MRB is one of the largest metropolitan areas in the country with a total area of 1,819 km 2 .The 2010 census released by the Brazilian Institute of Geography and Statistics (IBGE), indicated that approximately 2,197,807 people reside in this region, which corresponds to approximately 30% of the population of the state of Pará.
The hydrological supply complex of the MRB's is composed of the Guamá river, along with the Água Preta and Bolonha lakes.Figure 2 shows that this supply system begins with the capture of water from the Guamá to maintain the levels of the lakes.Through pipelines, these waters are then dumped in Água Preta lake, and then are taken to the Bolonha lake through an artificial channel.Following this step, waters are conducted to the water treatment plant (WTP) Bolonha, and then distributed to the MRB.This complex is located in Belém city.
Água Preta lake is considered the main lake which serves as a source of supply of the MRB (Lima et al. 2014).It is formed by the Catu and Água Preta streams and has three springs.It has a free surface area of approximately 3,116,868 m 2 and estimated volume of 10,550,000 m 3 with maximum depth of approximately 8.5 m.The output control of its waters to the Bolonha lake through the channel is done by gravity effect with a system of locks.Vasconcelos and Souza (2011) indicate that the waste produced in MRB is dumped, untreated, directly into the Guajará Bay and Guamá river, which in turn has its water transported to the wellspring of Utinga where the Água Preta lake is.In addition, the resident population on the outskirts of this pond also contribute directly to the increased pollution of the waters, since it not having basic sanitation, dump their sewage SIMULATION OF THE PHOSPHORUS TRANSPORT IN A REAL LAKE directly into the waters dammed in this source (Holanda et al. 2011).Because of this, the Água Preta lake is not exempt from the problem of eutrophication, as can be seen in Figure 2, probably due to the proximity of inhabited areas.In this figure one can see the water supply complex.Figure 3 shows the location of three selected points for analyzing the behavior of velocity field within the domain; the first is over the entry contour, the second point is located inside the path of the flow region between upstream and downstream and the last point is situated at the outflow contour.Additionally, three localized sources were inserted, representing the discharge of wastes derived from domestic sewage, as indicated also in Figure 3.The approximate locations of those sources were obtained during visits to the site and their concentration variations were modeled by rainfall distribution, since these variables are correlated.In order for this information to be utilized in each time step, it was continuously interpolated, both in terms of the contour and the sources, using the cubic spline method during processing.Therefore, due to the relevance of this reservoir, monitoring and management of the water in this lake is extremely important for the wellbeing of the population it serves.One of the possible tools to aid this management lies in the computational simulation of pollutant transport in aquatic environments.In this context, this paper examined total phosphorus transport within Água Preta lake, through a computerized simulation using data collected in loco over a twelve-month period, and generated a hydrodynamic models.These were used for numerical solution, with the θ /SUPG method, of a simple model of an advection-diffusion equation in the presence of a reaction term.This method consists in solving the mathematical model in two parts.The first one is the temporal component, this is done using the general formulation of the finite difference (the θ family of methods) for integrating first order differential equations (Donea and Huerta 2003).Several methods are obtained for different values of the θ parameter.For values of θ < 1/2 these schemes are conditionally stable, the Euler method, obtained for θ = 0, is the well-known.For values of θ ≥ 1/2 these schemes are unconditionally stable, the Backward Euler (θ = 1) and Crank-Nicolson (θ = 1/2), are the most usual ones.The second part of the solution consists in the discretization of the spatial components for the stabilized Finite Element method called Streamline-upwind Petrov-Galerkin (SUPG).This is done due to the deficiencies of the standard Galerkin formulation in solving convective-dominated problems.For a detailed explanation of the θ/SUPG method we suggest (Lima et al. 2014, Donea andHuerta 2003).

GOVERNING EQUATIONS AND NUMERICAL METHOD
Due to the horizontal dimensions (width and length) of the lake being greater than its vertical one (depth), a bidimensional model is the most suitable to simulate the phenomena of the transport of a substance in a fluid.So, for the computer simulation of phosphorus transport, we used Equation (1), known as 2D advectiondiffusion-reaction equation.It represents the distribution of the concentration of a matter inside a domain as a function of time and space.
it is subject to initial condition presented in Equation (2) (Marchuk 1986, Loucks et al. 2005, Tsanis et al. 2006, Ji 2008): where, in the present study, C is the phosphorus concentration, u is the bidimensional velocity field, D is the diffusion coefficient, t is the temporal variable, k is the reaction coefficient and the term S represents the sources.In this work dispersion and reaction coefficients are assumed to be constant in time and space.The values adopted for these coefficients were obtained by adjusting of the simulation to the real data.
To obtain a well posed problem it is necessary to impose conditions on the contour of domain, denoted by Γ.At solid boundaries no flux must be adopted, this can be expressed as a second type boundary condition: where n is the normal vector on the solid boundaries.On the free flow boundaries of the domain, it is necessary to impose the condition at inflow and outflow contour.To the inlet boundary, and the other entrances, it was assumed that the value of the scalar was known, in other words, Variable g in this equation represents the concentration that came from the sources (river and sewage), so it is a function of time.For the water from the river, g is the polynomial interpolation through the data collected at this point.In sewages, no concentration data was taken, so we assumed them as percentages of the rainfall, since they are directly correlated.
The boundary condition for the outlet section is represented by a mass balance equation of the form (Massabó et al. 2011) where B is the outflow boundary, and in and out means, respectively, the flux inside and outside the domain.However, C x 1 = out is not known a priori.A way of solve this problem is to assume that the concentrations are continuous across this boundary (Massabó et al. 2011), namely, Substituting this in Eq.( 5) results in (6) The system formed by Eqs.( 1)-( 4) and ( 6), is the problem to be solved by the θ /SUPG-method.In this paper we first explained the discretization of the problem in temporal variable.However, during the computational implementation the spatial discretization was chosen first, since it was assumed that the mesh does not change during the temporal evolution.This equation was discretized by the generalized method of finite differences known as θ family, for the temporal component and by the SUPG method of finite elements for the spatial components (Donea and Huerta 2003), which produces the following linear system (details can be found in Lima et al. 2014), where θ ∈ [0,1] is the parameter that determines the scheme of finite differences that will be used, Δt is the temporal discretization interval, n is the rate of times, α = τ∆t/2, with τ being the stabilization parameter in the SUPG method, given by Donea and Huerta (2003), John and Schmeyer ( 2008): (8) where h is the local size of the mesh element, in which N is the vector of the base functions for interpolating finite elements, C is the vector of the nodal values of the concentration, F is the vector for the sources and RÔMULO C. LIMA et al.
For simulating the transport of total phosphorus, from the waters of the Guamá river, the data in the input model were used as boundary condition during temporal evolution.To build this dynamic simulation, using monthly averages for height and flow through the pipeline that brings water from the river, and so generates velocity fields for each month of the period analyzed, we used the Modeleur/Hydrosim package, which consists of the joining of geographical information discretized in two dimensional quadratic finite elements (Modeleur) in addition to a module to solve the 2D hydrodynamics equations (Hydrosim) (Secretan and Leclerc 1998).

EXPERIMENTAL INPUT DATA
For the computer simulation of the dispersion of total phosphorus into the Água Preta lake, groups of data collected on site or in the region where the lake is located were used.The first of these groups was data concerning the region's rainfall variation, which is primarily responsible for the reservoir's dynamics, since it determines the need for the increase or decrease in the volume of water pumped from the Guamá river, as well as the variation of concentration of pollutants from external sources.
Figure 4 shows the rainfall average for the months of November 2008 to October 2009.One can observe that the period from December to May presents rainfall average above 150 mm, reaching values around 450 mm in March, thus characterizing the local rainy season.On the other hand, the interval from July to November has average values for accumulated precipitation below 150 mm, representing the dry season.
The second set of data, refers to the hydrological variables in Água Preta lake, height of the water level and pumping flow of water from the Guamá river to the lake, as shown in Figure 5.One may note the importance of rainfall precipitation in these variables, given that on the months with lower rainfall it is necessary to increase the pumping flow of water from the Guamá in order to maintain the water level.This can be seen mainly during the first four and last three months of collection, which indicate below-average precipitation for those periods, and thus, higher values for water flow being pumped.On the other hand, the period from March to July 2009 recorded average rainfall above normal, notably in the months of March and June, which have an accentuated reduction in the flow.However, this period marks the transition from the rainy to the dry season, and so there is a gradual increase in the flow pumped, from a minimum value of 3.052 m 3 /s until it again reaches higher numbers.
The third group of date reports the concentration of total phosphorus collected at two points on the lake, at the exit of pipelining of pumped water from the Guamá river (entry model) and one situated at the entrance of the canal that carries water to the lake Bologna (model output), as indicated in Figure 3.The samples were collected on a single day each month during the period from November 2008 to October 2009.The information obtained at these points is presented graphically in Figure 6.
Unfortunately, due to the growth of aquatic plants in the lake, collecting data was difficult.Thus, the points shown are the only two where it was possible to collect information on the concentration of phosphorus into the lake over the entire period.
It is possible to observe the influence of the flow of water pumped from the Guamá river on phosphorus concentration recorded at the entry of the Água Preta lake, since a comparison of these two sets of information indicates the same behavior.For example, from November 2008 to February 2009, there was an increase in the admission of river water, going from approximately 5 m 3 /s to 6 m 3 /s, with the concentration following the same tendency and rising from 0.08 mg/L to 0.29 mg/L.Later, due to the rainier period, there was a reduction in water being injected into the lack, which during the first four months went from an average of 5.5 m 3 /s to slightly more than 4.0 m 3 /s.From March to July 2009, the average concentration showed a decline from 0.65 mg/L to 0.188 mg/L, and then reached its lowest value (0.06 mg/L).Next, during the final quarter of data collection, the flow again rose, climbing to an average of 5.1 m 3 /s, and the concentration showed the same tendency for growth in this period; however, its average was 0.177 mg/L, due to the low values recorded for this period when compared to the previous four months.
Cubic spline interpolation was applied in concentration of phosphorous at the input and output of the water in the lake, Figure 6.
With regards to the data obtained in the outlet part of the model, one may note the presence of higher values of the concentrations than in the months of January, February, May and June 2009, with a greater difference between the concentrations upstream and downstream, which occurred in the month of May, and was equal to 0.17 mg/L.This behavior indicates the existence of variable sources in time, and possibly in space as well.As previously indicated, the dynamics of the water within the lake suffers great influence of rainfall precipitation, modifying the volume of water that must be pumped from the Guamá river.Each of these fields was maintained invariant on the simulation period of 30 days, in this way, the simulation of transport can be considered "semi-instationary", when considering the entire period of modeling.
The simulated velocity fields present a small interval between their maximum and minimum values, as shown in Figure 7, which shows the model used at the beginning of the simulation, obtained with the flow and height of water from November 2008.In this case, the velocity at the entrance of the model is approximately equal to 0.031 m/s, while at the outlet this value is slightly below 0.24 m/s.Furthermore, the model presents a flow with an average velocity of 0.01 m/s in the western portion of the model, as well as a large area where there are recirculation, located near the outline of admission of the waters of the Guamá river.On the other hand, the eastern portion of the model has very low, approximately null, velocities.
The image on the left of Figure 7, is the vector field and the current lines (light curves), while the image on the right represents the velocity module for seepage.In this work only this velocity model was showed for brevity (for more details, please see Lima et al. 2014).Due to the very small dynamics of the lake, there are no notable changes among the velocity fields.
Figure 8 shows the flow velocity at the points indicated in Figure 3.It is possible to note that the velocity at the outlet of the pipeline shows small variations throughout the hydrodynamic simulation between consecutive months, with the largest difference between the maximum and minimum values equal to 0.02 m/s on the months of February and March 2009.The outlet of the model is the site that presents the greatest variations in velocity between adjacent months, with the difference being most accentuated for the period of March to April 2009, when the velocity went from approximately 0.26 m/s to 0.2 m/s.However, SIMULATION OF THE PHOSPHORUS TRANSPORT IN A REAL LAKE the maximum velocity obtained was 0.27 m/s in February 2009, and the minimum was slightly below 0.2 m/s in June 2009.With regard to velocity at the internal point in the model, there is an average of 0.006 m/s, with the greatest variation in its module at 0.0032 m/s, from March to May 2009.The initial condition for distribution of phosphorus, presented in Figure 9, was produced by simulation, assuming that, initially, there was no phosphorus concentration in the lake area.So we adopted constant injection from the sources (river (0.08 mg/L) and sewages (2.164 mg/L, 1.701 mg/L and 1.712 mg/L, respectively)), corresponding to phosphorus data for November 2008.Due to the velocity field being almost null in the right branch of the lake, sewage dumps located on it did not have influence on the phosphorus concentration in the water that leaves the lake, so they were not included in the simulation.The velocity field was also kept constant, with reference to that month, and diffusion coefficient was kept equal to 0.09 m 2 /s.The simulation took place until it reached the steady-state at one point in the output, as can be seen in Figure 10, which occurred in a range of 87 h (approximately 3.5 days).Furthermore, a value of 3x10 -6 day -1 was adopted for the reaction coefficient.One can see that, in this scenario, the phosphorus from the Guamá river has no direct influence on the concentration that abandons the reservoir, due to reactions within the lake.
Figure 10 presents the behavior curve of concentration versus time, obtained at a point on the model output during the simulation to get the initial condition.Note that the element reaches the exit after 10 h.From there, the concentration increases until it reaches the stationary condition in approximately 0.044 mg/L, which occurs after 50 h, generating the initial model shown in Figure 9.The distribution of total phosphorus in Água Preta lake from November 2008 to October 2009, is depicted in Figure 11.Adopted coefficients of diffusion and reaction are the same used for getting the initial condition, i.e., 0.09 m 2 /s and 3x10 -6 day -1 , respectively.In addition, the scale of the concentration ranges from zero to the maximum value on the output of the model, approximately 0.38 mg/L, such that the behavior of the element can be evaluated in this contour.Note that the main sources of pollution are the point sources in the vicinity of the channel connecting the two lakes, since, the largest part of the concentration resulting from the Guamá river is consumed before it reaches this channel.The greatest contribution of the river for the phosphorus concentration occurs on March 2009 (120 days), due to the higher incidence rainfall, that rise the values of the phosphorus transported by river, as well as rises the velocity field inside the reservoir, causing a greater amount of this substance than the consumption capacity of the reactions.Even so, this contribution is of the order of 0.1 mg/L, which corresponds to nearly one third of the value obtained at the output for that month.
The advective effect can be noticed due to the presence of phosphorus in the west portion of the model and in the recirculation area, where speeds are higher, which causes the drag of this element before its consumption.Another area where this effect can be observed is the neighborhood of the output where there is an acceleration of the waters, causing the conduct of phosphorus from the sources 1 and 2 for the channel between the lakes.On the other hand, the diffusive effect is indicated by the spread of phosphorus in regions where speeds are lower, as, for example, in the surrounding areas of the stream connecting the input and output of the model and, mainly, in the eastern portion near the entrance to the lake.Also, it is due to diffusive effect, the conduction of the phosphorus from source 3, once the velocity field in that region is approximately null.
The reactive component corresponds to all of the forms of the phosphorus consumption, from its interaction with sediments and deposition on the botton, to its consumption by primary producers.The reaction rate, considered constant, does not reflect reality, given that this coefficient is a function of the temperature in the environment.However, it is possible to use it as an average value for analyzing the behavior of the simulation results.That being the case, as had been observed earlier, a large share of the element derived from the Guamá river is depleted before reaching the downstream portion of the model, and its lowest contribution occurs in June 2009 (240 days), when the injection of phosphorus by the river in the lake is smallest.Furthermore, this period recorded the lowest concentration value at the outlet, given that the source points also present the minimum values of phosphorus insertion.Silva (2010) affirms that despite the strong environmental impacts to which the catchment system and its associated lake are subjected, it was observed that it maintains its capacity for self-purification and renewal.Thus, the results of this simulation are in agreement with that statement, given that most of the lake area presents low values of concentration during the majority of the period simulated.Furthermore, observations of these results indicate that the areas with the highest rates of this element throughout the period and domain modeled, are due to the sewerage system.This is in agreement with the satellite image presented in Figure 2 (or Figure 3), indicating that this result is correct, at least in qualitative terms, given that these areas are covered by floating aquatic plants.
Figure 12 represents the comparison of the real observations (points) versus computer simulation (line).One may see some discontinuities that cause breaks in the smoothness of the curvilinear line.That is caused by the abrupt passage of one velocity field to the next, despite this, one may note that the curve have the tendency of the data, indicating that the computational model recovers the information about the concentration at the assessment point.However, this same curve suggests the need for more knowledge in order to have a refined calibration of the result, such as enhanced models of production and consumption of the total phosphorus components within the lake area, more precise data on the injected discharges, etc.
One may affirm that the greatest differences obtained among these data lie in December 2008 and May 2009, the latter one is due to the value out of the trend in the real data.Indeed, it was an unusual concentration for that month, caused for unknown sources, and calculated solution shows the normal trend.As noted before, this may correspond to an anomalous variation in the sources existing in the lake.In fact, one possible explanation for that behavior is the appearance of a diffused source in the lake area, such as resuspension of phosphorus deposited on the bottom or decomposition of organic matter.To conclude, a realistic simulation contaminant transport requires the description of several complex relations among the components in the environment in which this element is being conducted, as well as the need for collecting specific data.On the other hand, one should take into account the balance between the computing cost, the goal of the study and the data available for such a simulation.That being the case, the test presented in this paper indicates that the θ/SUPG method applied in the numerical solution of the simple advectiondiffusion-reaction equation, that serves as the basis for the transport models in fluid environments, may be employed in a qualitative assessment of the dispersion of pollutants in aquatic environments at a reasonable computational cost.To obtain better results for simulated data, it was necessary more information inserted in the mathematical model, for example, the mass variation of dead material since it produces phosphorus inside the lake, the amount of living aquatic plant since it absorbs the element from the lake, etc.However, these data were not available.

CONCLUSIONS
This paper analyzed the behavior of the transport of phosphorus in Água Preta lake, one of the reservoirs that is part of the water supply system in the metropolitan region of Belém, Amazonia, Brazil, through computer simulation obtained by numerical solution of the advection-diffusion-reaction equation, which is the basis of the model of pollutant transport as well as other transport phenomena in hydrodynamic environments.Particularly, the solution consisted of applying the SUPG finite elements method for spatial discretization and a generalized method of finite differences, known as θ family, to integrate the temporal component of the equation.In this manner, this procedure was applied to real data collected during November 2008 up to October 2009.Moreover, in this simulation the velocity fields were generated with information on the height of the water level and flow at the lake entrance for each month in the aforementioned time interval and all of them were considered to be invariable in the 30 day time interval.The phosphorus concentration, for its part, had measurements taken at the entrance and outlet of the lake.From that point, using the data upstream of the model, together with hypothetical information for source points, phosphorus transport downstream was simulated, applying a constant rate of consumption of this element within the domain.As a consequence, the result obtained indicates the need for denser information for a more accurate quantitative analysis of the process in question.However, they indicate that phosphorus inserted by the Guamá river is consumed in an area near the upstream portion, before it can reach the reservoir outlet.Furthermore, the simulation indicated that the main supplier of phosphorus in the water, leaving the lake, is domestic sewage discharged by the population living close by.This result is in agreement with images obtained via satellite, given that in this area there is a great abundance of floating aquatic plants, which indicates the existence of abundant nutrients, including phosphorus.To conclude, the deposits of the phosphorus at the bottom of the lake must to be monitored to verify their future behavior, to prevent the deterioration of the water quality inside of it.RÔMULO C. LIMA et al.

Figure 1 -
Figure 1 -Metropolitan Region of Belém (MRB) location within the state of Pará.At the bottom right the location of this state on the map of Brazil.Source: Adapted from Wikipédia 2012.

Figure 2 -
Figure 2 -Water supply complex of the metropolitan region of Belém.

Figure 4 -
Figure 4 -Average rainfall in the Água Preta lake area for the period from November 2008 to October 2009.Source: SECTAM 2012.

Figure 5 -
Figure 5 -Height of the water level in the lake (H) and flow of pumping of water from the Guamá river (Q).

Figure 6 -
Figure 6 -Graphic representation of the concentration values for total phosphorus (mg/L) at the entrance and exit of reservoir waters.

Figure 7 -
Figure 7 -Velocity field generated by Hydrosim hydrodynamic simulator, to the left, and its module.TOTAL PHOSPHORUS DISPERSION

Figure 8 -
Figure 8 -Velocity module at the three sampling points.

Figure 9 -
Figure 9 -Initial condition for the distribution of the total phosphorus.

Figure 10 -
Figure 10 -Temporal evolution of the concentration in the output of the model.

Figure 11 -
Figure 11 -Temporal evolution of the transport simulation of phosphorus.