Interactive comment on “ Analysing spatio-temporal patterns of the global NO 2-distribution retrieved from GOME satellite observations using a generalized additive model ”

This manuscript presents a novel approach for analysis of satellite NO2 observations which leads to interesting results concerning long-term trends, annual and weekly cycles and the effects of changing wind directions. Although I do not fully understand how a general additive model works (how it determines the functional form of each term), the results appear to be reasonable. For some of the aspects of the analysis (e.g., weekly cycle) the results of the GAM are similar to those from prior literature. However, the GAM provides interesting new information, such as the apparent propagation of the weekly signal from western to eastern Europe. I do not believe that a wind directional


Introduction
Nitrogen oxides -NO and NO 2 , often referred to as NO xbelong to the most important atmospheric pollutants.NO 2 is poisonous by inhalation (World Health Organization, 2005) and NO x plays an important role in the atmospheric ozone budget (Jacob, 1999;Seinfeld and Pandis, 1997).NO x is influencing chemical and biological processes both locally (Uno et al., 1996;Wakamatsu et al., 1998) and globally (Stohl et al., 2003;Wenig et al., 2003), and its occurrence is closely related to human activities.Tropospheric NO x is a major contributor to tropospheric ozone smog in urban areas, and even at a global scale a disproportionally high amount of the NO x originates from anthropogenic sources (Olivier et al., 1990;Seinfeld and Pandis, 1997).
With space-borne instruments, such as the Global Ozone Monitoring Experiment (GOME), global time series of NO 2 and other tropospheric trace gases are becoming increasingly available, with a considerable resolution both in time and space (European Space Agency, 1995;Burrows et al., Published by Copernicus Publications on behalf of the European Geosciences Union.1999; Bovensmann et al., 1999;Wagner et al., 2008;Leue et al., 2001;Richter and Burrows, 2002;Martin et al., 2002;Beirle et al., 2003;Beirle, 2004a,b;Boersma et al., 2004Boersma et al., , 2008Boersma et al., , 2009)).Also, with an increasing amount of observations from other space-borne sensors and high level products derived from them, such as global 3-dimensional wind fields, the joint analysis of observational data from multiple sources is becoming more and more attractive.Here, an exploratory data-driven analysis of the remote sensing data may be of particular interest, since only a little prior knowledge about local and global interactions between the different observational variables may be available for specific questions of interest.
A number of earlier studies focused on different temporal patterns in global NO 2 observations.Examples are the weekly cycle found to correlate well with anthropogenic sources (Beirle et al., 2003) and the analysis of long term trends and seasonal cycles (Richter et al., 2005;van der A et al., 2006, 2008;Stavrakou et al., 2008).All these studies used parametric models for the seasonal variation of the temporal trends and ad hoc extensions, e.g.averaging the residual over monthly windows (Beirle, 2004b).
In the present study we followed a different approach and used a generalized additive model (GAM) to analyze spatiotemporal dynamics of the observed NO 2 in a non-parametric fashion (Hastie and Tibshirani, 1986;Wood, 2006).In addition to a parametric linear trend and a discrete weekly cycle, we used non-parametric terms for annual cycles and another non-parametric term for the interaction of the observed NO 2 and the local wind fields.Our choice of this generalized additive model was motivated by two different incentives: -Firstly, we test the data for temporal trends which are more complex than the parametric formulations used so far for global studies of trends in the NO 2 distribution.
-Secondly, the GAM allows us to approach the modeling problem of the interaction of NO 2 and local wind, where no parametric relation had to be known.
The alternative approach for modeling the dependence on the wind direction, that consists of discretizing the wind directions, would decrease the angular resolution for this term.The same is true for the alternative of estimating monthly means instead of using the GAM for modeling the annual cycle.
Non-parametric approaches are in common use to model the time series of different trace gases to pursue, for example, such tasks as identifying relationships between air pollution and public health (Dominici et al., 2002;Smith et al., 1999), but also to increase the sensitivity in the monitoring of local trace gas distributions (Aldrin et al., 2005;Kim et al., 2005).We aimed at generalizing these local approaches for the analysis of global trace gas distributions.
In general, the influence of the wind component to the observed NO 2 is particularly high close to strong and con-tinuous point sources like power plants or individual cities.Such regions are characterized by strong fluctuations of the tropospheric NO 2 concentration, which can easily be visualized by forming the ratio of the standard deviation and its counterpart from robust statistics -the median of absolute deviation -of the time series of the observed tropospheric NO 2 .High values of this ratio are, for example, found close to Hong Kong and Johannesburg (see Fig. 1), indicating that few observations contributing to large part to the overall observed NO 2 , and suggesting that these events might be due to wind-related transport from the nearby sources.So, focusing on these prominent examples may allow us to understand the contribution of local transport processes to the observed NO 2 distributions.
In our study, we used data derived from the GOME instrument on board the ERS-2 satellite, which provided one of the longest global records on tropospheric NO 2 observations at the beginning of our work.(Meanwhile also SCIAMACHY on board ENVISAT and OMI on board AURA may provide time series of similar lengths, with even better spatial resolution and/or coverage.)In the following we will detail on GOME and the satellite data used (Sect.2), and the generalized additive model adapted to our task (Sect.3).We will present global maps of spatial and temporal trends from the application of the model to the GOME data, and discuss the influence of the local wind field on the NO 2 trace gas dynamics observed (Sect.4).

