Mapping and Dilution Estimation of Wastewater Discharges Based on Geostatistics Using an Autonomous Underwater Vehicle

two vertical thrusters, one at the front and the other at the rear. This conﬁguration allows for small operational speeds and high maneuverability, including pure vertical motions. It is equipped with an omnidirectional acoustic transducer and an electronic system that allows for long baseline navigation. The vehicle can be programmed to follow predeﬁned trajectories while


Introduction
Wastewaters are often discharged into coastal waters through outfall diffusers that efficiently dilute effluent and usually restrict any environmental impact within a small area. However, predicting this impact is difficult because of the complexity of the hydrodynamic processes that mix the wastewater and also because of the variability in oceanic conditions. Despite great improvements over the years in the understanding of these mixing processes, since models are now available that can make reasonable predictions under steady-state conditions (Hunt et al., 2010), many aspects remain unknown and unpredictable. For this reason, much effort has been recently devoted to improve ways of monitoring and characterizing sewage plumes under a variety of oceanographic conditions.

MARES AUV
Autonomous Underwater Vehicles (AUVs) have been used efficiently in a wide range of applications. They were first developed with military applications in mind, for example for mine hunting missions. Later on, scientists realized their true potential and started to use them as mobile sensors, taking measurements in difficult scenarios and at a reasonable cost (Bellingham, 1997;Bellingham et al., 1992;Fernandes et al., 2000;Nadis, 1997;Robinson et al., 1999;Yu et al., 2002). MARES (Modular Autonomous Robot for Environment Sampling) AUV has been successfully used to monitor sea outfalls discharges 2011a;;c) (see Fig. 1). MARES is 1.5 m long, has a diameter of 8-inch and weighs about 40 kg in air. It features a plastic hull with a dry mid body (for electronics and batteries) and additional rings to accommodate sensors and actuators. Its modular structure simplifies the system's development (the case of adding sensors, for example). It is propelled by two horizontal thrusters located at the rear and two vertical thrusters, one at the front and the other at the rear. This configuration allows for small operational speeds and high maneuverability, including pure vertical motions. It is equipped with an omnidirectional acoustic transducer and an electronic system that allows for long baseline navigation. The vehicle can be programmed to follow predefined trajectories while

Mapping and Dilution Estimation of Wastewater
Discharges Based on Geostatistics Using an Autonomous Underwater Vehicle 11 collecting relevant data using the onboard sensors. A Sea-Bird Electronics 49 FastCAT CTD had already been installed onboard the MARES AUV to measure conductivity, temperature and depth. MARES' missions for environmental monitoring of wastewater discharges are conducted using a GUI software that fully automates the operational procedures of the campaign . By providing visual and audio information, this software guides the user through a series of steps which include: (1) real time data acquisition from CTD and ADCP sensors, (2) effluent plume parameter modeling using the CTD and ADCP data collected, (3) automatic path creation using the plume model parameters, (4) acoustic buoys and vehicle deployment, (5) automatic acoustic network setup and (6) real time tracking of the AUV mission.

Data processing
Data processing is the last step of a sewage outfall discharge monitoring campaign. This processing involves the ability to extrapolate from monitoring samples to unsampled locations. Although very chaotic due to turbulent diffusion, the effluent's dispersion process tends to a natural variability mode when the plume stops rising and the intensity of turbulent fluctuations approaches to zero (Hunt et al., 2010). It is likely that after this point the pollutant substances are spatially correlated. In this case, geostatistics appears to be an appropriate technique to model the spatial distribution of the effluent. In fact, geostatistics has been used with success to analyze and characterize the spatial variability of soil properties, to obtain information for assessing water and wind resources, to design sampling strategies for monitoring estuarine sediments, to study the thickness of effluent-affected sediment in the vicinity of wastewater discharges, to obtain information about the spatial distribution of sewage pollution in coastal sediments, among others. As well as giving the estimated values, geostatistics provides a measure of the accuracy of the estimate in the form of the kriging variance. This is one of the advantages of geostatistics over traditional methods of assessing pollution. In this work ordinary block kriging is used to model and map the spatial distribution of temperature and salinity measurements gathered by an AUV on a Portuguese sea outfall monitoring campaign. The aim is to distinguish the effluent plume from the receiving waters, characterize its spatial variability in the vicinity of the discharge and estimate dilution.

