Advances in Spatio-Temporal Modeling and Prediction for Environmental Risk Assessment

1. Introduction

Meteorological readings, hydrological parameters and many measures of air, soil and water pollution are often collected for a certain span, regularly in time, and at different survey stations of a monitoring network. Then, these observations can be viewed as realizations of a random function with a spatio-temporal variability. In this context, the arrangement of valid models for spatio-temporal prediction and environmental risk assessment is strongly required. Spatio-temporal models might be used for different goals: optimization of sampling design network, prediction at unsampled spatial locations or unsampled time points and computation of maps of predicted values, assessing the uncertainty of predicted values starting from the experimental measurements, trend detection in space and time, particularly important to cope with risks coming from concentrations of hazardous pollutants. Hence, more and more attention is given to spatio-temporal analysis in order to sort out these issues.

Spatio-temporal geostatistical techniques provide useful tools to analyze, interpret and control the complex evolution of various variables observed by environmental monitoring networks. However, in the literature there are no specialized monographs which contain a thorough presentation of multivariate methodologies available in Geostatistics, especially in a spatio-temporal context. Several authors have developed different multivariate models for analyzing the spatial and spatio-temporal behavior of environmental variables, as it is clarified in the following brief review.

In multivariate spatial analysis, direct and cross correlations for the variables under study are quantified by estimating and modeling the matrix variogram. The difficulty in modeling this matrix function, especially the off diagonal entries of the same matrix, has been first faced by using the linear coregionalization model (LCM), proposed by [Matheron, 1982].

For matrix covariance functions, [Gneiting et al., 2010] constructed a parametric family of symmetric covariance models for stationary and isotropic multivariate Gaussian spatial random fields, where both the diagonal and off diagonal entries are of the Matérn type. In the bivariate case, they provided necessary and sufficient conditions for the positive definiteness of the second-order structure, whereas for the other multivariate cases they suggested a parsimonious model which imposes restrictions on the scale and cross-covariance smoothness parameters. In the bivariate case, where the smoothness parameter is the same for both covariance functions, the Gneiting model is a simplified LCM. The Gneiting cross-covariance model also assumes that the scale parameter is the same for all the covariance functions and the cross-covariance functions. Both the LCM and Gneiting constructions for cross-covariances result in symmetric models; however, no distributional assumptions are required for using a LCM, which can easily incorporate components with compact support and multiple ranges and an unbounded variogram component.

Although models for multivariate spatial data have been extensively explored [Wackernagel, 2003, Gelfand et al., 2004, Ver hoef Barry, 1998], models for multivariate spatio-temporal data have received relatively less attention. In the literature, it is common to use classical techniques for multivariate spatial and temporal analysis [Wackernagel, 2003, De Iaco et al., 2001b]. Recently, canonical correlation analysis was combined with space-time geostatistical tools for detecting possible interactions between two groups of variables, associated with pollutants and atmospheric conditions [De Iaco, 2011]. In the dynamic modeling framework, there are some results in studying the spatio-temporal variability of several correlated variables: [Gelfand et al., 2005], for example, extended univariate spatio-temporal dynamic models to multivariate dynamic spatial models. Moreover, [Li et al., 2008] proposed a methodology to evaluate the appropriateness of several common assumptions, such as symmetry, separability and linear model of coregionalization, on multivariate covariance functions in the spatio-temporal context, while [Choi et al., 2009] proposed a spatio-temporal LCM where the multivariate spatio-temporal process was expressed as a linear combination of independent Gaussian processes in space-time with mean zero and a separable spatio-temporal covariance. [Apanasovich Genton, 2010] considered some solutions to the symmetry problem; moreover, they proposed a class of cross-covariance functions for multivariate random fields based on the work of [Gneiting, 2002]. The maximum likelihood estimation of heterotopic spatio-temporal models with spatial LCM components and temporal dynamics was developed by [Fassò Finazzi(2011)]. A GSLib [Deutsch Journel, 1998] routine for cokriging was properly modified in [De Iaco et al., 2010] to incorporate the spatio-temporal LCM, previously developed using the generalized product-sum variogram model [De Iaco et al., 2003]. Recently, in [De Iaco et al., 2012a] an automatic procedure for fitting the spatio-temporal LCM using the product-sum variogram model has been presented and some computational aspects, analytically described by a main flow-chart, have been discussed. In [De Iaco et al., 2012 b] simultaneous diagonalization of the sample matrix variograms has been used to isolate the basic components of a spatio-temporal LCM and it has been illustrated how nearly simultaneous diagonalization of the cross-variogram matrices simplifies modeling of the matrix variogram.

In the following, after an introduction of the theoretical framework of the multivariate spatio-temporal random function and its features (Section 2), a review of recent techniques for building admissible models is proposed (Section 3). Successively, the spatio-temporal LCM, its assumptions and appropriate statistical tests are presented (Section 4) and techniques for prediction and risk assessment maps are introduced (Section 5). Some critical aspects regarding sampling, modeling and computational problems are discussed (Section 6). Finally, a case study concerning particle pollution ($P{M}_{10}$) and two atmospheric variables (Temperature and Wind Speed) in the South of Apulian region (Italy), has been presented (Section 7). Before using the spatio-temporal LCM to describe the spatio-temporal multivariate correlation structure among the variables under study, its adequacy with respect to the data has been analyzed; in particular, the assumption of symmetry of the cross-covariance function, has been properly tested [Li et al., 2008]. By using a recent fitting procedure [De Iaco et al., 2012 b], based on the simultaneous diagonalization of several symmetric real-valued matrix variograms, the basic structures of the spatio-temporal LCM which describes the spatio-temporal correlation among the variables, have been easily detected. Predictions of the primary variable ($P{M}_{10}$) are obtained by using a modified GSLib program, called COK2ST [De Iaco et al., 2010]. Then, risk maps showing the probability that the particle pollution exceeds the national law limit have been associated to predition maps and the estimation of the probability distributions for two sites of interest have been produced.

2. Multivariate spatio-temporal random function

Let $$Z(u)={[Z_1(u),\dots ,Z_p(u)]}^T,$$ be a vector of $p$ spatio-temporal random functions (STRF) defined on the domain$D\times T\subseteq {ℝ}^{d+1}$, with$(d\le 3)$, then

$\{Z(u),u=(s,t)\in D\times T\subseteq {ℝ}^{d+1}\},$E1

represents a multivariate spatio-temporal random function ( MSTRF), where $s=(s_1,\dots ,s_d)$ are the coordinates of the spatial domain $D\subseteq {ℝ}^d$ and $t$ the coordinate of the temporal domain$T\subseteq ℝ$.

