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

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.


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 [45].
For matrix covariance functions, [28] 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 [25,54,55], 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 [8,54].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 [6].In the dynamic modeling framework, there are some results in studying the spatio-temporal variability of several correlated variables: [26], for example, extended univariate spatio-temporal dynamic models to multivariate dynamic spatial models.Moreover, [38] 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 [4] 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.[1] 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 [27].The maximum likelihood estimation of heterotopic spatio-temporal models with spatial LCM components and temporal dynamics was developed by [22].A GSLib [19] routine for cokriging was properly modified in [12] to incorporate the spatio-temporal LCM, previously developed using the generalized product-sum variogram model [10].Recently, in [15] 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 [16] 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 (PM 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 [38].By using a recent fitting procedure [16], 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 (PM 10 ) are obtained by using a modified GSLib program, called "COK2ST" [12].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.

Multivariate spatio-temporal random function
Let Z(u) = [Z 1 (u), . . ., Z p (u)] T , be a vector of p spatio-temporal random functions (STRF) defined on the domain represents a multivariate spatio-temporal random function (MSTRF), where s = (s 1 , . . ., s d ) are the coordinates of the spatial domain D ⊆ R d and t the coordinate of the temporal domain T ⊆ R. Afterwards, the MSTRF will be denoted with Z and its components with Z i .The p STRF Z i , i = 1, . . ., p, are the components of Z and they are associated to the spatio-temporal variables under study; these components are called coregionalized variables [29].The observations z i (u α ), i = 1, . . ., p, α = 1, . . ., N i , of the p variables Z i , at the points u α ∈ D × T, are considered as a finite realization of a MSTRF Z.

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 , ( 1 ) • the second-order moments, 1. the variance of each component Z i , 2. the cross-covariance for each pair of STRF (Z i , Z j ), i = j, u, u ∈ D × T, i, j = 1, . . ., p, i = j; 3. the cross-variogram for each pair of STRF (Z i , Z j ), i = j, u, u ∈ D × T, i, j = 1, . . ., p, i = j.
Note that for i = j, we obtain: • the covariance of the STRF Z i , called direct covariance, or simply covariance, 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.

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, . . ., p, and let {u 1 , . . ., u N } a set of N points of a spatio-temporal domain D × T; the direct and cross-covariances of the MSTRF must satisfy the following inequality: for any choice of the N points u α and for any choice of the weights λ αi .Using the matrix notation, the (p of the direct and cross-covariances of the STRF Z i (u α ) and Z j (u β ) will be admissible if they satisfy the following condition: where λ α = λ α1 , . . ., λ αp T is a (p × 1) vector of weights λ αi .
As in the univariate case, the (p × p) matrices Γ(u α − u β ) = γ ij (u α − u β ) of the direct and cross-variograms of the STRF Z i (u α ) and Z j (u β ) will be admissible if, for any choice of the N points u α , they satisfy the following condition under the constraint:

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.

Second-order stationarity
A MSTRF Z, with p components, is second-order stationary if: • for any STRF Z i , i, . . ., p, • for any pair of STRF Z i and Z j , i, j = 1, . . ., 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: where u, u + h ∈ D × T, i, j = 1, . . ., 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.

Intrinsic hypotheses
A MSTRF Z, with p components, satisfies the intrinsic hypotheses if: • for any pair of STRF Z i and Z j , i, j = 1, . . ., p, the cross-variogram exists and it depends only on the spatio-temporal separation vector h: where u, u + h ∈ D × T, i, j = 1, . . ., p.
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: u, u + h ∈ D × T, i, j = 1, . . ., p.For i = j, the direct variogram of the STRF Z i is obtained.

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: as well as it is not invariant with respect to the sign of the vector h: However, the cross-covariance is invariant with respect to the joint exchange of the variables and the sign of the vector h:

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: 2. The cross-variogram is invariant with respect to the exchange of the variables: 3. The cross-variogram is invariant with respect to the sign of the vector h: From ( 15) and ( 16) follows that the cross-variogram is completely symmetric, as it will be pointed out in the next sections.

Separability for a MSTRF
The cross-covariance C ij for a second-order stationary MSTRF Z is separable if: where a ij are the elements of a (p × p) positive definite matrix and ρ(•) is a correlation function.In this case, it results: 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: where a ij are the elements of a (p × p) positive definite matrix, ρ S (•) is a spatial correlation function and ρ T (•) is a temporal correlation function.In the literature, many statistical tests for separability have been proposed and are based on parametric models [2,32,51], likelihood ratio tests and subsampling [46] or spectral methods [23,50].

Symmetry for a MSTRF
The cross-covariance C ij of a second-order stationary MSTRF Z, with p components, is or, equivalently, if: The cross-covariance C ij of a second-order stationary MSTRF Z, with p components, is fully symmetric if: or, equivalently, Atmospheric, environmental and geophysical processes are often under the influence of prevailing air or water flows, resulting in a lack of full symmetry [18,27,52].
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 [38][39][40]50].

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, . . ., p [44], that is: where the coefficients a ij are the elements of a (p × p) positive definite matrix, and ρ(•, •) 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.
2. In the kernel convolution method [55] the cross-covariance functions is represented as follows: where the k i are square integrable kernel functions and ρ is a valid stationary correlation function.This approach assumes that all the variables Z i , i = 1, . . ., 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.
3. In the covariance convolution for stationary processes [24,43] the cross-covariance functions is represented as follows: 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.

4.
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 [1].

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 [33].
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 [29][30][31].
Let Z be a second-order stationary MSTRF with p components, Z i , i = 1, . . ., p, the direct and cross-covariances for the spatio-temporal LCM are defined as follows: where c l are covariances, called basic structures, and the non-negative coefficients b l ij satisfy the following property: b l ij = b l ji , i, j = 1, . . ., p; hence, in the LCM, it is assumed that: with i, j = 1, . . ., p, i = j.
The matrix C for the second-order stationary Z is built as follows: 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: where each basic structure g l is a variogram and the L matrices B l of the coefficients b l ij , 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 Γ of the direct and cross-variograms is: The necessary and sufficient conditions because the model defined in ( 17) and ( 20) are admissible are: 1. c l (g l ) must be covariances (variograms), 2. the matrices B l = b l ij , called coregionalization matrices, must be positive definite.
The necessary, but not sufficient conditions, for the coefficients b l ij are the following: 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 [20], Cressie-Huang models [5], Gneiting models [27] and many others [9,35,41,42,49,52].
As discussed in [10], each basic spatio-temporal structure, g l (h s , h t ), l = 1, . . ., L, can be modelled as a generalized product-sum variogram [7], that is where • γ l (h s , 0) and γ l (0, h t ) are, respectively, the marginal spatial and temporal variograms at lth scale of variability; • k l , defined as follows is the parameter of generalized product-sum variogram model and it is such that The inequality (23) represents a necessary and sufficient condition in order that each basic structure g l (h s , h t ), with l = 1, . . ., 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 [13,14].Substituting (21) in (20), the spatio-temporal LCM can be defined through two marginals: one in space and one in time, i.e.: 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, . . ., 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.

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 [16], 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 [54], 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 [53].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 [47] of the standardized variables Zi , i = 1, . . ., p, i.e., γij (h) = 0.5 Var[ Zi (u) − Zj (u + h)], i, j = 1, . . ., p, i = j.If the differences between the estimated pseudo cross-variograms γij (h) and γji (h), i, j = 1, . . ., 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 [38], based on the asymptotic joint normality of the sample spatio-temporal cross-covariances estimators.Given a set Λ of user-chosen spatio-temporal lags and the cardinality c of Λ, let Then the tests for symmetry properties can be based on the following statistics 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 ≥ 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., where b ij are the entries of a (p × p) positive definite matrix B and ρ(•) is a spatio-temporal correlation function.Hence, as in the spatial context [54], 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 [28].
• 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 [54].

Prediction and risk assessment in space-time
For prediction purposes, various cokriging algorithms can be found in the literature [3,33].
As a natural extension of spatial ordinary cokriging to the spatio-temporal context, the linear spatio-temporal predictor can be written as where u = (s, t) ∈ D × T is a point in the spatio-temporal domain, u α = (s, t) α ∈ D × T, α = 1, . . ., N, are the data points in the same domain and Λ α (u), α = 1, . . ., N, are (p × p) matrices of weights whose elements λ ij α (u) are the weights assigned to the value of the jth variable, j = 1, . . ., p, at the αth data point, to predict the ith variable, i = 1, . . ., p, at the point u ∈ D × T. The predicted spatio-temporal random vector Z(u) at u ∈ D × T, is such that each component Z i (u), i = 1, . . ., p, is obtained by using all information at the data points u α = (s, t) α ∈ D × T, α = 1, . . ., N.
The matrices of weights Λ α (u), α = 1, . . ., N, are determined by ensuring the unbiased condition for the predictor Z(u) and the efficiency condition, obtained by minimizing the error variance [29].
The new GSLib routine "COK2ST" [12] 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 [11].
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 ), . . ., I p (u, z p )] T , be a vector of p spatio-temporal indicator random functions (STIRF) defined on the domain D × T ⊆ R d+1 , with (d ≤ 3), as follows In other words, for each coregionalized variable Z i , with i = 1, . . ., 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, . . ., 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" [12] 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: which depends solely on the lag vector h = (h s , h t ), (s, s + h s ) ∈ D 2 and (t, t + h t ) ∈ T 2 .After fitting a model for γ 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" [17] can be used for prediction purposes in space and time.

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 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 2. partial heterotopy, where the sets of the sampled points are not pairwise disjoint, that is 3. isotopy, where the sets of the sampled points coincide, that is 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 [54].
The use of multivariate correlation models well-known in the literature, such as the LCM [54], the class of non-separable and asymmetric cross-covariances, proposed by [1], or the parametric family of cross-covariances, where each component is a Matérn process [28], 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.[30] and [31] 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, [56] and [57] 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; [36] proposed estimating the range parameters of a LCM using a non-linear regression method to fit the range parameters; [37] used simulated annealing to minimize a weighted sum of squares of differences between the empirical and the modelled variograms; [48] 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; [58] 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 [21].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 [19] 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 [12].This routine could be merged with the automatic procedure for fitting the spatio-temporal LCM using the product-sum variogram model, presented in [15], 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.