GOME instrument and data retrieval
The GOME instrument is one of several instruments aboard the European research satellite ERS-2 (European Space Agency, 1995;Burrows et al., 1999).It consists of a set of four spectrometers that simultaneously measure sunlight scattered and reflected from the Earth's atmosphere and ground in a total of 4096 spectral channels, covering the wavelength range between 240 and 790 nm, with moderate spectral resolutions.While GOME was primarily designed for the observation of the ozone layer, many other trace gases can also be retrieved from the GOME spectra (several of them for the first time from space, see e.g.Burrows et al., 1999).The satellite operates in a nearly polar, sunsynchronous orbit at an altitude of 780 km and crosses the equator at approximately 10:30 local time.While the satellite orbits the Earth in an almost north-south direction, the GOME instrument scans the surface of the Earth in the perpendicular east-west direction.During one scan, three individual ground pixels are observed, each covering an area of 320 km east to west by 40 km north to south.They lie side by side: a western, a center, and an eastern pixel.The Earth's surface is totally covered within 3 days (poleward from about 70 • latitude within 1 day).During this three day orbital repetition pattern, the local equator crossing time varies by about 35 min.Over the considered period ( [1996][1997][1998][1999][2000][2001], these patterns of equator crossing times stayed constant.The ERS-2 repeat cycle of 35 days hast to be kept in mind for the analysis of the weekly cycle since this period is a multiple of one week (see Sect. 4.3).
In the raw spectra from GOME, the NO 2 absorption around 430 nm was used and analyzed by differential optical absorption spectroscopy (DOAS, see Platt and Stutz (2008); Wagner et al. (2008)).Besides the NO 2 cross section, also those for O 3 , H 2 O, and the oxygen dimer O 4 , as well as a Ring spectrum were included in the analysis (for details of the spectral analysis see also Leue et al. (2001); Beirle (2004a)).Output of the DOAS analysis was the NO 2 slant column density (SCD), the NO 2 concentration integrated along the atmospheric absorption path.
The vertical column density (VCD), the amount of molecules in a vertical column, was calculated from the SCD by means of a modeled Air Mass Factor (AMF) (Solomon et al., 1987;Leue et al., 2001).For simplicity, we used an AMF for a purely stratospheric NO 2 concentration profile (see Sect. 2.1 and also discussion in Leue et al. (2001); Velders et al. (2001)).To obtain the tropospheric NO 2 VCD from the total VCD, the stratospheric part of the total VCD had to be subtracted.In this study the stratospheric NO 2 VCD was estimated over a reference sector over the Pacific Ocean and then subtracted from the total NO 2 VCD measured at any location at the same latitude (see also Richter and Burrows (2002)).The difference was then used as the estimate of the tropospheric NO 2 vertical column density, referred to as NO 2 TVCD in the following.The data used in this study were from the period beginning on 17 January 1996 and ending on 31 December 2001.We re-sampled the data at a spatial resolution of 0.5×0.5 degrees of latitude and longitude, respectively, with 0.5 • being the approximate north-south range of the scan.This may be a reasonable compromise between obtaining high resolution maps and working with a resolution which can be supported by the satellite data.

Radiative transfer effects on the retrieval of tropospheric NO 2 VCD
In contrast to higher layers of the atmosphere, the radiative transfer in the troposphere, especially in the boundary layer, is rather complex.The increased air density, and also the higher probability for the presence of aerosols and clouds, cause many photons to be scattered more than once.In addition, the reflection and absorption at the surface becomes more important.Thus, knowledge on these parameters is required for the accurate determination of the tropospheric NO 2 VCD from the satellite observations.In our study, we did not account for all of these effects in detail (instead we used a stratospheric AMF), because of several reasons.Firstly, the detailed consideration of these parameters requires a huge effort, which we wanted to avoid.Secondly, accurate information on these parameters (especially on the vertical profiles of clouds, aerosols and NO 2 ) is usually not available.Even if external information (e.g. from chemical models) is used, systematic biases remain.Thirdly, most of our results describe the relative variation of the observed tropospheric NO 2 VCD.The relative variation of most of the model components of our GAM -like weekly cycle, linear trend, or wind direction (see below) -is not affected by this simplification.
We may expect, however, a systematic seasonal variation of the parameters influencing the radiative transfer in the troposphere.As a consequence, the retrieved seasonal cycle of the tropospheric NO 2 VCD may not only contain the variation of the tropospheric NO 2 concentrations, but also variations originating from parameters like the boundary layer height, or the surface albedo.
In order to give the reader an indication on the effect on the retrieved tropospheric NO 2 VCD, we calculated correction factors for different assumed profile heights and solar zenith angles (see Table 1).These correction factors are the ratio of the stratospheric AMF (which was used in our study) to the tropospheric AMF, for the different scenarios.The correction factors are between about 2 and 5.They describe the underestimation of the true values by the retrieved NO 2 VCD in this study.The smallest underestimation occurs for high layer altitude and low SZA.
In contrast to the phenomena taking place at short time scales (like weekly cycle, and wind influence) or long time scales (like linear trends), many atmospheric parameters change systematically with the season.Typically, the boundary layer height is larger in summer, thus the underestimation is smaller.In addition, especially in the mid and high latitudes, the solar zenith angle is smaller in summer, leading again to a smaller underestimation.In contrast, the surface albedo in winter is more frequently affected by snow, leading to a smaller underestimation.Finally, the situation is complicated by the seasonal cycles of clouds and aerosols.To reduce this bias, we confined ourselves to use observations with an effective cloud fraction of less than 0.3.This resulted in time series with an average length of about 400 data points per 0.5 • ×0.5 • grid box.Regarding the fact that GOME's scannings reach global coverage after three days, an unfiltered time series of six years would have the approximate length of 730 data points at mid-latitudes.Hence about half of the data were filtered out.The information about the cloud fractions was also obtained from GOME observations using the HICRU algorithm (Grzegorski et al., 2006).While part of the mentioned effects will cancel each other out, the effects of changing layer height and SZA, especially, will lead to a stronger underestimation of the true tropospheric NO 2 VCD in winter.
It should be noted that also due to the simple stratospheric corrections, components of the annual cycle of the stratosphere might be artificially transferred to the estimated tropospheric NO 2 TVCD, especially in the presence of longitudinal gradients of the stratospheric NO 2 distributions.

Wind data from ECMWF
In addition to the GOME NO 2 measurements, freely available wind data of the European Center for Medium range Weather Forecasting (ECMWF) were used (Kållberg, 2004).The resolution of the wind data is 2.5×2.5 degree in latitude and longitude.They are temporally sampled every 6 h (00:00, 06:00, 12:00, 18:00 UTC).Since the satellite crosses the equator at approximately 10:30 local time each orbit, depending on the longitude and latitude, the times of the satellite observations and of the modeled data can deviate more or less.In our eyes, in order to reduce the systematic differences caused by differences in the matches of the sampling periods of GOME data and wind data, applying the wind speeds averaged over the last 24 h is a good compromise with the temporal resolution.Due to the rather low lifetime of tropospheric NO 2 , this selection should be well representative for the transport processes to the measurement location.We will leave the question of in what detail -in spatial and temporal resolution -the wind information should be used to study the interaction of wind and trace gas open for further studies.

The Generalized Additive Model
The Generalized Additive Model (GAM) provides a general statistical framework to model the interaction between a specific feature of interest Y and a set of q (potentially) explanatory variables X = X 1 , . . ., X q .The methodology behind the GAM follows a data-driven, non-parametric approach and has greater flexibility than traditional parametric modeling.The observable Y is modeled as a superposition of separate functions f j on the features X i .Few restrictions apply to these, and either (linear) parametric models, or a non-parametric smoothing function may be chosen for f j .Certain constraints may apply to f j , such as cyclic boundary conditions.This approach has two desirable features in the exploration of large, unstructured data sets.Firstly, a GAM is able to identify those variables X i in X that are relevant to Y , even in a large set of potential candidates.Secondly, a GAM does not require the structural relationship between Y and X to be defined right from the outset when using nonparametric model terms, but is able to "learn" it from the data, individually for each variable X i .Hence their use might be indicated when no prior knowledge about the relationship is available, and one would like the data to "suggest" the appropriate functional form; or when the functional form is expected to be complex -with threshold effects or other nonlinearities -and cannot easily be represented in a parametric model.
In the analysis of the global distribution of NO 2 (our modeled observable Y ), we were interested in characterizing the structural relationship between NO 2 and the signal of selected relevant features X i -representing temporal cycles of predefined length, or the aforementioned wind fields.

Learning the structural relationship
Assume we have for a given location a time series of m measurements of Y available -in our case a set of approximately 300-700 observations of NO 2 taken in the time from 1996 to 2001 -which can be represented by a vector y of length m, i.e. with y ∈ R m .The number of available observations vary mainly because of two reasons: The first reason is the cloud cover filtering, individual for each pixel.The second is the spatially varying observational coverage, with the higher latitudes being more frequently visited by the satellite than equatorial regions.
We assume that the measurement y arises from a true value η ∈ R m , superposed by a measurement error ε ∈ R m : y=η+ε. (1) The true value η can be, for example, the time course: η = f ann (t), where t is the time of the year.In this case, X stands for the seasonal time t.We chose to estimate the functional form of f ann from the data using a univariate spline model which can be fitted to the m data pairs of (X, Y ) under the assumption of a sufficiently "smooth" behavior of f ann along time: (2) Spline basis {b k } and spline coefficients {β k } have to be either predefined (type and number of basis functions b k ), or estimated from the data (spline coefficient {β k }).For the annual signal f ann , trigonometric basis functions would be appropriate basis functions, since they fulfill periodic boundary conditions.The spline coefficients for the function f ann are found as follows: We look for the spline f ann minimizing with f i being the second derivative of f i and λ the smoothing parameter governing β k in (2) for a specific spline representation b.The latter formula (3) can be presented in matrix notation: )) building the model matrix and L building the Cholesky factor of the covariance matrix The parameters β minimizing this expression can be found in several ways, e.g.
using QR-decomposition of the system matrix O √ λL , or using iterative approaches.
The smoothing parameter λ allows us to trade the fit-error from overfitting the noise in the available observations (first term in 3), which results in a rough function f ann , with the model-error f ann from an unrealistic, overly strongly smoothed spline (second term in 3).It is adjusted to minimize the sum of both fit-and model-error in (3) -the expected prediction error -using generalized cross-validation (GCV) (Craven and Wahba, 1979).Next to the GCV method, there exist several more methods for the selection of the smoothing parameters.Methods are listed, for example, in the comparison study of Lee (2003).
The minimization problem described with (3) and ( 4) is optimal under the assumption that the residuals ε are independent realizations of a normally distributed random variable.For simplicity we assume this to be true for NO 2 .Nevertheless, the GAM potentially allows us to use different distributional models for the model residuals -for example a binomial distribution in a binary detection task, or a Poisson distribution when measuring rare events (in accordance with McCullagh et al. (1989) and Wood (2006)).
We can expand our model ( 1) and ( 2) by assuming the true value η to be a superposition of several processes, represented by a set of functional terms The functions f j may depend on one variable X j only, with f j = f j (X j ) as a univariate function, or may depend on several variables, with f j = f j (X i , X k , X l , . ..) 1≤i,k,l≤q; 1≤j ≤p .
In the present study we confined ourselves to univariate terms, and obtain a model for X and Y with p = q.y = f 1 (X 1 ) + . . .+ f p (X p ) + ε (6) www.atmos-chem-phys.net/9/6459/2009/Atmos.Chem.Phys., 9, 6459-6477, 2009 In the implementation of the "mgcv" package (Wood, 2006) used in our work (and available from cran.r-project.org),fitting Eq. ( 5) is performed by optimizing Eq. ( 3) individually for each model term f j , and iterating the procedure over all functional models f j until convergence -a procedure referred to as backfitting (Hastie and Tibshirani, 1986).The individual model terms f j are assumed to be additive.Furthermore, while the sum of the model terms, resulting from the backfitting algorithm, is unique, the model terms themselves are not, since dependence among the covariates X j can lead to more than one representation for the same sum (Hastie and Tibshirani, 1986): a seasonal model component f ann , modeling dependency on the X ann = day in the year, for example, would be assumed to be unrelated from the effects a weekly trend f week with X week = day in the week would have on η, as would be the wind component f wind with X wind = wind direction.We note that there are little restrictions on the type of model components to be used for the terms in Eq. ( 6) -we could, for example, easily introduce a rain component f rain that has X rain = amount of rainf all during last 24 h as an argument without any specific normalization of X rain to match the other arguments X i .