Geostatistical analysis 2.1 Stationary random function models
The most widely used geostatistical estimation procedures use stationary random function models. A random function is a set of random variables that have some spatial locations and whose dependence on each other is specified by some probabilistic mechanism. A random function is stationary if all the random variables have the same probability distribution and if any pair of random variables has a joint probability distribution that depends only on the separation between the two points and not on their locations. If the random function is stationary, then the expected value and the variance can be used to summarize the univariate behavior of the set of random variables. The parameter that is commonly used to summarize the bivariate behavior of a stationary random function is its covariance function, its correlogram, and its variogram. The complete definition of the probabilistic generating mechanism of a random function is usually difficult even in one dimension. Fortunately, for many of the problems we typically encounter, we do not need to know the probabilistic generating mechanism. We usually adopt a stationary random function as our model and specify only its covariance or variogram (Isaaks & Srivastava, 1989;Kitanidis, 1997;Wackernagel, 2003).

Ordinary kriging
Ordinary kriging method is often referred with the acronym BLUE which stands for "Best Linear Unbiased Estimator". "Linear" because its estimates are weighted linear combinations of the available data; "Unbiased" since it tries to have the mean error equal to 0; and "Best"because it aims at minimizing the variance of the errors. Let us then see how the concept of a random function model can be used to decide how to weight the nearby samples so that our estimates are unbiased. For any point at which we want to estimate the unknown value, our model is a stationary random function that consists of n random variables, one for the value at each of the n sample locations, Z(x 1 ), Z(x 2 ), . . . , Z(x n ), and one for the unknown value at the point we are trying to estimate Z(x 0 ). Each of these random variables has the same probability law; at all locations, the expected value of the random variable is m and the variance is σ 2 . Every value in this model is seen as an outcome (or realization) of the random variable. Our estimate is also a random variable since it is a weighted linear combination of the random variables at the n sampled locations (Cressie, 1993;Goovaerts, 1997;Isaaks & Srivastava, 1989;Kitanidis, 1997;Stein, 1999;Wackernagel, 2003;Webster & Oliver, 2007):Ẑ The estimation error is defined as the difference between the random variable modeling the true value and the estimate: The estimation error is also a random variable. Its expected value, often referred to as the bias, is Setting this expected value to 0, to ensure an unbiasedness estimate results in: This is known as the condition of unbiasedness (Isaaks & Srivastava, 1989;Kitanidis, 1997;Wackernagel, 2003). The expression of the variance of the modeled error is Since we have already assumed that all of the random variables have the same variance σ 2 , then var [Z(x 0 )] = σ 2 . The second term in Equation 5 can be written as The third term of Equation 5 can be expressed as Then, the expression of the error variance comes in the following way: Equation 8 expresses the error variance as function of the n weights, once chosen the random model function parameters, namely the variance σ 2 and all the covariances C ij . The minimization of σ 2 ε(x 0 ) is constrained by the unbiasedness condition imposed earlier, which can be solved using the method of Lagrange multipliers. We start by introducing a new parameter μ, called the Lagrange multiplier, in Equation 8 in the following way: Then we minimize σ 2 ε(x 0 ) by calculating the n + 1 partial first derivatives of Equation 9 with respect to the n weights and the Lagrange multiplier, and setting each one to 0, which produces the following system of equations: which can also be written in a compact way as This system of equations, often referred to as the ordinary kriging system, can be written in matrix notation as The set of weights and the Lagrange multiplier that will produce an unbiased estimate of Z(x 0 ) with the minimum error variance are then given by The value of σ 2 ε(x 0 ) can be obtained in a quicker way using an alternative expression to Eq. 8. Multiplying each of the n equations given in Eq. 10 by w i and summing these n equations leads to the following: Substituting this into Equation 8 the minimized error variance comes as follows: or, in terms of matrices as The minimized error variance is usually called the ordinary kriging variance.

Block kriging
A consideration in many environmental applications has been that ordinary kriging usually exhibits large prediction errors (Bivand et al., 2008). This is due to the larger variability in the observations. When predicting averages over larger areas, i.e. within blocks, much of the variability averages out and consequently block mean values have lower prediction errors. If the blocks are not too large the spatial patterns do not disappear. The block kriging system is similar to the point kriging system given by Equation 11. The matrix C is the same since it is independent of the location at which the block estimate is required. The covariances for the vector D are point-to-block covariances. Supposing that the mean value over a block V is approximated by the arithmetic average of the N point variables contained within that block (Goovaerts, 1997;Isaaks & Srivastava, 1989), i.e.
the point-to-block covariances required for vector D are The block kriging variance is where C VV is the average covariance between pairs of points within V: An equivalent procedure, that can be computationally more expensive than block kriging, is to obtain the block estimate by averaging the N kriged point estimates within the block (Goovaerts, 1997;Isaaks & Srivastava, 1989).