Case study
In the present case study, the environmental data set, with a multivariate spatio-temporal structure, involves PM 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 PM 10 survey stations, 60 and 54 atmospherical stations for monitoring Temperature and Wind Speed, respectively.As it is highlighted in Fig. 2 Note that, during the period of interest, the PM 10 threshold value fixed by National Laws for the human health protection (i.e.50 μ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 PM 10 daily averages have Spatio-temporal modeling and prediction techniques have been applied in order to assess PM 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 [16] based on the nearly simultaneous diagonalization of several sample matrix variograms, has been applied and the product-sum variogram model [7] has been fitted to the basic components; (2) spatio-temporal cokriging based on the estimated model, in order to obtain prediction maps for PM 10 pollution levels during the period 24-29 November 2009 and indicator kriging [34] in order to construct risk maps related to the probability that predicted PM 10 concentrations exceed the threshold value (50 μ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 PM 10 daily concentrations exceed some risk levels during the period 24-29 November 2009.

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 [38].
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 PM 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 [16] 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).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: 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.

Prediction maps and risk assessment
In order to obtain the prediction maps for PM 10 pollution levels for the period 24-29 November 2009, spatio-temporal cokriging has been applied, using the routine "COK2ST" [12].Risk assessment maps have been associated to the prediction maps.Spatial indicator kriging has been applied to assess the probability that predicted PM 10 daily concentrations exceed the PM 10 threshold value fixed by National Laws for the human health protection (i.e.50 μg/m 3 ), during the period of interest.Fig. 6 and Fig. 7 show contour maps of the predicted PM 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 PM 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 PM 10 daily concentrations exceed 50 μg/m 3 is high during the predicted week.It is worth noting that on Saturday and Sunday the predicted values of PM 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 PM 10 conditions.

Probability distributions for different sites
After producing predicted maps of PM 10 daily concentrations, it is also interesting to estimate the probability distribution that PM 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 PM 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 PM 10 daily concentrations exceed the threshold fixed by National Laws (50 μ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 PM 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 PM 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 PM 10 exceeding.

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 (PM 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 PM 10 , which is still considered one of the most hazardous pollutant for human health.

Figure 2 .
Figure 2. Posting map and box plots of PM 10 daily concentrations classified by typology of survey stations evident that low (high) values of Temperature and Wind Speed are associated with high (low) values of PM 10 .

Figure 3 .
Figure 3.Time series plots of PM 10 , Temperature and Wind Speed daily averages 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 PM 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.

Figure 4 .
Figure 4. Spatio-temporal variogram and cross-variogram surfaces of PM 10 , Temperature and Wind Speed daily averages

Figure 5 .Figure 7 .Figure 8 .
Figure 5. Spatio-temporal variograms and cross-variograms fitted to variogram and cross-variogram surfaces of PM 10 , Temperature and Wind Speed daily averages p} be a vector of cp 2 cross-covariances at spatio-temporal lags k = (h s , h t ) in Λ.Moreover, let 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 × T n , where D represents the spatial domain and T n = {1, . . ., n} the temporal one, and define { C ji (h s , h t ) : (h s , h t ) ∈ Λ, i, j = 1, . . ., p}.Under the assumptions given in the above paper,