Model definition and identification of relevant functional terms
Functions f i of the structural relationship between Y and X i in (6) can easily be estimated from the data in a nearly automated fashion: Optimizing Eq. (3) over λ i according to the cross-validated prediction error allows us to generate a large amount of hypotheses on the functional form of f i , and to test them accordingly (Wood, 2006).The more fundamental definition of the additive model, however, i.e. the identification of the relevant predictors X i and the specification of the additive terms f i in Eq. ( 6), will require some amount of user-interaction: Although a concise model (Eq.6) with few explanatory features X i is preferred over a model with too many predictors, it is possible to start with a model incorporating all available sources of information, and the maximal amount of available observables X i .After fitting such a model, a visual inspection of the functional model terms f i will allow us to identify irrelevant parameters X i and a successive optimization of the model (compare Fig. 2): A functional term f i , which models an irrelevant variable X i , can be removed and the set of potential explanatory features X can be reduced successively.An quantitative approach in such a recursive elimination of irrelevant terms in Eq. ( 6) is provided by an analysis of variance (ANOVA).This procedure compares predictions of a model including a specific term f i , and predictions of a model without this term, for example also using predictions in a generalized cross-validation.
A subsequent statistical test on the significance of the difference between the two distributions allows us to score the importance of the tested model term by a p-value (Toutenburg and Heumann, 2008).If predictions of the complete model are significantly better than the reduced model -measured, for example, by a paired parametric t-test (Toutenburg and Heumann, 2008), or a non-parametric Cox-Wilcoxon test (Toutenburg and Heumann, 2008) -the tested model term will be retained.If the reduced model outperforms the complete model, or does not differ significantly from the full model, the model term can be dropped.Providing a quantitative score, the outcome of the test can also be used to compare the relevance of a specific term f i for the observations Y at different locations, and, hence, to map the local relevance of the different additive terms f i of the functional model (Eq.6).