Spatial continuity
Spatial continuity exists in most earth science data sets. When we look at a contour map, or anything similar, the values do not appear to be randomly located, but rather, low values tend to be near other low values and high values tend to be near other high values. I.e. two measurements close to each other are most likely to have similar values than two measurements far apart (Isaaks & Srivastava, 1989). To compute the set of weights and the Lagrange multiplier, that will produce each estimate and the resulting minimized error variance, we need to know the covariances of C and D matrices. As we said before, since our random function is stationary, all pairs of random variables separated by a distance and direction h (known as lag) have the same joint probability distribution. The covariance function, C(h) is the covariance between random variables separated by a lag h (Isaaks & Srivastava, 1989;Kitanidis, 1997;Wackernagel, 2003). For a stationary random function, the covariance function C(h) is: The covariance between random variables at identical locations is the variance of the random function: The semivariogram, or simply variogram, is half the expected squared difference between random variables separated by a lag h: The quantity γ(h) is known as the semivariance at lag h. The "semi"refers to the fact that it is half of a variance. The variogram between random variables at identical locations is zero, i.e. γ(0) = 0. Using Equations 19, 20 and 21, we can relate the variogram with the covariance function as: In practice, the pattern of spatial continuity chosen for the random function is usually taken from the spatial continuity evident in the sample data set. Geostatisticians usually define the spatial continuity of the sample data set through the variogram and solve the ordinary kriging system using covariance (Isaaks & Srivastava, 1989). The maximum value reached by the variogram is called the sill. The distance at which the sill is reached is called the range. The vertical jump from zero at the origin to the value of semivariance at extremely small separation distances is called the nugget effect. The estimator of the variogram usually used, known as Matheron's method-of-moments estimator (MME) is (Matheron, 1965;Webster & Oliver, 2007) where z(x i ) is the value of the variable of interest at location x i and N(h) is the number of pairs of points separated by the particular lag vector h. Cressie and Hawkins (Cressie & Hawkins, 1980) developed an estimator of the variogram that should be robust to the presence of outliers and enhance the variogram spatial continuity, having also the advantage of not spreading the effect of outliers in computing the maps. This estimator (CRE) is defined as follows (Cressie & Hawkins, 1980): Once the sample variogram has been calculated, a function (called the variogram model) has to be fit to it. First, because the matrices C and D may need semivariance values for lags that are not available from the sample data. And second, because the use of the variogram does not guarantee the existence and uniqueness of solution to the ordinary kriging system. The most commonly used variogram models are the spherical model, the exponential model, the Gaussian model and the Matern model (Isaaks & Srivastava, 1989).

Cross-validation
Cross-validation is a procedure used to compare the performance of several competing models (Webster & Oliver, 2007). It starts by splitting the data set into two sets: a modelling set and a validation set. Then the modelling set is used for variogram modelling and kriging on the locations of the validation set. Finally the measurements of the validation set are compared to their predictions (Bivand et al., 2008). If the average of the cross-validation errors (or Mean Error, ME) is close to 0, we may say that apparently the estimates are unbiased (Z(x i ) andẐ(x i ) are, respectively, the measurement and estimate at point x i and m is the number of measurements of the validation set). A significant negative (positive) mean error can represent systematic overestimation (underestimation). The magnitude of the Root Mean Squared Error (RMSE) is particularly interesting for comparing different models (Wackernagel, 2003;Webster & Oliver, 2007): The RMSE value should be as small as possible indicating that estimates are close to measurements. The kriging standard deviation represents the error predicted by the estimation method. Dividing the cross-validation error by the corresponding kriging standard deviation allows to compare the magnitudes of both actual and predicted error (Wackernagel, 2003;Webster & Oliver, 2007). Therefore, the average of the standardized squared cross-validation errors (or Mean Standardized Squared Error, MSSE) should be about one, indicating that the model is accurate. A scatterplot of true versus predicted values provides additional evidence on how well an estimation method has performed. The coefficient of determination R 2 is a good index for summarizing how close the points on the scatterplot come to falling on the 45-degree line passing through the origin (Isaaks & Srivastava, 1989). R 2 should be close to one.

