Spatial data collection increasingly turns to vector valued measurements at spatial locations. An example is the observation of pollutant measurements. Typically, several different pollutants are observed at the same sampled location, referred to as a monitoring station or gauged site. Usually, interest lies in the modeling of the joint process for the levels of the different pollutants and in the prediction of pollutant levels at ungauged sites. In this case, it is important to take into account not only the spatial correlation but also the correlation among the different variables at each gauged site. Since, conceptually, there is a potentially observable measurement vector at every location in the study region, a multivariate spatial process becomes a natural modeling choice. In using a Gaussian process, the main challenge is the specification of a valid and flexible cross-covariance function. This paper proposes a rich class of covariance functions developed through the so-called linear coregionalization model [see, e.g., Wackernagel, 1998] for multivariate spatial observations. Following the ideas in the work of, for example, Royle and Berliner [1991], we can reparameterize a multivariate spatial model using suitable univariate conditional spatial processes, facilitating the computation. We provide explicit details, including the computation of the range associated with the different component processes. As an example, we fit our model to a particular day average of CO, NO, and NO2 for a set of monitoring stations in California, USA. Copyright 2003 by the American Geophysical Union.