Application to the spatio-temporal distributions of NO 2
The model used in the following consists of four terms f i : An additive constant µ, the linear trend f lin = s•t (thus, from a purist point of view, one might refer to our realization of Eq. ( 5) also as a mixed model), the annual cycle f ann , the weekly cycle f week , and a component f wind modeling the dependence of the NO 2 TVCD y = η + ε (see Eq. 1) from the wind direction θ: f wind (θ (t)) While the focus of our study is on the consideration of the wind direction as a new influencing parameter, the other influencing parameters -linear trend, seasonal and weekly cycle -have to be included in the study to ensure convergence.Their influence on the TVCD was worked out in several studies on these parameters (cited later in the subsections for the discussion of the results of the respective model terms), which partly include more complex data retrieval schemes and/or observations from sensors with higher spatial resolution.
The annual cycle f ann is modeled using smoothing splines with periodic boundary conditions: The weekly cycle f week is a discrete function of the day in the week.Although each of the first three terms is a function of the time, we can assume independence between the explaining variables as we expect different processes to be responsible for changes at the time scales of f lin , f week and f ann .The term f wind is a cyclic spline over the wind direction θ with the additional border conditions: Both information about wind direction and wind speed were available from the ECMWF wind data.For reasons of simplicity we confined ourselves to consider the direction of the wind at the surface only.Of course, the wind speed and its vertical variation can also have an influence on the observed NO 2 patterns, but using the wind direction at the surface alone already has a strong influence on the observed NO 2 fields, especially close to strong sources near the surface (Sect.4.4).In future studies using data sets of trace gas concentrations and winds with higher spatial resolution and better coverage, more properties of the wind field might be included.
It should also be noted that considering the wind speed in addition to the wind direction is not a trivial task, since the components modeled by GAM have a given (additive) form.One solution would be to apply a model with splines in two variables (direction and magnitude of the wind).We leave this as an extension for further studies.