Study site
A map of the study site is shown in Fig. 2(a). Foz do Arelho outfall is located off the Portuguese west coast near Óbidos lagoon. In operation since June 2005, is presently discharging about 0.11 m 3 /s of mainly domestic wastewater from the WWTPs of Óbidos, Carregal, Caldas da Rainha, Gaeiras, Charneca and Foz do Arelho, but it can discharge up to 0.35 m 3 /s. The total length of the outfall, including the diffuser, is 2150 m. The outfall pipe, made of HDPE, has a diameter of 710 mm. The diffuser, which consists of 10 ports spaced 8 or 12 meters apart, is 93.5 m long. The ports, nominally 0.175 m in diameter, are discharging upwards at an angle of 90 • to the pipe horizontal axis; the port height is about 1 m. The outfall direction is southeast-northwest (315.5 • true bearing) and is discharging at a depth of about 31 m. In that area the coastline itself runs at about a 225 • angle with respect to true north and the isobaths are oriented parallel to the coastline. A seawater quality monitoring program for the outfall has already started in May 2006. Its main purposes are to evaluate the background seawater quality both in offshore and nearshore locations around the vicinity of the sea outfall and to follow the impacts of wastewater discharge in the area. During the campaign the discharge remained fairly constant with an average flowrate of approximately 0.11 m 3 /s. The operation area specification was based on the outputs of a plume prediction model (Hunt et al., 2010) which include mixing zone length, spreading width, maximum rise height and thickness. The model inputs are, besides the diffuser physical characteristics, the water column stratification, the current velocity and direction, and the discharge flowrate. Information on density stratification was obtained from a vertical profile of temperature and salinity acquired in the vicinity of the diffuser two weeks before the campaign (see Fig. 3). The water column was weakly stratified due to both low-temperature and salinity variations. The total difference in density over the water column was about 0.13 σ-unit. The current direction of 110 • was estimated based on predictions of wind speed and direction of the day of the campaign. A current velocity of 0.12 m/s was estimated based on historic data. The effluent flowrate consider for the plume behavior simulation was 0.11 m 3 /s. According to the predictions of the model, the plume was spreading 1 m from the surface, detached from the bottom and forming a two-layer flow. The end of the mixing zone length was predicted to be 141 m downstream from the diffuser. Fig. 2

Exploratory analysis
In order to obtain elementary knowledge about the temperature and salinity data sets, conventional statistical analysis was conducted (see the results in Table 1 and Table 2). At the depth of 1.5 m the temperature ranged from 15.359ºC to 15.562ºC and at the depth of 3 m the temperature ranged from 15.393ºC to 15.536ºC. The mean value of the data sets was 15.463ºC and 15.469ºC, respectively at the depths of 1.5 m and 3 m, which was very close to the median value that was respectively 15.466ºC and 15.472ºC. The coefficient of skewness is relatively low (-0.309) for the 1.5 m data set and not very high (-0.696) for the 3 m data set, indicating that in the first case the histogram is approximately symmetric and in the second case that distribution is only slightly asymmetric. The very low values of the coefficient of variation (0.002 and 0.001) reflect the fact that the histograms do not have a tail of high values. At the depth of 1.5 m the salinity ranged from 35.957 psu to 36.003 psu and at the depth of 3 m the salinity ranged from 35.973 psu to 36.008 psu. The mean value of the data sets was 35.991 psu and 35.996 psu, respectively at the depths of 1.5 m and 3 m, which was very close to the median value that was respectively 35.990 and 35.998 psu. The coefficient of skewness is not to much high in both data sets (-0.63 and -1.1) indicating that distributions are only slightly asymmetric. The very low values of the coefficient of variation (0.0002 and 0.0001) reflect the fact that the histograms do not have a tail of high values. The ordinary kriging method works better if the distribution of the data values is close to a normal distribution. Therefore, it is interesting to see how close the distribution of the data values comes to being normal. Fig. 4 shows the plots of the normal distribution adjusted to the histograms of the temperature measured at depths of 1.5 m and 3 m, and Fig. 5 shows the plots of the normal distribution adjusted to the histograms of the salinity measured at depths of 1.5 m and 3 m. The density value in the histogram is the ratio between the number of samples in a bin and the total number of samples divided by the width of the bin (constant). Apart from some erratic high values it can be seen that the histograms are reasonably close to the normal distribution.

