Reliability Engineering and System Safety 79 (2003) 239–244 www.elsevier.com/locate/ress
Sensitivity analysis as a tool for the implementation of a water quality regulation based on the Maximum Permissible Loads policy R. Pastresa,*, S. Ciavattaa, G. Cossarinib, C. Solidorob a
Dipartimento di Chimica Fisica, University of Venice, Dorsoduro 2137, Venezia 30123, Italy b Osservatorio Geofisico Sperimentale, P.O. Box 2011, 34016 Opicina, Trieste, Italy
Abstract This paper shows how local sensitivity analysis, in respect of the parameters which specify the boundary conditions, can be used for relating the total load of non-conservative pollutants to their distributions within a water body. The method is applied to the estimation of the Maximum Permissible Load of inorganic nitrogen in the lagoon of Venice, that is of the maximum load of nitrogen which keeps its average yearly concentration below a prescribed threshold. The use of the spatial distributions of sensitivity coefficients in order to rank the sources of pollution and to forecast the effect of a reduction in the pollution is also discussed. q 2002 Elsevier Science Ltd. All rights reserved. Keywords: Sensitivity analysis; Maximum Permissible Loads; Dissolved inorganic nitrogen; Lagoon of Venice
1. Introduction Since 1999, in Italy the pollutant loads which enter the water-bodies have been disciplined by means of the socalled maximum-permissible-loads (MPLs) policy. Within this framework, Local Authorities should make an inventory of the sources of pollution and then fix the level of emission of each activity, in order to maintain the concentrations of potentially dangerous substances below prescribed thresholds, called ‘Quality Targets’ (QTs). The implementation of this policy may clearly benefit from the use of mathematical models, which can be used as tools for both estimating the MPLs, by solving the so-called’ inverse problem’, and exploring the consequences of different input scenarios. In fact, in mathematical terms, the loads are specified by a set of boundary conditions: numerical models can then be used for determining a functional relationship between the set of parameters which specify the boundary conditions and the output variables which one decides to compare with the QTs. Once this task has been accomplished, one can invert this function in order to estimate the MPLs which are compatible with the targets. These problems are investigated in this paper using the lagoon of Venice as a case-study and the sensitivity analysis * Corresponding author. Tel.: þ 39-41-2578-674; fax: þ39-41-2578-594. E-mail address:
[email protected] (R. Pastres).
in respect of each source of pollution as a tool. In fact, because of its peculiarity, this system was thoroughly investigated and a 3D reaction –diffusion water quality model is already available [1]. Advective transport is not included in the model, as numerical simulations suggested that the residual currents in the lagoon of Venice are very small [2 – 4]. However, diffusivities embody information about the contribution of the tide to the dispersion, as they are computed at each grid point on the basis of a statistical analysis of the results of a Lagrangian particle dispersion model [5]. The reaction– diffusion Eq. (1) is solved using a finitedifference scheme
›cðx; y; z; tÞ=›t ¼ 7ðKðx; y; zÞ7cðx; y; z; tÞÞ þ fðcðx; y; zÞ; b; tÞ ð1Þ In Eq. (1), c is the state vector; K, the tensor of eddy diffusivities; f, the reaction term; b is the set of site-specific parameters. The model simulates the dynamic of the ecosystem up to the second trophic level by using 12 state variables. The state vector includes the concentrations of the two main forms of inorganic nitrogen, ammonium and nitrate, as well as the inorganic reactive phosphorous one. These chemicals are considered to be the main cause of the eutrophication and, therefore, the current legislation fixes their QT for the lagoon of Venice. In this paper, we have applied the method outlined below to the estimation of
0951-8320/03/$ - see front matter q 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 9 5 1 - 8 3 2 0 ( 0 2 ) 0 0 2 3 5 - 1
240
R. Pastres et al. / Reliability Engineering and System Safety 79 (2003) 239–244
the MPLs of ammonium and nitrate: the sum of their concentration gives the dissolved inorganic nitrogen (DIN), as the concentrations of nitrite is very low. At present, the concentration of DIN is above the target, while the concentration of reactive phosphorous is very close to it, as its use in detergents was prohibited in 1989. Ammonium, NH4, and nitrate, NO3, are carried into the lagoon by the rivers, and are directly released from the Industrial area of Porto Marghera, on the edge of the lagoon, and from the city of Venice and the nearby islands. The yearly evolutions of these inputs were modelled using Von Neumann-type time-dependent boundary conditions: the fluxes Fi were specified using a set of trigonometric polinomia 4 4 FNH ðtÞ ¼ aNH i i;0 þ
3 X
4 ½aNH i;2j21 cosð2pjt=365Þ
j¼1
£
4 aNH i;2j
senð2pjt=365Þ
3 3 FNO ðtÞ ¼ aNO i i;0 þ
3 X
ð2Þ
3 ½aNO i;2j21 cosð2pjt=365Þ
j¼1
£
3 aNO i;2j
senð2pjt=365Þ
ð3Þ
3 4 where FNH ðtÞ and FNO ðtÞ are the daily fluxes of i i ammonium of nitrate released by the ith source at time t, expressed in days. NO3 NO3 NH4 4 The parameters aNH were i;0 ; …; ai;6 ; …; ai;0 ; …; ai;6 estimated for each source by means of a least squares regression of the monthly data reported in Ref. [6] and on total loads data presented in Ref. [7]. The exchanges with the Adriatic sea at the three inlets were described by means of Dirichlet-type boundary conditions: the concentrations of ammonia and nitrate at the boundaries were taken from Ref. [8]. On the basis of the available data, it
was possible to define 16 sources, which are shown in Fig. 1 together with the sampling stations which were monitored from 1995 to 1999. Therefore, 2 £ 7 £ 19 parameters had to be considered as potential input factors, if one wishes to take into consideration also the fluxes through the three inlets.
2. Method Global sensitivity analysis methods, based on Monte Carlo techniques [9,10], could hardly be used in this case for exploring the dependence of the state of the system on the variation of the set a of parameters which specifies the loads, because each run of the model requires about 5400 s to simulate 1 year on a 2-CPU Digital AU533MHz WS. Therefore, we decided to use local sensitivity analysis in order to † estimate the role of each source regarding the determination of the concentration of DIN in a given portion of the water body; † estimate the MPL for DIN, MPLDIN, given the water quality target, QTDIN; † obtain a first approximation of the evolution of the system within the MPLDIN input scenario; † analyse the consequences of different policies concerning the reduction of the load of nitrogen. In order to simplify the problem, we supposed that the QTDIN had to be respected for the yearly average concentration of the system: this interpretation of the law was arbitrary, but, at the moment, no precise statement exists about the type of average which has to be compared with the QT. In this case, one has to take into consideration only the NO3 NO3 NO3 NH4 NH4 4 first parameters aNH 1;0 ; …; ai;0 ; ans;0 ; a1;0 ; …; ai;0 ; ans;0 ; for the ns manageable sources ðns ¼ 16Þ: To be concise, in the following development the second subscript index is dropped and the parameters are regrouped in a single vector, 3 by setting aNO ¼ ansþi : i The sensitivity equation for a generic parameter ai reads as
›Si ðx; y; z; tÞ=›t ¼ 7ðKðx; y; zÞ7Si ðx; y; z; tÞÞ þ JSi
Fig. 1. Sources of DIN, white circles, and the monitoring network, filled circles.
ð4Þ
where Si ¼ ›c=›ai is the sensitivity vector, which components are the sensitivities of each state variables with respect to the parameter ai ; and J ¼ ›f=›c is the Jacobian matrix of the vector function f. The sensitivity vector depends on the choice of the reference state around which the linearization is performed, as the reaction term of Eq. (4) is non-linear. The set of 2 £ ns vector Eq. (4) is solved by means of the direct method [11]. The partial derivatives which form the Jacobian matrix are calculated using symbolic calculus: this may appear to be a limit, with regard to extending this approach to other problems, but such calculations are now performed automatically by
R. Pastres et al. / Reliability Engineering and System Safety 79 (2003) 239–244
241
a number of software packages, which also give the corresponding piece of Fortran code as an output. The estimation of the MPLDIN is based on the linearization of the state Eq. (1): this hypothesis is certainly a critical one, but it allows one to obtain a quick first guess. Under this assumption, the effect of the simultaneous variation of all the parameters ai ; on the state vector can be expressed as: Dcðx; y; z; tÞ <
2£ns X
Si ðx; y; z; tÞDai
ð5Þ
i¼1
Eq. (5) provides an estimation of the effect of the variation P of the yearly load of DIN, which is given by the sum 2£ns i¼1 Dai ; on each state variable. In accordance with Eq. (5), the variation of DIN ¼ NH4 þ NO3 ; can be estimated as: DcDIN ðx; y; z; tÞ <
2£ns X
4 3 ½SNH ðx; y; z; tÞ þ SNO ðx; y; z; tÞDai i i
ð6Þ
i¼1
The effect on the average yearly concentration, DINav is obtained by integrating over the spatial and time domains DDINav < ð1=VÞð1=TÞM
2£ns X
Dai
ð ð T
i¼1
V
4 ½SNH ðx; y; z; tÞ i
3 ðx; y; z; tÞdV dt þ SNO i
ð7Þ
which can be related to our problem and rewritten in a more compact form DDINav < QTDIN 2 DINav ðaÞ <
2£ns X
NO3 4 ½SNH av;i þ Sav;i Dai
ð8Þ
i¼1
Ð Ð iv where Siv av;i ¼ ð1=VÞð1=TÞ T V Si ðx; y;z; tÞdV dt; iv ¼ NH4, NO3 and DINav(a) is the average concentration obtained by using the current estimates of the loads a. Eq. (8) can then be used as a constraint to assess the MPLDIN for the 2 £ ns unknowns: in order to obtain a unique solution, one has to state in mathematical form the other f ¼ ð2 £ ns 2 1Þ constraints.
Fig. 2. Yearly evolution of the average concentration of DIN in the monitoring network of Fig. 1: comparison between model output and observed values.
store only the maps of the yearly average sensitivities S, with respect to the 16 input parameters. The results of this simulation are compared in Fig. 2 with a set of field observations collected during the years 1995– 1999 in the sampling sites shown in Fig. 1 [12]. The figure shows the yearly evolution of the average DIN concentration in the monitoring network which was estimated using the model, continuous line, and using the field data, dots, which were collected each month. This comparison demonstrates that the model is capable of simulating the yearly evolution of DIN in the central part of the basin both at a qualitative and at a quantitative level. In particular, the model correctly reproduces the decrease in DIN which occurs during the summer months and which is due to a decrease in the river discharges and to micro- and macro-algae uptake. An example of the spatial distribution of the annual average sensitivities ST;i ðx; y; zÞ; (T ¼ 1 year) which are used as inputs by a post-processing module, is shown in Fig. 3. The map represents the relative sensitivities, normalised on QTDIN, of the most important source: sT;i; ðx; y; zÞ ¼ ðST;i ðx; y; zÞai Þ=QTDIN
3. Results The 3D model described in Section 1 was run using the current estimates of the loads of nitrogen, thus obtaining the spatial distributions of the yearly average concentrations of all state variables: in the following, such distributions will be called ‘current distributions’. The implementation of the routines for the solution of Eq. (4) produces a very large amount of additional information: in relation to our choice of comparing the QT with the yearly average in the system, we decided to
Fig. 3. Distribution of the yearly average relative sensitivities.
ð9Þ
242
R. Pastres et al. / Reliability Engineering and System Safety 79 (2003) 239–244
the problem of the optimization of the interventions aimed at reaching the MPL is a complex one and, in principle, could be dealt with by estimating the costs of the remediation plan as a function of the expected reduction of each source, and by minimizing the total cost under the constraint represented by Eq. (8). However, we decided to test the method by assuming that the load of each source is decreased by the same fraction 0 , d , 1; thus setting.Dai ¼ 2dai : By substituting the above expression in Eqs. (8) and (11) one obtains: ) (2£ns X NH NO3 DIN 4 ð10Þ d < 2½QT 2 DINav ðaÞ ½Sav;i þ Sav;i ai Fig. 4. Importance of the sources, as determined by their relative squared sensitivities, averaged over the whole system.
In this way, the distribution in Fig. 3 can be read as the expected improvement, in terms of QT units, that would follow a reduction in the load of a given fraction d. For example, the complete removal of DIN from the emission of the source S5, would lead to a decrease in the concentration of DIN ranging from 0.75 to 1 QT-unit in the darkest area of Fig. 3. Furthermore, the distribution of the sensitivities shows that the effects of the reduction are less pronounced in the central and southern sub-basins: this is due to the watershed between Lido and Malamocco inlets, Fig. 1, which in our model is characterized by a low-diffusivity zone. This information, which must be regarded as an initial guess because it is based on a linear approximation, is extremely useful, because it shows immediately the area which is most affected by a given source. Sensitivity maps were used for ranking the input parameters, according to their influence on the space-time average DIN concentration, which represent our output, to be compared with QTDIN. The results are presented in Fig. 4, which shows that the loads discharged by the urban sewage system of the city of Venice and of the other islands, S11–S16, are, by far, less important than the loads from the drainage basin: this means that the improvement of Venice and its islands sewage systems would have little effect on DIN levels.
4. Discussion The results presented in Section 3 are based on a single run, which took about 25 h on a 2-CPU Digital AU533 MHz WS, and on the post-processing of the sensitivity maps, which was performed using a post-processing code: these further numerical treatments take a few seconds. In this section, we discuss two possible uses of the large amount of information stored in the sensitivity maps. 4.1. Estimation of MPL of DIN As mentioned above, Eq. (8) does not yield a unique solution in presence of two or more sources. In fact,
i¼1
MPLDIN ¼
2£ns X
ð1 2 dÞai
ð11Þ
i¼1
The substitution of the results of our simulation in Eq. (11) yields d ¼ 0:38; that is a reduction of 38% in the current load. This estimate relies on the linearization of the state equation and, therefore, in order to test this hypothesis, we performed a set of four simulations, in which the load of DIN was reduced by 20, 40, 60 and 80% of its nominal value. The results are summarised in Table 1, in which the average DIN concentration, which were estimated by using Eq. (8), are compared with the values which were actually obtained by re-running the model. The results presented in Table 1 indicates that the linearization leads to a very accurate estimation of DINav, at least for reductions within the range 0 –40%, which includes the percentual decrease in the DIN load which leads to the respect of the water QT. The spatial distributions of DIN and of phytoplankton were then computed by setting a reduction of 38% in the current loads and running the model again. The distributions thus obtained are compared with the current ones in Fig. 5. The comparison between Fig. 5(a) and (b) shows that the concentration of DIN, which is ‘measured’ in QTDIN units ðQTDIN ¼ 0:35 gN=m3 Þ; is in both cases below the threshold in a large part of the southern lagoon. However, even though the reduction computed using Eq. (6) leads to the respect of QT on the average, most of the central and northern part of the basin are still above it. As given in Section 1, DIN cannot be treated as a conservative tracer, as it interacts with the biologic communities. In particular, inorganic nitrogen is an essential Table 1 Comparison between the estimated DIN yearly average values, and the computed ones Percentual reductions in the DIN load
Estimates of DINav, based on Eq. (8) (gN/m3)
Computed DINav (gN/m3)
20 40 60 80
0.458 0.338 0.222 0.104
0.454 0.340 0.234 0.140
R. Pastres et al. / Reliability Engineering and System Safety 79 (2003) 239–244
243
Fig. 5. Current distribution of DIN (a), and phytoplankton (c), and distribution of the same variables for the input scenario which correspond to a reduction of 38% in the current load, (b) and (d).
macronutrient for algae: its abundance can lead to the eutrophication of the water body, but its lack may limit the primary production, and, therefore, the amount of biomass which is available for the other ‘consumers’, such as filter feeders (clams and mussels), and fish. The effects of a 38% reduction of the DIN load on the density of phytoplankton is illustrated by Fig. 5(c) and (d): as one can see, the highest density class, darkest grey, is not represented in Fig. 5(d).
4.2. Effects of different scenarios of nitrogen load Sensitivity maps can also be used for a quick analysis of the effects of different remediation policies which are aimed at reaching the QTs, or in relation to different choices of the space-time averages to be compared with the QTs. In fact, as we saw in Section 4.1, the respect of the threshold ‘on the average’ does not imply that the limit is respected everywhere, particularly when dealing with complex
Fig. 6. Fraction of the lagoon below QT as a function of a stepwise reduction of the load from the seven most important sources of DIN.
244
R. Pastres et al. / Reliability Engineering and System Safety 79 (2003) 239–244
transition water bodies, such as lagoons and estuaries. An alternative criterion would be that of requiring the respect of the QT in a given fraction of the waterbody: such a criterion certainly seems much more sensible when dealing with toxicants, such as heavy metals and organic micropollutants. We analysed such an alternative by estimating the effects of a stepwise reduction of the DIN emission from the seven sources S1– S7, which were found to be the most important ones. The results are summarised in Fig. 6, which represents the fraction of the lagoon below QT, left axis, and the average concentration of DIN, right, obtained by progressively reducing the load from the seven sources using a ‘step’ of 5%. This graph shows that, when the average concentration is equal to the threshold, about half of the lagoon is, above it: in order to reduce this fraction to 20%, one should plan a much more severe load reduction, equal to 50% of the current one, as it is shown in Fig. 6.
5. Conclusion Local sensitivity analysis, LSA, with respect to the parameters which characterise the sources of pollution in waterbodies appears to be a flexible and effective tool for dealing with problems which are posed by environmental policies which are based on the respect of water quality thresholds. Such approach seems particularly suitable for complex transition systems, as the simulation of their dynamics requires the use of time-consuming reaction transport models. However, it must be kept in mind that LSA relies on the linearization of the state-equation and, therefore, the results may not always be accurate. In this paper, LSA was used for obtaining a quick first estimation of the Maximum Permissible Load, which was based on the nominal simulation: this finding was then checked by taking as input the MPL which were determined and re-running the model. Such a check, which in this case-study gave satisfactory results, should always be made. The possibility of computing very quickly the distributions of the concentration of the pollutants from the results of a single run, even though with some degree of inaccuracy, is extremely interesting in cases when the remediation policies must be agreed upon in a public contexts, when the local authorities, the population and the stakeholders concerned are encouraged to express their
views. In such occasions, the post-processing module could be included in a GIS system and used interactively to show the effects of alternative remediation strategies.
Acknowledgements This work was partly supported by CORILA, Consortium for Coordination of Research Activities Concerning the Venice Lagoon System. References [1] Pastres R, Franco D, Pecenik G, Solidoro C, Dejak C. Using parallel computers in environmental modelling: a working example. Ecol Modelling 1995;80:69 –85. [2] Dejak C, Mazzei Lalatta I, Meregalli L, Pecenik G. Development of a mathematical eutrophication model of the lagoon of Venice. Ecol Modelling 1987;37:1– 20. [3] Dejak C, Mazzei Lalatta I, Messina E, Pecenik G. A two dimensional diffusion model of the lagoon of Venice and relative open boundary conditions. Ecol modelling 1987;37:21–45. [4] Dejak C, Mazzei Lalata I, Molin M, Pecenik G. Tidal threedimensional diffusion in a model of the lagoon of Venice and reliability conditions for its numerical integration. Ecol Modelling 1987;37:81–101. [5] Pastres R, Solidoro C, Cossarini G, Melaku Canu D, Dejak C. Managing the rearing of tapes philippinarum in the lagoon of Venice: a decision support system. Ecol Modelling 2001;138:231–45. [6] Regione Veneto. Piano per la prevenzione dell’inquinamento e il risanamento delle acque del bacino idrografico immediatamente sversante nella laguna di Venezia (Piano direttore 2000); 2000. [7] Collavini F, Zonta R, Bettiol C, Fagarazzi OE, Zaggia L. Metal and nutrient loads from the drainage basin to the Venice Lagoon. Preliminary Proceedings of the DRAIN Project Workshop. Venice, Italy, June 14–15; 2001. [8] Bergamasco A, Zago C. Exploring the nitrogen cycle and macroalgae dynamics in the lagoon of Venice using a multibox model. Estuarine, Coast Shelf Sci 1999;48:155–75. [9] Saltelli A, Tarantola S, Chan K. A quantitative, model independent method for global method for global sensitivity analysis of model output. Technometrics 1999;41:39–56. [10] Pastres R, Chan K, Solidoro C, Dejak C. Global sensitivity analysis of a shallow-water 3D eutrophication model. Comput Phys Commun 1999;117:62–74. [11] Koda M, Drogu AH, Seinfeld JH. Sensitivity analysis of partial differential equations with applications to reaction and diffusion processes. J Comput Phys 1979;30:259–82. [12] Pastres R, Solidoro C. Trattamento ed analisi dei dati anche mediante un modello di qualita` dell’acqua. Technical Report. Consorzio Venezia Nuova; 1999.