Results and discussion
All additive terms in Eq. ( 7) can be visualized.Inspecting the regression coefficients for each of the terms allows to gain insight into the functional relationship modeled -for example between observed NO 2 and time of the year -and also to check the regression results for plausibility.In particular one may compare the non-parametric regression function obtained with the shape of the corresponding parametric models to see whether the latter would have been an appropriate alternative.In Fig. 2, we show the different additive functional relationships f i for four exemplary locations.Here, the first row shows the model coefficients of the different f i for a location close to Beijing (37.75 • N, 114.75 • E), representing an example with both a strong annual cycle and a strong linear trend.The second row provides an example for south Australia (22.25 • S, 125.25 • E), where the seasonal cycle is the only significant component.In the third row results for a location in the North Sea, (54.25 • N, 7.75 • E) are given, with both a dependency to f ann and f wind .In the bottom row the result for a location close to Milano (44.75 • N, 11.25 • E) provides an example for a significant linear trend, weekly and seasonal cycle, but a non-significant wind component (here even with a constant function f wind ).
Unfortunately, an inspection of all functional terms for all pixels in a global map may be impossible (or at least impractical).So, we chose specific features of the model term highlighting properties of interest in the different f i and we mapped these scalar features instead of the full set of regression coefficients: Fig. 3 shows the logarithm of the p-values of the different model terms, indicating the significance for the different components.Amplitudes of the different f i are shown in Fig. 4 as an indicator of sign and magnitude of the respective functional term.Figures 5 to 8, map individual coefficients of the weekly cycle, seasonal cycle and the wind term.For the wind term we also map the local wind directions correlating with the maximum tropospheric NO 2 pollution, both locally for the Johannesburg example (Fig. 9), and at a regional level for South Africa and Europe (Fig. 10).
In the following we will discuss observations arising from these maps in detail, focusing on understanding the different components of the additive model (Eq.7).It should be kept in mind, however, that the purpose of these maps is only to highlight specific features such as "significance", "amplitude" or "extremum" of the model terms, and that patterns identified in the maps are to be checked with the actual functional form of the different f i , i.e., by inspecting regression functions of the different terms as shown in Fig. 2.

Linear trend
In Fig. 4, we show the spatial distribution of those areas where the linear trends contributed significantly (with pvalues less than 0.001) to the NO 2 time series for the period 1996-2001.Significant linear trends can be found in areas with dominant anthropogenic sources.We also observe low p-values for linear trends with almost constant behavior -i.e.negligible slope -such as in wide areas of the Indian and Pacific Ocean.This is supported by the fact that these areas show low seasonal significance (Fig. 3).
Significant positive trends appear in China (up to approx.15% year −1 ), in the western USA (up to approx.11% year −1 ) and in the Middle East (up to approx.11% year −1 ).Negative slopes occur less often and with lower statistic significance.Examples are some European re-gions (up to approx.−12% year −1 , see also Fig. 2) and the east of USA (up to approx.−8% year −1 ).
The linear trends are in reasonable qualitative agreement with those of other studies (Richter et al., 2005;van der A et al., 2008;Stavrakou et al., 2008) in the sense that the trends show the same signs.However, the two main areas detected in van der A et al. (2008); Richter et al. (2005) that have a strong negative trend (Europe, east coast of USA) do not show a trend of high significance in our study (Fig. 3).Conversely, we detected other regions with a significant trend, such as the west coast of the USA and the western Middle East (the north of the Arabian peninsula, mainly Iraq), which were not reported so far.Explanations for the differences in the results can be the different time periods studied.

Annual cycle
The p-values for the annual cycle (Fig. 3) vary in space on continental scales mainly.This indicates that climatic conditions are expected to have the major impact instead of human sources.
In Fig. 5, the months with the maximum NO 2 TVCD are shown.According to Jaeglé et al. (2004); van der A et al. ( 2008), the major climatic influences are biomass burning (central Africa, large parts in Brazil), lightning activity (equatorial Africa) and soil emissions (Australia, Canada and the area ranging from the Sahara over the Arabian peninsula into the south of Asia till anthropogenic sources start to dominate in eastern China).In the respective regions, we observe maxima at months that agree quite well with the months reported in Jaeglé et al. (2004); van der A et al. (2008).
In regions with high anthropogenic emissions, the highest values are typically observed in winter time, indicating strong emissions (due to heating) and long atmospheric lifetimes.Examples are the US east coast, central Europe, east Asia and the track along the Trans-Siberian Railway connecting the biggest cities in Siberia.
An interesting example for the relation to wind speed and direction is the observed high significance of the annual model component over the central Indian Ocean (Fig. 3).Here, the Monsoon causes a seasonal reversal of winds (Kunhikrishnan et al. ( 2004)) and we find that the total maximum of the seasonal component in these regions is during the Monsoon transition period of September/October (Fig. 5a).
For most of the central areas of the Indian Ocean (latitudes between 20 • and 30 • S), we observe mostly two local maxima in the annual signal, again in good agreement with Kunhikrishnan et al. (2004).In this example, an apparent wind pattern ("Monsoon") has been absorbed into the annual cycle.
As discussed in Sect.2.1, the seasonal variation of the tropospheric NO 2 VCD derived in this study probably not only contains the signal of the tropospheric NO 2 concentrations, but also those of several other parameters, such as the height of the NO 2 layer or longitudinal gradients of the stratospheric Atmos.Chem.Phys., 9, 6459-6477, 2009 www.atmos-chem-phys.net/9/6459/2009/NO 2 distribution.In addition, usually the wind direction (and speed) depends systematically on season, and part of the seasonal signal might be also attributed to the influence of the wind direction as visible from the Monsoon example.In the same way, the significant annual cycle in the source free regions over the north Atlantic Ocean may be both attributed to transport processes, or to secondary effects from the preprocessing mentioned above.The latter is in accordance to results from Richter and Burrows (2002), who observed a large variability over the Atlantic Ocean at latitudes of 55 • -60 • N mainly in winter and spring.Their explanation for this artifact (that the reference sector method fails at high latitudes) could also explain the significant annual cycles in polar and sub-polar regions in Fig. 3. Also, in some highly polluted regions, where highest values are found in other than in winter months (see also Fig. 2, North Sea), this might be due to the preprocessing and/or due to the influence of the wind speed and direction.Here, we find (Fig. 10) this region to be strongly influenced by transport processes due to wind (Sect.4.4).Thus the results derived for the seasonal cycle should be treated with caution, especially when compared to the results from other studies (van der A et al., 2008;Jaeglé et al., 2004).
However, in some regions, systematic differences were found, which can be explained not by the differences in the data observation, but by the different applied models.So, in our studies, the region in the western part of the USA showing an annual signal with a maximum in summer is significantly more extended to the east than in the study of van der A et al. (2008).In these regions (in the middle of the USA), we found mainly annual cycles with two or more local maxima, and an unimodal sinusoidal model of the annual cyclewith the one annual maximum reported in (van der A et al., 2008) -may have been a too coarse approximation.Here, the seasonal cycle exhibits a good example for using a nonparametric formulation of the seasonal term in Eq. ( 7), because the functional relationship of the seasonal cycle might deviate strongly from simple parametric forms -the standard sinusoid (Fig. 2).

