Within the last couple of decades much effort has been put into monitoring and analysing air pollution levels in an attempt to improve both our understanding of the scientific mechanisms involved and our ability to make predictions of future levels. In this paper we use extreme value methods to produce a statistical model for the joint distribution of surface-level ozone (O3), nitric oxide (NO) and nitrogen dioxide (NO2) daily maxima, observed at a single urban location in the UK. Much recent work on the statistical analysis of extreme values has focused on methods for multivariate extremes, however, for all of the existing methods, it is unclear how to model non-stationary data. By extending the pre-processing method for the analysis of the extremes of non-stationary univariate processes, we propose a hierarchical modelling approach for non-stationary multivariate processes. This method allows prediction of the probabilities of any marginal or joint extreme events for non-stationary multivariate data. We illustrate this by predicting marginal return levels for each of the pollutants of interest and then looking at the bivariate distribution of NO and NO2, conditional on ozone achieving a given marginal return level.