Variogram modeling
For the purpose of this analysis, the temperature and the salinity measurements were divided into a modeling set (comprising 90% of the samples) and a validation set (comprising 10% of the samples). Modeling and validation sets were then compared, using Student's-t test, to check that they provided unbiased sub-sets of the original data. Furthermore, sample variograms for the modeling sets were constructed using the MME estimator and the CRE estimator. This robust estimator was chosen to deal with outliers and enhance the variogram's spatial continuity. An estimation of semivariance was carried out using a lag distance of 2 m. Table 3 and Table 4 show the parameters of the fitted models to the omnidirectional sample variograms constructed using MME and CRE estimators. All the variograms were fitted to Matern models (for several shape parameters ν) with the exception to the salinity data measured at the depth of 3 m. The range value (in meters) is an indicator of extension where autocorrelation exists. The variograms of salinity show significant differences in range. The autocorrelation distances are always larger for the CRE estimator which may demonstrate the enhancement of the variogram's spatial continuity. All variograms have low nugget values which indicates that local variations could be captured due to the high sampling rate and to the fact that the variables under study have strong spatial dependence. Anisotropy was investigated by calculating directional variograms. However, no anisotropy effect could be shown.

Cross-validation
The block kriging method was preferred since it produced smaller prediction errors and smoother maps than the point kriging. Using the 90% modeling sets of the two depths, a two-dimensional ordinary block kriging, with blocks of 10 × 10 m 2 , was applied to estimate temperature at the locations of the 10% validation sets. The validation results for both parameters measured at depths of 1.5 m and 3 m depths are shown in Table 5 and Table 6. At both depths temperature was best estimated by the variogram constructed using CRE. Salinity at the depth of 1.5 m was best estimated by the variogram constructed using CRE and at the depth of 3 m was best estimated using the Gaussian model with the MME. The  Table 5. Cross-validation results for the temperature maps at depths of 1.5 and 3 m. difference in performance between the two estimators: block kriging using the MME estimator (MBK) or block kriging using the CRE estimator (CBK) is not substantial. Fig. 6 shows the omnidirectional sample variograms for temperature at the depth of 1.5 m and 3 m fitted by the preferred models. Fig. 7 shows the omnidirectional sample variograms for salinity at the depth of 1.5 m and 3 m fitted by the preferred models. Fig. 8 and Fig. 9 show the scatterplots of true versus estimated values for the most satisfactory models. The dark line is the 45º line passing through the origin and the discontinuous line is the OLS (Ordinary Least Squares) regression line. These plots show that observed and predicted values are highly positively correlated. The R 2 value for the temperature at the depth of 1.5 m was 0.9211 and the RMSE was 0.0088248ºC, and at the depth of 3 m was 0.8827 and the RMSE was 0.0058316ºC (Table 5). The R 2 value for the salinity at the depth of 1.5 m was 0.9513 and the RMSE was 0.0016435 psu, and at the depth of 3 m was 0.8982 and the RMSE was 0.0019793 psu (    Fig. 9. Predicted versus observed salinity at the depths of 1.5 m (left) and 3 m (right) using the preferred models.

Mapping
Fig . 10 shows the block kriged maps of temperature on a 2 × 2 m 2 grid using the preferred models. Fig. 13 shows the block kriged maps of salinity on a 2 × 2 m 2 grid using the preferred models. In the 1.5 m kriged map the temperature ranges between 15.407ºC and 15.523ºC and the average value is 15.469ºC (the measured range is 15.359ºC-15.562ºC and the average value is 15.463ºC). In the 3 m kriged map the temperature ranges between 15.429ºC and 15.502ºC and the average value is 15.467ºC (the measured range is 15.393ºC-15.536ºC and the average value is 15.469ºC). We may say that estimated values are in accordance with the measurements since their distributions are similar (identical average values, medians, and quartiles). The difference in the ranges width is due to only 5.0% of the samples in the 1.5 m depth map (2.5% on each side of the distribution) and only 5.3% of the samples in the 3.0 m depth map (3.1% on the left side and 2.2% on the rigth side of the distribution). These samples should then be identified as outliers not representing the behaviour of the plume in the established area. In the 1.5 m kriged map the salinity ranges between 35.960 psu and 36.004 psu and the average value is 35.992 psu, which is in accordance with the measurements (the measured range is 35.957psu -36.003psu and the average value is 35.991 psu). In the 3 m kriged map the salinity ranges between 35.977 psu and 36.004 psu and the average value is 35.995 psu, which is in accordance with the measurements (the measured range is 35.973psu -36.008psu and the average value is 35.996 psu). As predicted by the plume prediction model, the effluent was found dispersing close to the surface. From the temperature and salinity kriged maps it is possible to distinguish the effluent plume from the background waters. It appears as a region of lower temperature and lower salinity when compared to the surrounding ocean waters at the same depth. At the depth of 1.5 m the major difference in temperature compared to the surrounding waters is about -0.116ºC while at the depth of 3 m this difference is about -0.073ºC. At the depth of 1.5 m the major difference in salinity compared to the surrounding waters is about -0.044 psu while at the depth of 3 m this difference is about -0.027 psu. It is important to note that these very small differences in temperature and salinity were detected due to the high resolution of the CTD sensor. (Washburn et al., 1992) observed temperature and salinity anomalies in the plume in the order, respectively of -0.3ºC and -0.1 psu, when compared with the surrounding waters within the same depth range. The small plume-related anomalies observed in the maps are evidence of the rapid mixing process. Due to the large differences in density between the rising effluent plume and ambient ocean waters, entrainment and mixing processes are vigorous and the properties within the plume change rapidly (Petrenko et al., 1998;Washburn et al., 1992). The effluent plume was found northeast from the diffuser beginning, spreading downstream in the direction of current. Using the navigation data, we could later estimate current velocity and direction and the values found were, respectively, 0.4 m/s and 70ºC, which is in accordance with the location of the plume.

Dilution estimation
Environmental effects are all related to concentration C of a particular contaminant X.
Defining C a as the background concentration of substance X in ambient water and C 0 as the concentration of X in the effluent discharge, the local dilution comes as follows (Fischer et al., 1979): which can be rearranged to give C = C a S−1 S + 1 S C 0 . In the case of variability of the background concentration of substance X in ambient water the local dilution is given by where C a0 is the background concentration of substance X in ambient water at the discharge depth. This expression in 29 can be arranged to give C = C a + 1 S (C 0 − C a0 ), which in simple terms means that the increment of concentration above background is reduced by the dilution factor S from the point of discharge to the point of measurement of C. Using salinity distribution at depths of 1.5 m and 3 m we estimated dilution using Equation 29 (see the contour maps in Fig. 13). We assumed C 0 = 2.3 psu, C a0 = 35.93 psu, C a = 36.008 psu at 1.5 m depth and C a = 36.006 psu at 3 m depth. The minimum dilution estimated at the depth of 1.5 m was 705 and at the depth of 3.0 m was 1164 which is in accordance with Portuguese legislation that suggests that outfalls should be designed to assure a minimum dilution of 50 when the plume reaches surface (INAG, 1998). (Since dilution increases with the plume rising we should expect that the minimum values would be greater if the plume reached surface (Hunt et al., 2010)).

Conclusion
Through geostatistical analysis of temperature and salinity obtained by an AUV at depths of 1.5 m and 3 m in an ocean outfall monitoring campaign it was possible to produce kriged maps of the sewage dispersion in the field. The spatial variability of the sampled data has been analyzed and the results indicated an approximated normal distribution of the temperature and salinity measurements, which is desirable. The Matheron's classical estimator and Cressie and Hawkins' robust estimator were then used to compute the omnidirectional variograms that were fitted to Matern models (for several shape parameters) and to a Gaussian model. The performance of each competing model was compared using a split-sample approach.
In the case of temperature, the validation results, using a two-dimensional ordinary block kriging, suggested the Matern model (ν = 0.5 − 1.5 m and ν = 0.7 − 3.0 m) with semivariance estimated by CRE. In the case of salinity, the validation results, using a two-dimensional ordinary block kriging, suggested the Matern model (ν = 0.6 − 1.5 m and ν = 0.8 − 3.0 m) with semivariance estimated by CRE, for the depth of 1.5 m, and with semivariance estimated by MME, for the depth of 3 m. The difference in performance between the two estimators was not substantial. Block kriged maps of temperature and salinity at depths of 1.5 m and 3 m show the spatial variation of these parameters in the area studied and from them it is possible to identify the effluent plume that appears as a region of lower temperature and lower salinity when compared to the surrounding waters, northeast from the diffuser beginning, spreading downstream in the direction of current. Using salinity distribution at depths of 1.5 m and 3 m we estimated dilution at those depths. The values found are in accordance with Portuguese legislation. The results presented demonstrate that geostatistical methodology can provide good estimates of the dispersion of effluent that are very valuable in assessing the environmental impact and managing sea outfalls.

Acknowledgment
This work was partially funded by the Foundation for Science and Technology (