Weekly cycle
The significance of the weekly cycle changes on much smaller spatial scales than, for example, the annual cycle (Fig. 3).We find a highly significant contribution of this term in densely settled and highly industrialized regions.A closer look confirms that many significant points in the global map coincide with large cities.We also note a light swath pattern, in particular in regions of low latitudes.Also in Beirle et al. (2003), stripe-like structures parallel to the satellite tracks were recognized in maps concerning the weekly cycle.According to Beirle et al. (2003), the swath pattern arises due to the 35 day periodicity of the ERS-2 flight tracks, introducing a systematic weekly pattern in the viewing geometry.
In Fig. 2 (river Po) the minimum of the weekly cycle is a prominent feature.As expected, we find here a minimum on Sunday.In Fig. 6, we summarize the f week term by mapping the weekday with maximal and minimal NO 2 TVCD.We can confirm a number of earlier findings: In the USA, the minima occur mostly on Sunday, when traffic and industrial activity is reduced.The same is true for most regions in Europe and Japan.In cities in the Middle East, the minimum day is on Friday.In Israel, a significant minimum is on Saturday.No weekend effect could be detected in the large anthropogenic sources of China and South Africa.These findings agree very well with the results of former publications (Beirle et al., 2003;Boersma et al., 2009).
It is interesting to take a closer look at the day of the minimum NO 2 TVCD in the plume of large anthropogenic sources with a strong weekly cycle, as in Europe (see Fig. 8).Here, the original Sunday minimum in western Europe (5 • E-15 • E) moves eastwards with the dominating west winds: The minimum in east Poland (20 • E-23 • E) is observed on Monday, in the Ukraine (23 • E-28 • E) on Tuesday, and in western Russia (30 • E-38 • E) on Wednesday.Thus, the signal travels about 24 • , or 1700 km (at 50 • N), within 3 days, corresponding to a day-of-minimum velocity of 6.6 m/s.
Given the artificial swath pattern, caused by the 35 days repeat cycle, one may argue that these observations may not be attributed to a west-east translation of the weekly NO 2 minimum, but result from the temporal sampling of the satellite passes.This, however, might be ruled out for several reasons: Firstly, the artificial swath pattern -visible over the oceans in the tropical regions, mostly in regions which are free of anthropogenic NO 2 -sources -cannot be found for the northern Atlantic Ocean (at the latitudes of the region affected by the studied plume).Secondly, the weekly signal for areas shown in Fig. 8 is much stronger than the weekly signal along the equatorial sampling artifact.Thirdly, we do not find a systematic pattern -as in Fig. 8 -when mapping the frequency of the passes for eastern Europe.Furthermore, stripes of the width, observed in the plume, could not be found in the swath pattern.
From a more conceptual point of view one should note the very difference of 1. Tracking single "negative" plumes from Sunday to Thursday -similar to tracking plumes of unique atmospheric events, for example from volcanic outbursts, as reported for example in Schneider et al. (1999).This may be in fact impossible given the temporal and spatial coverage of the given observations.
2. The patterns visible from model terms in the GAM.These are estimated from time series of lengths of years, here averaging over multiple negative plumes.
So, Fig. 8 represents the average pattern of several hundred plumes over the whole time of observations.Such a systematic drift of the NO 2 plume has -to the best of our knowledge -not been reported before.The detection of a significant weekly cycle with a shift of the day of minimum up to 3 days is only possible for a quite high NO 2 lifetime of the order of a day or more, which however could be reached in wintertime.The analysis of the downwind evolution of the weekly pattern thus holds information on the signal travel velocity, which generally provides lifetime information if combined with the NO 2 loss along the track.However, quantitative studies of such kind should be performed with new datasets with improved spatial resolution, data preprocessing, and separately for different seasons.
When applying this model, we obtain information not only about the temporal behavior of sources, but also about the area influenced by a source that shows a characteristic temporal behavior.Such results might be important input for model simulations of the source strengths.As illustrated for the west-European NO 2 plume, the GAM may be able to give a first estimate on the distance up to which sources have to be taken into account for atmospheric models.A more general approach for identifying areas influenced by specific emission sources may be obtained from the wind term in Eq. ( 7) discussed in the following.