Afterwards, the MSTRF will be denoted with $Z$ and its components with$Z_i$. The $p$ STRF $Z_i,i=1,\dots ,p,$ are the components of $Z$ and they are associated to the spatio-temporal variables under study; these components are called coregionalized variables [Goovaerts, 1997].

The observations $z_i(u_{\alpha }),i=1,\dots ,p,\alpha =1,\dots ,N_i,$ of the $p$ variables$Z_i$, at the points$u_{\alpha }\in D\times T$, are considered as a finite realization of a MSTRF Z.

2.1. Moments of a MSTRF

Given a MSTRF Z, with $p$ components, we define, if they exist and they are finite:

• the expected value, or first-order moment of each component$Z_i$,

$E\left[Z_i(u)\right]=m_i(u),u\in D\times T,i=1,\dots ,p;$E2

• the second-order moments,

1. the variance of each component$Z_i$,

$Var\left[Z_i(u)\right]=E{\left[Z_i(u)-m_i(u)\right]}^2,u\in D\times T,i=1,\dots ,p;$E3

2. the cross-covariance for each pair of STRF $(Z_i,Z_j),$$i\ne j$,

$Cov[Z_i(u),Z_j(u\text{'})]=C_{ij}(u,u\text{'})=E\left[(Z_i(u)-m_i(u))(Z_j(u\text{'})-m_j(u\text{'}))\right],$E4

$u,u\text{'}\in D\times T,i,j=1,\dots ,p,i\ne j;$E5

3. the cross-variogram for each pair of STRF $(Z_i,Z_j),$$i\ne j$,

$2{\gamma }_{ij}(u,u\text{'})=Cov\left[(Z_i(u)-Z_i(u\text{'})),(Z_j(u)-Z_j(u\text{'}))\right],$E6

$u,u\text{'}\in D\times T,i,j=1,\dots ,p,i\ne j.$E7

Note that for$i=j$, we obtain:

• the covariance of the STRF$Z_i$, called direct covariance, or simply covariance,

$C_{ii}(u,u\text{'})=E\left[(Z_i(u)-m_i(u))(Z_i(u\text{'})-m_i(u\text{'}))\right],$E8

with

$u,u\text{'}\in D\times T;$E9

• the direct variogram of the STRF$Z_i$,

$2{\gamma }_{ii}(u,u\text{'})=Var\left[(Z_i(u)-Z_i(u\text{'})\right],u,u\text{'}\in D\times T.$E10

These moments describe the basic features of a MSTRF, such as the spatio-temporal correlation for each variable and the cross-correlation among the variables.

2.2. Admissibility conditions

In multivariate Geostatistics, admissibility conditions concern both the cross-covariances and the cross-variograms, as described in the following.

Let $Z$ be a MSTRF, with components$Z_i,i=1,\dots ,p$, and let $\{u_1,\dots ,u_N\}$ a set of $N$ points of a spatio-temporal domain$D\times T$; the direct and cross-covariances of the MSTRF must satisfy the following inequality:

${\displaystyle \sum }_{i=1}^p{\displaystyle \sum }_{j=1}^p{\displaystyle \sum }_{\alpha =1}^N{\displaystyle \sum }_{\beta =1}^N{\lambda }_{\alpha i}{\lambda }_{\beta j}{C}_{ij}(u_{\alpha }-u_{\beta })\ge 0,$E11

for any choice of the $N$ points $u_{\alpha }$ and for any choice of the weights${\lambda }_{\alpha i}$. Using the matrix notation, the $(p\times p)$ matrices $C(u_{\alpha }-u_{\beta })=\left[C_{ij}(u_{\alpha }-u_{\beta })\right]$ of the direct and cross-covariances of the STRF $Z_i(u_{\alpha })$ and $Z_j(u_{\beta })$ will be admissible if they satisfy the following condition:

$${\displaystyle \sum _{\alpha =1}^N{\displaystyle \sum _{\beta =1}^N{\overrightarrow{\lambda }}_{\alpha }^{T}C(u_{\alpha }-u_{\beta })}}{\overrightarrow{\lambda }}_{\beta }\ge 0,$$E12

where ${\overrightarrow{\lambda }}_{\alpha }={\left[{\lambda }_{\alpha 1},\dots ,{\lambda }_{\alpha p}\right]}^T$ is a $(p\times 1)$ vector of weights${\lambda }_{\alpha i}$.

As in the univariate case, the $(p\times p)$ matrices $\Gamma (u_{\alpha }-u_{\beta })=\left[{\gamma }_{ij}(u_{\alpha }-u_{\beta })\right]$ of the direct and cross-variograms of the STRF $Z_i(u_{\alpha })$ and $Z_j(u_{\beta })$ will be admissible if, for any choice of the N points$u_{\alpha }$, they satisfy the following condition

$$-{\displaystyle \sum _{\alpha =1}^N{\displaystyle \sum _{\beta =1}^N{\overrightarrow{\lambda }}_{\alpha }^T\Gamma }}(u_{\alpha }-u_{\beta }){\overrightarrow{\lambda }}_{\beta }\ge 0,$$E13

under the constraint: 3mm $${\displaystyle \sum _{\alpha =1}^N{\overrightarrow{\lambda }}_{\alpha }=0.}$$

2.3. Stationarity hypotheses

Stationarity hypotheses allow to make inference on the MSTRF. In particular, second-order stationarity and intrinsic hypotheses concern the first and second-order moments of the MSTRF.

2.3.1. Second-order stationarity

A MSTRF Z, with p components, is second-order stationary if:

• for any STRF $Z_i,i,\dots ,p,$

$E[Z_i(u)]=m_i,u\in D\times T,i=1,\dots ,p;$E14

• for any pair of STRF $Z_i$ and $Z_j,i,j=1,\dots ,p,$ the cross-covariance $C_{ij}$ depends only on the spatio-temporal separation vector $h=(h_s,h_t)$ between the points u and u + h:

$C_{ij}(h)=E[(Z_i(u+h)-m_i)(Z_j(u)-m_j)]=E[Z_i(u+h)Z_j(u)]-m_i{m}_j,$E15

where$u,u+h\in D\times T,i,j=1,\dots ,p$. For$i=j$, the direct covariance function of the STRF $Z_i$ is obtained.

There exist several physical phenomena for which neither variance, nor the covariance exist, however it is possible to assume the existence of the variogram.

2.3.2. Intrinsic hypotheses

A MSTRF Z, with $p$ components, satisfies the intrinsic hypotheses if:

• for any STRF $Z_i,i=1,\dots ,p,$

$E\left[Z_i(u+h)-Z_i(u)\right]=0,u,u+h\in D\times T,i=1,\dots ,p;$E16

• for any pair of STRF $Z_i$ and $Z_j,i,j=1,\dots ,p,$ the cross-variogram exists and it depends only on the spatio-temporal separation vector h:

$2{\gamma }_{ij}(h)=Cov[(Z_i(u+h)-Z_i(u)),(Z_j(u+h)-Z_j(u))],$E17

where

$$u,u+h\in D\times T,i,j=1,\dots ,p.$$E18

Second-order stationarity implies the existence of the intrinsic hypotheses, however the converse is not true. Intrinsic hypotheses imply that the cross-variogram can be expressed as the expected value of the product of the increments:

$${\gamma }_{ij}(h)=12E\{[Z_i(u+h)-Z_i(u)][Z_j(u+h)-Z_j(u)]\},$$E19

$u,u+h\in D\times T,i,j=1,\dots ,p.$For$i=j$, the direct variogram of the STRF $Z_i$ is obtained.

2.3.3. Properties of the cross-covariance for second-order stationary MSTRF

Given a second-order stationary MSTRF, the cross-covariance satisfies the properties listed below.

The cross-covariance is not invariant with respect to the exchange of the variables:

$C_{ij}(h)\ne C_{ji}(h),i\ne j,$E20

as well as it is not invariant with respect to the sign of the vector$h$:

$C_{ij}(-h)\ne C_{ij}(h),i\ne j.$E21

However, the cross-covariance is invariant with respect to the joint exchange of the variables and the sign of the vector h:

$C_{ij}(h)=C_{ji}(-h).$E22

(13)

2.3.4. Properties of the cross-variogram for intrinsic MSTRF

Afterwards, the main properties of the cross-variogram for intrinsic MSTRF are given.

1. The cross-variogram vanishes at the origin, that is:

${\gamma }_{ij}(0)=0.$E23

1. The cross-variogram is invariant with respect to the exchange of the variables:

${\gamma }_{ij}(h)={\gamma }_{ji}(h).$E24

1. The cross-variogram is invariant with respect to the sign of the vector h:

${\gamma }_{ij}(-h)={\gamma }_{ij}(h).$E25

From (15) and (16) follows that the cross-variogram is completely symmetric, as it will be pointed out in the next sections.

2.3.5. Separability for a MSTRF

The cross-covariance $C_{ij}$ for a second-order stationary MSTRF Z is separable if:

$C_{ij}(h)=\rho (h)a_{ij},h=(h_s,h_t)\in D\times T,i,j=1,\dots ,p,$E26

where $a_{ij}$ are the elements of a $(p\times p)$ positive definite matrix and $\rho (\cdot )$ is a correlation function. In this case, it results:

$C_{ij}(h)/C_{ij}(h\text{'})=\rho (h)/\rho (h\text{'}),h,h\text{'}\in D\times T,i,j=1,\dots ,p,$E27

hence the changes of the cross-covariance functions, with respect to the changes of the vector h, do not depend on the pair of the STRF$Z_i,Z_j$.

The cross-covariance $C_{ij}$ for a second-order stationary MSTRF Z is fully separable if:

$C_{ij}(h_s,h_t)={\rho }_S(h_s){\rho }_T(h_t)a_{ij},(h_s,h_t)\in D\times T,i,j=1,\dots ,p,$E28

where $a_{ij}$ are the elements of a $(p\times p)$ positive definite matrix, ${\rho }_S(\cdot )$is a spatial correlation function and ${\rho }_T(\cdot )$ is a temporal correlation function. In the literature, many statistical tests for separability have been proposed and are based on parametric models [Brown et al., 2000, Guo, 1998, Shitan Brockwell, 1995], likelihood ratio tests and subsampling [Mitchell et al., 2005] or spectral methods [Fuentes, 2006, Scaccia Martin, 2005].

2.3.6. Symmetry for a MSTRF

The cross-covariance $C_{ij}$ of a second-order stationary MSTRF Z, with $p$ components, is symmetric if:

$C_{ij}(h)=C_{ij}(-h),h\in D\times T,i,j=1,\dots ,p,$E29

or, equivalently, if:

$C_{ij}(h)=C_{ji}(h),h\in D\times T,i,j=1,\dots ,p.$E30

The cross-covariance $C_{ij}$ of a second-order stationary MSTRF Z, with $p$ components, is fully symmetric if:

$C_{ij}(h_s,h_t)=C_{ij}(h_s,-h_t),(h_s,h_t)\in D\times T,i,j=1,\dots ,p,$E31

or, equivalently,

$C_{ij}(h_s,h_t)=C_{ij}(-h_s,h_t),(h_s,h_t)\in D\times T,i,j=1,\dots ,p.$E32

Atmospheric, environmental and geophysical processes are often under the influence of prevailing air or water flows, resulting in a lack of full symmetry [de Luna Genton, 2005, Gneiting, 2002, Stein, 2005].

Fig. 1 summarizes the relationships between separability, symmetry, stationarity and the LCM in the general class of the cross-covariance functions of a MSTRF Z. If a cross-covariance is separable, then it is symmetric, however, in general, the converse is not true. Moreover, the hypothesis of full separability is a special case of full symmetry.

Several tests to check symmetry and separability of cross-covariance functions can be found in the literature [Scaccia Martin, 2005, Lu Zimmerman, 2005a, Lu Zimmerman, 2005b, Li et al., 2008].

Figure 1.

Relationships among different classes of spatio-temporal covariance functions

3. Techniques for building admissible models

In the following, a brief review of the most utilized techniques to construct admissible cross-covariance models is presented.

1. For the intrinsic correlation model, the matrices $C$ are described by separable cross-covariances$C_{ij}$, $i,j=1,\dots ,p$[Mardia Goodall, 1993], that is:

$C_{ij}(u_{\alpha },u_{\beta })=\rho (u_{\alpha },u_{\beta })a_{ij},$E33

where the coefficients $a_{ij}$ are the elements of a $(p\times p)$ positive definite matrix, and $\rho (\cdot ,\cdot )$ is a correlation function. However, this model is not flexible enough to handle complex relationships between processes, because the cross-covariance function between components measured at each location always has the same shape regardless of the relative displacement of the locations. As it will be discussed in the next section, the LCM is a straightforward extension of the intrinsic correlation model.

1. In the kernel convolution method [Ver hoef Barry, 1998] the cross-covariance functions is represented as follows:

$$C_{ij}(u_{\alpha },u_{\beta })=\underset{{ℝ}^{d+1}}{{\displaystyle \int }}\underset{{ℝ}^{d+1}}{{\displaystyle \int }}{k}_i(u_{\alpha }-u)k_j(u_{\beta }-u\text{'})\rho (u-u\text{'})dudu\text{'},$$E34

where the $k_i$ are square integrable kernel functions and $$\rho$$ is a valid stationary correlation function. This approach assumes that all the variables$Z_i,i=1,\dots ,p$, are generated by the same underlying process, which is very restrictive. Moreover, this model and its parameters lack interpretability and, except for some special cases, it requires Monte Carlo integration.

1. In the covariance convolution for stationary processes [Gaspari Cohn, 1999, Majumdar Gelfand, 2007] the cross-covariance functions is represented as follows:

$C_{ij}(h)=\underset{{ℝ}^{d+1}}{{\displaystyle \int }}{C}_i(h-h\text{'})C_j(h\text{'})dh\text{'},$E35

where $C_i$ are second-order stationary covariances. The motivation and interpretation of the resulting cross-dependency structure is rather unclear. Although some closed-form expressions exist, this method usually requires Monte Carlo integration.

1. Recently, an approach based on latent dimensions and existing covariance models for univariate random fields, has been proposed; the idea is to develop flexible, interpretable and computationally feasible classes of cross-covariance functions in closed form [Apanasovich Genton, 2010].

4. Linear coregionalization model

LCM is based on the hypothesis that each direct or cross-variogram (covariogram) can be represented as a linear combination of some basic models and each direct or cross-variogram (covariogram) must be built using the same basic models [Journel Huijbregts, 1981].

LCM is utilized in several applications because of its flexibility, moreover it encouraged the development of algorithms able to estimate quickly the parameters of the selected model, assuring the admissibility conditions, even in presence of several variables [Goovaerts, 1997, Goulard, 1989, Goulard Voltz, 1992].

Let Z be a second-order stationary MSTRF with p components, $Z_i,i=1,\dots ,p$, the direct and cross-covariances for the spatio-temporal LCM are defined as follows:

$C_{ij}(h)=Cov\left[Z_i(u),Z_j(u+h)\right]={\displaystyle \sum }_{l=1}^L{b}_{ij}^l{c}_l(h),i,j=1,\dots ,p,$E36

where $c_l$ are covariances, called basic structures, and the non-negative coefficients $b_{ij}^l$ satisfy the following property:

$b_{ij}^l=b_{ji}^l,i,j=1,\dots ,p;$E37

hence, in the LCM, it is assumed that:

$C_{ij}(h)=C_{ij}(-h),$E38

$C_{ij}(h)=C_{ji}(h),$E39

with

$i,j=1,\dots ,p,i\ne j.$E40

The matrix C for the second-order stationary Z is built as follows:

$C(h)={\displaystyle \sum }_{l=1}^L{B}_l{c}_l(h).$E41

Analogously, it is also possible to introduce the LCM for a MSTRF which satisfies the intrinsic hypotheses. In such a case, the direct and cross-variograms are built as follows:

${\gamma }_{ij}(h)={\displaystyle \sum }_{l=1}^L{b}_{ij}^l{g}_l(h),i,j=1,\dots ,p,$E42

where each basic structure $g_l$ is a variogram and the L matrices $$B_l$$ of the coefficients$b_{ij}^l$, corresponding to the sill values of the basic models$g_l$, are positive definite.

Then, for a MSTRF which satisfies the intrinsic hypotheses the matrix $\Gamma$ of the direct and cross-variograms is:

$$\Gamma (h)={\displaystyle \sum }_{l=1}^L{B}_l{g}_l(h).$$E43

The necessary and sufficient conditions because the model defined in (17) and (20) are admissible are:

1. $c_l(g_l)$

2. the matrices$B_l=\left[b_{ij}^l\right]$, called coregionalization matrices, must be positive definite.

The necessary, but not sufficient conditions, for the coefficients $b_{ij}^l$ are the following:

a)

$b_{ii}^l\ge 0i=1,\dots ,p,l=1,\dots ,L;$E44

b)

$\left|b_{ij}^l\right|\le \sqrt{b_{ii}^l{b}_{jj}^l}i,j=1,\dots ,p,l=1,\dots ,L.$E45

The basic structures $g_l(h)$ $=g_l(h_s,h_t)$ of the spatio-temporal LCM (20) can be modelled by using several spatio-temporal variogram models known in the literature, such as the metric model [Dimitrakopoulos Luo, 1994], Cressie-Huang models [Cressie Huang, 1999], Gneiting models [Gneiting, 2002] and many others [De Iaco et al., 2002, Kolovos et al., 2004, Ma, 2002, Ma, 2005, Porcu, 2008, Stein, 2005].

As discussed in [De Iaco et al., 2003], each basic spatio-temporal structure, $g_l(h_s,h_t)$, $l=1,\dots ,L$, can be modelled as a generalized product-sum variogram [De Iaco et al., 2001a], that is

$g_l(h_s,h_t)={\gamma }_l(h_s,0)+{\gamma }_l(0,h_t)-k_l{\gamma }_l(h_s,0){\gamma }_l(0,h_t),l=1,\dots ,L,$E46

where

${\gamma }_l(h_s,0)$and ${\gamma }_l(0,h_t)$ are, respectively, the marginal spatial and temporal variograms at $l$th scale of variability; $k_l$, defined as follows

$k_l=\frac{sill[{\gamma }_l(h_s,0)]+sill[{\gamma }_l(0,h_t)]-sill[g_l(h_s,h_t)]}{sill[{\gamma }_l(h_s,0)]sill[{\gamma }_l(0,h_t)]},$E47

is the parameter of generalized product-sum variogram model and it is such that

$$0<k_l\le \frac{1}{max\{sill[{\gamma }_l(h_s,0)];sill[{\gamma }_l(0,h_t)]\}}.$$E48

The inequality (23) represents a necessary and sufficient condition in order that each basic structure$g_l(h_s,h_t)$, with$l=1,\dots ,L$, is admissible. Recently, it was shown that strict conditional negative definiteness of both marginals is a necessary as well as a sufficient condition for the generalized product-sum (21) to be strictly conditionally negative definite [De Iaco et al., 2011a, De Iaco et al., 2011 b].

Substituting (21) in (20), the spatio-temporal LCM can be defined through two marginals: one in space and one in time, i.e.:

$\Gamma (h_s,0)={\displaystyle \sum }_{l=1}^L{B}_l{\gamma }_l(h_s,0),\Gamma (0,h_t)={\displaystyle \sum }_{l=1}^L{B}_l{\gamma }_l(0,h_t).$E49

Using the generalized product-sum variogram model it is possible:

1. to identify the different scales of variability and build the matrices $B_l,l=1,\dots ,L,$ by means of the direct and cross marginal variograms;

2. to describe the correlation structure of processes characterized by a different spatial and temporal variability.

4.1. Assumptions in the spatio-temporal LCM

Fitting a spatio-temporal LCM to the data requires the identification of the spatio-temporal basic variograms and the corresponding positive definite coregionalization matrices, however this is often a hard step to tackle. A recent approach [De Iaco et al., 2012b], based on the simultaneous diagonalization of a set of matrix variograms computed for several spatio-temporal lags, allows to determine the spatio-temporal LCM parameters in a very simple way.

In several environmental applications [Wackernagel, 2003], the cross-covariance function is not symmetric, as for example, in time series in presence of a delay effect, as well as in hydrology, for the cross-correlation between a variable and its derivative, such as water head and transmissivity [Thiebaux, 1990]. Hence, this assumption should be tested before fitting a spatio-temporal LCM.

A useful hint to verify the symmetry of the cross-covariance can be given by estimating all the pseudo cross-variograms [Myers, 1991] of the standardized variables${\tilde{Z}}_i,i=1,\dots ,p$, i.e., ${\tilde{\gamma }}_{ij}(h)=0.5Var[{\tilde{Z}}_i(u)-{\tilde{Z}}_j(u+h)],i,j=1,\dots ,p,i\ne j.$If the differences between the estimated pseudo cross-variograms ${\tilde{\gamma }}_{ij}(h)$ and ${\tilde{\gamma }}_{ji}(h),i,j=1,\dots ,p,$ are zero or close to zero, then it could be assumed that the cross-covariances are symmetric.

The appropriateness of the assumption of symmetry of a spatio-temporal LCM can be tested by using the methodology proposed by [Li et al., 2008], based on the asymptotic joint normality of the sample spatio-temporal cross-covariances estimators. Given a set $\text{Λ}$ of user-chosen spatio-temporal lags and the cardinality $c$ of$\text{Λ}$, let $G_n=\{C_{ji}(h_s,h_t):(h_s,h_t)\in \text{Λ},i,j=1,\dots ,p\}$ be a vector of $c{p}^2$ cross-covariances at spatio-temporal lags $k=(h_s,h_t)$ in$\text{Λ}$. Moreover, let ${\widehat{C}}_{ji}(h_s,h_t)$ be the estimator of $C_{ji}(h_s,h_t)$ based on the sample data in the spatio-temporal domain$D\times T_n$, where $D$ represents the spatial domain and $T_n=\{1,\dots ,n\}$ the temporal one, and define$\{{\widehat{C}}_{ji}(h_s,h_t):(h_s,h_t)\in \text{Λ},i,j=1,\dots ,p\}$. Under the assumptions given in the above paper, $|T_n{|}^{1/2}({\widehat{G}}_n-G)\stackrel{d}{\to }{N}_{c{p}^2}(0,\text{Σ})$, where $\left|T_n\right|\text{Σ}$ converges to$Cov({\widehat{G}}_n,{\widehat{G}}_n)$. Then the tests for symmetry properties can be based on the following statistics

$TS=\left|T_n\right|{(A{\widehat{G}}_n)}^T{(A\text{Σ}{A}^T)}^{-1}(A{\widehat{G}}_n)\stackrel{d}{\to }{\chi }_a^2,$E50

where a is the row rank of the matrix A, which is such that $AG=0$ under the null hypothesis.

Moreover, the choice of modeling the MSTRF Z by a spatio-temporal LCM is based on the prior assumption that the multivariate correlation structure of the variables under study is characterized by $L(L\ge 2)$ scales of spatio-temporal variability. On the other hand, if the multivariate correlation of a set of variables does not present different scales of variability ($L=1$), then the cross-covariance functions are separable, i.e.,

$C_{ij}(h)=\rho (h)b_{ij},i,j=1,\dots ,p,$E51

where $b_{ij}$ are the entries of a $(p\times p)$ positive definite matrix B and $\rho (\cdot )$ is a spatio-temporal correlation function. Hence, as in the spatial context [Wackernagel, 2003], a spatio-temporal intrinsic coregionalization model can be considered.

Obviously, this last model is just a particular case ($L=1$) of the spatio-temporal LCM defined in (17) and it is much more restrictive than the linear model of coregionalization since it requires that all the variables have the same correlation function, with possible changes in the sill values. Note that, if a cross-covariance is separable, then it is symmetric.

Remarks

• In the spatio-temporal LCM, each component is represented as a linear combination of latent, independent univariate spatio-temporal processes. However, the smoothness of any component defaults to that of the roughest latent process, and thus the standard approach does not admit individually distinct smoothness properties, unless structural zeros are imposed on the latent process coefficients [Gneiting et al., 2010].

• In most applications, the wide use of the spatio-temporal LCM is justified by practical aspects concerning the admissibility condition for the matrix variograms (covariances). Indeed, it is enough to verify the positive definiteness of the coregionalization matrices, $B_l$, at all scales of variability.

• The spatio-temporal LCM allows unbounded variogram components to be used [Wackernagel, 2003].

5. Prediction and risk assessment in space-time

For prediction purposes, various cokriging algorithms can be found in the literature [Journel Huijbregts, 1981, Chilés Delfiner, 1999]. As a natural extension of spatial ordinary cokriging to the spatio-temporal context, the linear spatio-temporal predictor can be written as

$\widehat{Z}(u)={\displaystyle \sum }_{\alpha =1}^N{\text{ Λ}}_{\alpha }(u)Z(u_{\alpha }),$E52

where $u=(s,t)\in D\times T$ is a point in the spatio-temporal domain, $u_{\alpha }={(s,t)}_{\alpha }\in D\times T,$$\alpha =1,\dots ,N,$ are the data points in the same domain and ${\text{Λ}}_{\alpha }(u),\alpha =1,\dots ,N,$ are $(p\times p)$ matrices of weights whose elements ${\lambda }_{\alpha }^{ij}(u)$ are the weights assigned to the value of the $j$th variable, $j=1,\dots ,p,$at the $\alpha$th data point, to predict the ith variable, $i=1,\dots ,p$, at the point $u\in D\times T.$

The predicted spatio-temporal random vector $\widehat{Z}(u)$ at $u\in D\times T,$ is such that each component ${\widehat{Z}}_i(u),i=1,\dots ,p,$ is obtained by using all information at the data points$u_{\alpha }={(s,t)}_{\alpha }\in D\times T,\alpha =1,\dots ,N$.

The matrices of weights ${\text{Λ}}_{\alpha }(u),\alpha =1,\dots ,N,$ are determined by ensuring the unbiased condition for the predictor $\widehat{Z}(u)$ and the efficiency condition, obtained by minimizing the error variance [Goovaerts, 1997].

The new.$GSLib$. routine COK2ST [De Iaco et al., 2010] produces multivariate predictions in space-time, for one or all the variables under study, using the spatio-temporal LCM (20) where the basic spatio-temporal variograms are modelled as generalized product-sum variograms. An application is also given in [De Iaco et al., 2005].

Similarly, for environmental risk assessment, the formalism of multivariate spatio-temporal indicator random function ($MSTIRF$) and corresponding predictor, have to be introduced. Let

$I(u,z)={[I_1(u,z_1),\dots ,I_p(u,z_p)]}^T,$E53

be a vector of p spatio-temporal indicator random functions (STIRF) defined on the domain$D\times T\subseteq {ℝ}^{d+1}$, with$(d\le 3)$, as follows

$$I_i(u,z_i)=\{\begin{array}{c}1\text{if }{Z}_i\text{ is not greater (or not less) than the threshold }{z}_i,\\ 0\text{otherwise}\end{array}$$E54

where$z={[z_1,\dots ,z_p]}^T$. Then

$\{I(u,z),u=(s,t)\in D\times T\subseteq {ℝ}^{d+1}\},$E55

represents a MSTIRF. In other words, for each coregionalized variable $Z_i,$ with $i=1,\dots ,p,$ a STIRF $I_i$ can be appropriately defined. Then the linear spatio-temporal predictor (27) can be easily written in terms of the indicator random variables$I_i,i=1,\dots ,p$. If the spatio-temporal correlation structure of a $MSTIRF$ is modelled by using the spatio-temporal LCM, based on the product-sum, the new $GSLib$ routine COK2ST [De Iaco et al., 2010] can be used to produce risk assessment maps, for one or all the variables under study. If$p=1$, the dependence of the indicator variable is characterized by the corresponding indicator variogram of I:$2{\gamma }_{{}_{ST}}(h)=Var[I(s+h_s,t+h_t)-I(s,t)]$, which depends solely on the lag vector$h=(h_s,h_t)$, $(s,s+h_s)\in D^2$and$(t,t+h_t)\in T^2$. After fitting a model for${\gamma }_{{}_{ST}}$, which must be conditionally negative definite, ordinary kriging can be applied to generate the environmental risk assessment maps. In this case, the $GSLib$ routine K2ST [De Iaco and Posa, 2012] can be used for prediction purposes in space and time.

6. Some critical aspects

Multivariate geostatistical analysis for spatio-temporal data is rather complex because of several problems concerning:

a. sampling,

b. the choice of admissible direct and cross-correlation models,

c. the definition of automatic procedures for estimation and modeling.

Sampling problems

There exist several sampling techniques for multivariate spatio-temporal data, as specified herein. Let

$U_i=\{u_{\alpha }\in D\times T,\alpha =1,\dots ,N_i\},i=1,\dots ,p,$E56

be the sets of sampled points in the spatio-temporal domain for the p variables under study. It is possible to distinguish the following situations:

1. total heterotopy, where the sets of the sampled points are pairwise disjoint, that is

$$\forall i,j=1,\dots ,p,U_i\cap U_j=\varnothing ,i\ne j;$$E57

1. partial heterotopy, where the sets of the sampled points are not pairwise disjoint, that is

$$\exists i\ne j|U_i\cap U_j\ne \varnothing ,i,j=1,\dots ,p;$$E58

1. isotopy, where the sets of the sampled points coincide, that is

$\forall i,j=1,\dots ,p,U_i\equiv U_j.$E59

A special case of partial heterotopy is the so-called undersampling: in such a case, the points where a variable, called primary or principal, has been sampled, constitute a subset of the points where the remaining variables, called auxiliary or secondary, have been observed. The secondary variables provide additional information useful to improve the prediction of the primary variable.

Modeling problems

In Geostatistics, the main modeling problems concern the choice of an admissible parametric model to be fitted to the empirical correlation function (covariance or variogram). In particular, in multivariate analysis for spatio-temporal data it is important to identify an admissible model able to describe the correlation among several variables which describe the spatio-temporal process. In this context, it is suitable to underline that

• it is not enough to select an admissible direct variogram to model a cross-variogram;

• the direct variograms are positive functions, on the other hand cross-variograms could be negative functions;

• only some necessary conditions of admissibility are known to model a cross-variogram, as the Cauchy-Schwartz inequality, however sufficient conditions cannot be easily applied [Wackernagel, 2003].

The use of multivariate correlation models well-known in the literature, such as the LCM [Wackernagel, 2003], the class of non-separable and asymmetric cross-covariances, proposed by [Apanasovich Genton, 2010], or the parametric family of cross-covariances, where each component is a Matérn process [Gneiting et al., 2010], requires the identification of several parameters, especially in presence of many variables.

Moreover, estimation and modeling the direct and cross-correlation functions could be compromised by the sampling plan.

Computational problems

For spatial multivariate data different algorithms for fitting the LCM have been implemented in software packages. [Goulard, 1989] and [Goulard Voltz, 1992] described an iterative procedure to fit a LCM using a weighted least-squares like technique: this requires first fitting the diagonal entries, i.e. the basic variogram structures must be determined first. In contrast, [Xie Myers, 1995] and [Xie et al., 1995] developed an alternative method for modeling the matrix-valued variogram by near simultaneously diagonalizing the sample variogram matrices, without assuming any model for the variogram matrix; [Kunsch et al., 1997] proposed estimating the range parameters of a LCM using a non-linear regression method to fit the range parameters; [Lark Papritz, 2003] used simulated annealing to minimize a weighted sum of squares of differences between the empirical and the modelled variograms; [Pelletier et al., 2004] modified the Goulard and Voltz algorithm to make it more general and usable for generalized least-squares or any other least-squares estimation procedure, such as ordinary least-squares; [Zhang, 2007] developed an algorithm for the maximum-likelihood estimation for the purely spatial LCM and proved that the EM algorithm gives an iterative procedure based on quasi-closed-form formulas, at least in the isotopic case. Significant contributions concerning estimation and computational aspects of a LCM can be found in [Emery, 2010]. Unfortunately, for multivariate spatio-temporal data there does not exist software packages which perform in an unified way a) the structural analysis, b) a convenient graphical representation of the covariance (variogram) models fitted to the empirical ones and c) predictions. One of the solutions could be to extend the above mentioned techniques to space-time. Indeed, some routines already implemented in the GSLib software [Deutsch Journel, 1998] or in the module gstat of R, can only be applied in a multivariate spatial context.

In the last years, a new GSLib routine, called COK2ST, can be used to make predictions in the domain under study, utilizing the spatio-temporal LCM, based on the generalized product-sum model [De Iaco et al., 2010]. This routine could be merged with the automatic procedure for fitting the spatio-temporal LCM using the product-sum variogram model, presented in [De Iaco et al., 2012 a], in order to provide a complete and helpful package for the analyst who needs to obtain predictions in a spatio-temporal multivariate context. This is certainly the first step for other developments and improvements in this field.

7. Case study

In the present case study, the environmental data set, with a multivariate spatio-temporal structure, involves $P{M}_{10}$ daily concentrations, Temperature and Wind Speed daily averages measured at some monitored stations located in the South of Apulia region (Italy), from the 1st to the 23rd of November 2009. In particular, there are 28 $P{M}_{10}$ survey stations, 60 and 54 atmospherical stations for monitoring Temperature and Wind Speed, respectively. As it is highlighted in Fig. 2(a), over the domain of interest, almost all the $P{M}_{10}$ monitoring stations are either traffic or industrial stations, depending on the area where they are located (close to heavy traffic area or to industrialized area). The remaining monitoring stations are called peripheral. In Fig. 2(b), box plots of $P{M}_{10}$ daily concentrations classified by typology of survey stations are illustrated.

Figure 2.

Posting map and box plots of $P{M}_{10}$ daily concentrations classified by typology of survey stations

Fig. 3 shows the temporal profiles of the observed values. It is evident that low (high) values of Temperature and Wind Speed are associated with high (low) values of$P{M}_{10}$.

Figure 3.

Time series plots of$P{M}_{10}$, Temperature and Wind Speed daily averages

Note that, during the period of interest, the $P{M}_{10}$ threshold value fixed by National Laws for the human health protection (i.e. 50 $\mu g/m^3$ which should not be exceeded more than 35 times per year) has been overcome the 13rd, the 14th, the 15th and the 23rd of November 2009. Indeed, as regards this last aspect, it is worth noting that the highest $P{M}_{10}$ daily averages have been registered at all kinds of monitored stations, even at stations located at peripheral areas, likely for the transport effects caused by wind. As shown in Fig. 2(b), although the maximum $P{M}_{10}$ values registered at stations close to heavy traffic area are greater than the threshold value, it is evident that these values are less than the ones measured at the other stations.

Spatio-temporal modeling and prediction techniques have been applied in order to assess $P{M}_{10}$ risk pollution over the area of interest for the period 24-29 November 2009. In particular, the following aspects have been considered:

1. estimating and modeling spatio-temporal correlation among the variables; in the fitting stage of a spatio-temporal LCM, the recent procedure [De Iaco et al., 2012b] based on the nearly simultaneous diagonalization of several sample matrix variograms, has been applied and the product-sum variogram model [De Iaco et al., 2001a] has been fitted to the basic components;

2. spatio-temporal cokriging based on the estimated model, in order to obtain prediction maps for $P{M}_{10}$ pollution levels during the period 24-29 November 2009 and indicator kriging [Journel, 1983] in order to construct risk maps related to the probability that predicted $P{M}_{10}$ concentrations exceed the threshold value (50$\mu g/m^3$) fixed by National Laws;

3. generating and comparing, for two sites of interest (one close to an industrial area and the other one close to a heavy traffic area), the probability distributions that $P{M}_{10}$ daily concentrations exceed some risk levels during the period 24-29 November 2009.

7.1. Modeling spatio-temporal LCM

Modeling the spatio-temporal correlation among the variables under study by using the spatio-temporal LCM, requires first to check the adequacy of such model. In particular, the symmetry assumption has been checked by

a. exploring the differences between the pseudo cross-variograms of the standardized variables (standardized by subtracting the mean value and dividing this difference by the standard deviation),

b. using the methodology proposed by [Li et al., 2008].

As regards point a), the largest absolute difference has been equal to 0.135 and has been observed among the differences between the two pseudo cross-variograms concerning $P{M}_{10}$ and Wind Speed standardized data. On the other hand, for the point b), three pairs composed by six stations, with consecutive daily average measurements have been selected for the test. These pairs of monitoring stations have been picked to be approximately along the $SE-NW$ direction, since the prevalent wind direction over the area under study and during the analyzed period was$SE-NW$. Moreover, the temporal lag has been selected in connection with the largest empirical cross-correlations for all variable combinations; hence $h_t=1$ has been considered for all testing lags. Three different variables and three pairs of stations generate 9 testing pairs in the symmetry test, and consequently the degree of freedom a for the symmetry test is equal to 9. The test statistic TS (25) has been equal to 0.34 with a corresponding p-value equal to 0.99. Hence, the results from both the procedures have highlighted that the spatio-temporal LCM is suitable for the data set under study.

By using the recent fitting procedure [De Iaco et al., 2012b] based on the nearly simultaneous diagonalization of several sample matrix variograms computed for a selection of spatio-temporal lags, the basic independent components and the scales of spatio-temporal variability have been simply identified. In particular, the spatio-temporal surfaces of the variables under study have been computed for 7 and 5 user-chosen spatial and temporal lags, respectively (Fig. 4).

Figure 4.

Spatio-temporal variogram and cross-variogram surfaces of$P{M}_{10}$, Temperature and Wind Speed daily averages

Then, the 35 symmetric matrices of sample direct and cross-variograms have been nearly simultaneous diagonalization in order to detect the independent basic components. In this way, 3 scales of spatial and temporal variability have been identified: 10, 18 and 31.5 km in space, and 2.5, 3.5 and 6 days in time.

Thus, the following spatio-temporal LCM has been fitted to the observed data:

$\Gamma (h_s,h_t)=B_1{g}_1(h_s,h_t)+B_2{g}_2(h_s,h_t)+B_3{g}_3(h_s,h_t),$E60

where the spatio-temporal variograms $g_l(h_s,h_t),l=1,2,3,$ are modelled as a generalized product-sum model, i.e.

$g_l(h_s,h_t)={\gamma }_l(h_s,0)+{\gamma }_l(0,h_t)-k_l{\gamma }_l(h_s,0){\gamma }_l(0,h_t).$E61

The spatial and temporal marginal basic variogram models, ${\gamma }_l(h_s,0)$and ${\gamma }_l(0,h_t),$ respectively, and the coefficients $k_l,l=1,2,3,$ previously defined in (22), are shown below:

${\gamma }_1(h_s,0)=86Exp(|\left|h_s\right||;10),{\gamma }_1(0,h_t)=165Exp(|h_t|;2.5),k_1=0.0057,$E62

${\gamma }_2(h_s,0)=0.95Exp(|\left|h_s\right||;18),{\gamma }_2(0,h_t)=3.7Exp(|h_t|;3.5),k_2=0.02418,$E63

${\gamma }_3(h_s,0)=0.29Gau(|\left|h_s\right||;31.5),{\gamma }_3(0,h_t)=0.83Exp(|h_t|;6),k_3=1.1633,$E64

where $Exp(\cdot ;a)$ and $Gau(\cdot ;a)$ denote the well known exponential and Gaussian variogram models, with practical range a [Deutsch Journel, 1998].

Finally, the matrices $B_l,l=1,2,3,$ of the spatio-temporal LCM (31), computed by the procedure described in [De Iaco et al., 2012 a], are the following:

$B_1=\left[\begin{array}{lll}0.982\hfill & -0.044\hfill & -0.024\hfill \\ -0.044\hfill & 0.018\hfill & 0.006\hfill \\ -0.024\hfill & 0.006\hfill & 0.007\hfill \\ \hfill & \hfill & \hfill \end{array}\right],B_2=\left[\begin{array}{lll}43.421\hfill & -2.632\hfill & -1.550\hfill \\ -2.632\hfill & 0.395\hfill & 0.118\hfill \\ -1.550\hfill & 0.118\hfill & 0.105\hfill \\ \hfill & \hfill & \hfill \end{array}\right],$E65

$B_3=\left[\begin{array}{lll}71.429\hfill & -10.119\hfill & -4.167\hfill \\ -10.119\hfill & 1.726\hfill & 0.414\hfill \\ -4.167\hfill & 0.414\hfill & 0.357\hfill \\ \hfill & \hfill & \hfill \end{array}\right].$E66

Fig. 5 shows the spatio-temporal variograms and cross-variograms fitted to the surfaces of$P{M}_{10}$, Temperature and Wind Speed daily averages. Then, the spatio-temporal LCM (31) has been used to produce prediction and risk assessment maps for $P{M}_{10}$ daily concentrations, as discussed hereafter.

7.2. Prediction maps and risk assessment

In order to obtain the prediction maps for $P{M}_{10}$ pollution levels for the period 24-29 November 2009, spatio-temporal cokriging has been applied, using the routine COK2ST [De Iaco et al., 2010]. Risk assessment maps have been associated to the prediction maps. Spatial indicator kriging has been applied to assess the probability that predicted $P{M}_{10}$ daily concentrations exceed the $P{M}_{10}$ threshold value fixed by National Laws for the human health protection (i.e. 50$\mu g/m^3$), during the period of interest. Fig. 6 and Fig. 7 show contour maps of the predicted $P{M}_{10}$ values and the corresponding risk maps, for the period 24-29 November 2009. The red points on the maps represent the monitoring stations.

It is important to highlight that the highest $P{M}_{10}$ values are predicted in the Eastern part of the domain of interest: this area corresponds to the boundary between Lecce and Brindisi districts which is strictly close to an industrial site, such as the thermoelectric power station Enel-Federico II, located in Cerano (Brindisi district). Moreover, in this area the probability that $P{M}_{10}$ daily concentrations exceed 50 $\mu g/m^3$ is high during the predicted week. It is worth noting that on Saturday and Sunday the predicted values of $P{M}_{10}$ daily concentrations show lower average levels than the ones estimated during the working days, when heavy traffic contributes to keep pollution concentrations high; consequently, the corresponding risk maps do not show hazardous $P{M}_{10}$ conditions.

Figure 5.

Spatio-temporal variograms and cross-variograms fitted to variogram and cross-variogram surfaces of$P{M}_{10}$, Temperature and Wind Speed daily averages

Figure 6.

Prediction maps of $P{M}_{10}$ daily concentrations and risk maps at the threshold fixed by National Laws, for the 24th, 25th and 26th of November 2009

Figure 7.

Prediction maps of $P{M}_{10}$ daily concentrations and risk maps at the threshold fixed by National Laws, for the 27th, 28th and 29th of November 2009

7.3. Probability distributions for different sites

After producing predicted maps of $P{M}_{10}$ daily concentrations, it is also interesting to estimate the probability distribution that $P{M}_{10}$ daily concentrations exceed some risk levels at sites characterized by sources of pollution.

Two different sites, one close to an industrialized area, located at Brindisi district, and the other one close to a heavy traffic area, located at Lecce district, have been considered in order to generate and compare the probability distributions that $P{M}_{10}$ daily concentrations overcome several risk levels during the period 24-29 November 2009. Fig. 8 shows that, all over the industrialized area of interest, the probability that $P{M}_{10}$ daily concentrations exceed the threshold fixed by National Laws (50$\mu g/m^3$) is very high (greater than 80%) during the period 24-27 November 2009 (working days); on the other hand, during the weekend (28-29 November 2009) such a probability decreases at 40-42%.

Note that at the site close to a heavy traffic area, the probability that $P{M}_{10}$ daily concentrations exceed the national law limit is very low during the analyzed period, except on the 25th of November; moreover, it drops off rapidly for values below the national threshold. This empirical evidence highlights that there is no critical $P{M}_{10}$ exceeding for the selected traffic site, especially during the last 3 days of the week. As it is shown, this is a very powerful tool since any action of environmental protection might be adopted in advance by taking into account the actual likelihood of dangerous $P{M}_{10}$ exceeding.

Figure 8.

Probability distributions that $P{M}_{10}$ daily concentrations overcome several risk levels during the period 24-29 November 2009

8. Conclusions

In this paper, some significant theoretical and practical aspects for multivariate geostatistical analysis have been discussed and some critical issues concerning sampling, modeling and computational aspects, which should be faced, have been pointed out. The proposed multivariate geostatistical techniques have been applied to a case study pertaining particle pollution ($P{M}_{10}$) and two atmospheric variables (Temperature and Wind Speed) in the South of Apulian region.

Further analysis regarding the integration of land use and possible sources of pollution through an appropriate geographical information system could be helpful to fully understand the dynamics of$P{M}_{10}$, which is still considered one of the most hazardous pollutant for human health.

Acknowledgement

The authors are grateful to the Editor and the reviewers, whose comments contribute to improve the present version of the paper. This research has been partially supported by the 5per1000 project (grant given by University of Salento in 2011).