Influence of wind
As discussed above, we find wind related processes absorbed in other terms -i.e., the annual trends (Sect.4.2), the weekly term (Sect.4.3) -both representing a consistent temporal correlation with the variable modeled in these terms.However, introducing f wind as an explicit term in Eq. ( 7) still increased the accuracy of the model over large areas over the whole globe (Fig. 3, bottom).As expected, we find areas with a highly significant contribution of the wind term close to strong continuous sources, such as the east coast of the USA and east Asia, the west coast of European countries and the Arabian peninsula, and in the environments of point sources like Johannesburg and Hong Kong (see also Fig. 1).In this study, we can now illustrate the size of the areas influenced by these point sources.
As is visible from Fig. 7, we find that for locations with significant wind component, the dominant wind direction matches the main sources in the close environment.We show www.atmos-chem-phys.net/9/6459/2009/Atmos.Chem.Phys., 9, 6459-6477, 2009 the dependence on the wind direction in more detail for several points located on a circle around the point of maximal mean NO 2 TVCD (Fig. 9).One may notice that the angle corresponding to the maximum NO 2 TVCD changes according to the position.The highest NO 2 TVCD are consistently observed from the direction of the strong emission source.For a larger area around Johannesburg (South Africa), these directions are shown in Fig. 10a (white/black arrows), together with an indicator of the significance of the wind component (blue color).It is interesting to note that the location of the source corresponds to low significance of the wind component, as the wind direction is not of importance at this point, and the local source dominates the observed NO 2 distribution.Around Johannesburg the directions of the maximally contributing winds form a highly ordered vector field, even for areas with a rather low significance of the wind term.Tracking the directions -the arrows in Fig. 10a -leads to large continuous "path ways", all leading to the Johannesburg region, and indicating those regions where the central emitter is the major source for non-local NO 2 .In the Johannesburg region these zone of wind-related transportation extends in particular over the sea, in western and eastern directions, over distances of more than 1000 km.Similar dependencies on the wind-related transport and characteristic emitting point sources in the center could be found for Hong Kong, Los Angeles, Ar Riyad (Saudi Arabia), Jakarta, although they were not as sharp as for Johannesburg.
Another example for an extended emitting source of NO 2 is western Europe, as already discussed above.Areas with a highly significant contribution of wind (Fig. 10b) are found for areas such as the North Sea and north-west France/southeast England, which are located close to the main industrial areas of western Europe.For the regions of the main emission sources e.g. in Benelux/Germany and in the Po valley, the wind component is not significant, similar to the source point of the Johannesburg example (compare also Fig. 2, Milano).Here the observed NO 2 is dominated by local (constant) sources.It should be noted that the wind directions for the maximum observed NO 2 TVCD do not necessarily indicate the directions of the sources contributing the most to the observed NO 2 .
In general, the question may be raised regarding which part of the transported NO 2 is modeled when modeling the dependence on only surface winds, but not winds of higher altitudes.For the example of Johannesburg, we made the estimation for the part of the TVCD affected by the data of the surface winds.Therefore, we investigated the correlations between the winds and wind directions at different levels.The winds inclusive ECMWF level two (925 hPa, 766 m height) correlate with more than 90% with the surface winds.At level three, we still have a correlation of about 80% (see Table 2).
Since we are interested in the wind direction, we also calculated the mean absolute deviation of the wind directions.Not including level three, we can approximate NO 2 up to a height of 800 m to be exposed to winds reasonably close to the winds we used for our studies.According to numerous ground based measurement studies (e.g.Schaub et al. (2006); Wittrock et al. (2004)), this corresponds to a height below which most of the tropospheric NO 2 occurs, especially in polluted regions.Overall, for large areas -in particular in vicinity to the major industrial areas in the Northern Hemisphere -the use of the wind component increased the prediction accuracy significantly (Fig. 3).We would like to emphasize that here a more accurate model may also help to reduce bias and to get a clearer signal from the other effects, discussed at Sects.4.2-4.3.Furthermore, with the correlation between wind and observed NO 2 resulting from a local transport of the trace gas, this model extension may help to uncouple the signal from transport and local generation, providing more information about the sources of locally observed NO 2 .

Conclusions
We have successfully used generalized additive models (GAM), so far primarily to analyse local time series in order to study spatio-temporal patterns in global satellite observations of tropospheric NO 2 (i.e., NO 2 TVCD) from 1996 to 2001.Since GAM do not require regularly sampled data (as is the case for Fourier analysis), it can be easily applied to satellite observations, which usually contain gaps (e.g.caused by clouds).Thus, our results are based on daily satellite observations.This high temporal resolution and the flexibility of the additive model allowed us to systematically investigate the influence of the wind direction (average during the last 24 h) on the observed tropospheric NO 2 VCD.
We focused on studying transport processes of NO 2 plumes for polluted regions and isolated point sources.The non-parametric modeling approach allowed us to visualize the average west-eastern drift of the west-European NO 2 plume, a characteristic so far not yet reported on.When including the wind direction from ECMWF data in the model, the wind component was found to be of high importance wherever strong sources border on regions free of permanent sources, like oceans.It was also possible to identify possible pathways of atmospheric pollution and to determine the extent of areas influenced by strong NO 2 emission sources.Such information is important for the correct description of emission sources in atmospheric models.
The inclusion of the wind influence is interesting in itself, but improves also the accuracy of the results for the other relevant model terms.
Besides the influence of the wind direction, we could detect a number of features of the temporal behavior of the time series of the global tropospheric NO 2 distribution.Significant linear trends could be found in industrialized regions.Annual cycles could be mainly attributed to natural sources or variations of the atmospheric lifetime.(Results for the seasonal cycle should be treated with care as they may also be affected by other seasonal factors which do not relate to the overall NO 2 concentration (Sect.2.1).)Clear weekly cycles appeared in urban areas indicating anthropogenic sources.The location of extremes during the cycles (months of maximal NO 2 TVCD or weekdays of minimal NO 2 TVCD) gave additional insight in the kind of sources or transport processes.Most results were in general agreement with former studies using parametric rather than non-parametric model terms.
A major incentive of our study was to introduce GAMs to the global analysis of trace gas dynamics.In addition, we may summarize our main findings as follows: -The wind term increases the accuracy of a model with linear, annual, and/or weekly terms as used by, for example (van der A et al., 2008).
-We find highly consistent flow fields spatially correlating with the NO 2 time series.Tracking these flow fields allows to estimate the areas influenced by sources.The question of how these areas correlate with average transport will be left open for further (simulation) studies.
-We may have identified an approach to uncouple the part of NO 2 signal correlating with local wind fieldsand possibly with short-term transport processes -from that part of the signal which can be attributed to locally generated NO 2 .
Our study was limited with respect to several aspects, which should be improved in future studies.The spatial resolution and coverage as well as the temporal sampling of the GOME observations is rather coarse.Also the radiative transfer effects and the stratospheric correction were addressed in a rather simplified way.In the future, increasing satellite data sets with improved spatio-temporal coverage, higher spatial resolution and improved cloud correction will become available.Using such data sets will allow much more detailed studies, taking into account, for example, the wind speed, vertical wind profiles or additional quantities like temperature and precipitation.Here generalized additive models will provide ideal means for a joint, exploratory analysis of observational data from multiple sources, allowing one to access information of complex spatio-temporal patterns easily and to visualize them at a global level.

Fig. 1 .
Fig. 1.Left: Mean tropospheric vertical column density (TVCD) of NO 2 (10 15 molec/cm 2 ) around Johannesburg (top) and Hong Kong (bottom).Right: Ratio of standard deviation over MAD (median absolute deviance) for the respective regions.The high ratios between the two variability measures around the local minimum of the sources are consequences of transport processes.High values indicate that the transport is dominated by relatively few, though large events.

Fig. 2 .
Fig. 2. Examples of the GAM components for China, Central Australia, the North Sea, and the Po Valley.(a) The original TVCD time series.(b) The trend component.(c) The annual cycle.(d) The weekly cycle.(e) The wind component.The dots in (b-e) illustrate the individual measurements corrected for all other GAM components.The numbers give the GAM p-values, indicating the significance of the respective component.p-values below 1‰ are considered as significant and marked in red.

Fig. 3 .
Fig. 3.The significance (log10 of the p-value) of the different GAM terms, i.e. trend, annual cycle, weekly cycle, and wind component (from top).Dark blue areas are affected significantly by the respective term.

Fig. 4 .
Fig. 4. The strength of the different GAM terms.(a) Absolute trend in 10 14 molec/cm 2 per year.(b) Amplitude of the annual cycle (peak to peak in molec/cm 2 ).(c) Amplitude of the weekly cycle (peak to peak in molec/cm 2 ).(d) Amplitude of the wind component (peak to peak in molec/cm 2 ).Only data of p-values less than 0.001 are shown.

Fig. 5 .
Fig. 5. Month with (a) maximum and (b) minimum TVCD.Only data of p-values less than 0.001 are shown.

Fig. 6 .
Fig. 6.Day of week with (a) maximum and (b) minimum TVCD.Only data of p-values less than 0.001 are shown.

Fig. 7 .
Fig. 7. Wind direction with (a) maximum and (b) minimum TVCD.Only data of p-values less than 0.001 are shown.

Fig. 8 .
Fig. 8. Day of week with minimum TVCD (see Fig. 6b) for Europe.Over most regions, the minimum occurs on Sunday.In eastern Europe the minimum is shifted to the beginning of the week indicating the transport of the temporal emission patterns of the strong western NO x sources with the dominating westerly winds.Only data of p-values less than 0.001 are shown.Colors coded as in Fig. 6.

Fig. 9 .
Fig. 9.The dependencies of the NO 2 TVCD (red circle indicates the fitted spline) on the wind direction for 8 points around Johannesburg (see image in the center).The wind direction of maximum NO 2 TVCD (pink pointer indicate the direction of the air flux) changes accordingly to the position relative to the center.

Fig. 10 .
Fig. 10.Wind direction of maximum TVCD as quiver-plot (compare Fig. 7a) for the region (a) around Johannesburg, South Africa and (b) over Europe.Also the significance as in Fig. 3d is shown (strength of the blue color).

Table 1 .
The factor by which the retrieved tropospheric NO 2 VCD has to be corrected, depending on the height of the NO 2 profile and the solar zenith angle.Results are for cloud-free situations and an aerosol optical depth of 0.3.

Table 2 .
Comparison of the winds at different ECMWF levels to the ECMWF surface level: δθ = the mean absolute deviation of wind direction, and r = correlation.