Impact of the COVID‐19 Economic Downturn on Tropospheric Ozone Trends: An Uncertainty Weighted Data Synthesis for Quantifying Regional Anomalies Above Western North America and Europe

Abstract This study quantifies the association between the COVID‐19 economic downturn and 2020 tropospheric ozone anomalies above Europe and western North America, and their impact on long‐term trends. Anomaly detection for an atmospheric time series is usually carried out by identifying potentially aberrant data points relative to climatological values. However, detecting ozone anomalies from sparsely sampled ozonesonde profiles (once per week at most sites) is challenging due to ozone's high temporal variability. We first demonstrate the challenges for summarizing regional trends based on independent time series from multiple nearby ozone profiling stations. We then propose a novel regional‐scale anomaly detection framework based on generalized additive mixed models, which accounts for the sampling frequency and inherent data uncertainty associated with each vertical profile data set, measured by ozonesondes, lidar or commercial aircraft. This method produces a long‐term monthly time series with high vertical resolution that reports ozone anomalies from the surface to the middle‐stratosphere under a unified framework, which can be used to quantify the regional‐scale ozone anomalies during the COVID‐19 economic downturn. By incorporating extensive commercial aircraft data and frequently sampled ozonesonde profiles above Europe, we show that the complex interannual variability of ozone can be adequately captured by our modeling approach. The results show that free tropospheric ozone negative anomalies in 2020 are the most profound since the benchmark year of 1994 for both Europe and western North America, and positive trends over 1994–2019 are diminished in both regions by the 2020 anomalies.

IAGOS data are more limited in 2020 due to the reduced aircraft operations during 249 the COVID-19 pandemic, with some data available at Frankfurt airport. Clark were probably linked to decreased NOx, but that there was an important impact 252 of exceptional meteorological conditions. In the free troposphere, they found that 253 ozone levels were slightly lower than seen over the previous 26 years. In the current 254 analysis we use IAGOS data for trend detection above Europe, and we also merge 255 these profiles with those from nearby ozonesonde and lidar sites to produce continuous 256 ozone records for anomaly detection above Western Europe and western North America.   282 where y t is the ozone time series with a monthly temporal index t, β 0 is the intercept, 283 β 1 is the linear trend estimate, β 2 is the coefficient associated with ENSO index, β 3 and 284 β 4 are coefficients associated with the monthly mean zonal wind at 30 and 50 hpa (Soukharev,  The free tropospheric (700-300 hPa) ozone trend estimate and associated 2-sigma uncertainty 294 for each station are provided in Table 2. We summarize the key findings as follows:    To introduce our statistical framework, we first review the concept of mixed modeling 358 for integrating a variety of sources of (potentially heterogeneous) data (Section 3.1). We  in atmospheric sciences, we typically have too few data sources and limited samples to 376 achieve a simple average that can reconcile all of the discrepancies. Therefore, the class 377 of mixed models is the preferred approach, because the potential discrepancy from different 378 monitoring stations can be treated separately.

379
In summary, the general relationship of the mixed models can be expressed as: 380 data sources = consensus process + discrepancy + random noise.
(2) 381 However, the type of data formulation also introduces additional complexity. For example,  and to enable anomaly detection (a different goal from trend detection, but it shares the 408 same consideration discussed above), further adjustments to this framework need to be 409 made for accommodating a broader range of data heterogeneity. 410 We summarize three major adjustments as follows:  To consistently detect ozone anomalies across the entire profile, we need to make

435
We then normalized the data by using the standard deviation (SD) of deseasonalized 436 series at each pressure layer (the deseasonalized mean is expected to be zero). Deseasonalization   with a temporal index t, from station s, then the statistical model can be expressed as: 483 where f (h, t) is the consensus process which represents the underlying ozone vertical distribution 484 evolving with time (only related to h and t); g s (h, t) represents the structured discrepancy 485 (if any) from station s in addition to the consensus process; and is the random noise. 486 At this point, the above consideration treats each station (or each source of data records) 487 as equally important. In order to take the sampling frequency and inherent uncertainty 488 in each data source into account, we derive the uncertainty estimate associated with each 489 monthly normalized mean deviation using the following procedures: period, determined by each data set), so these uncertainty estimates are more homogeneous 499 and can be compared across different layers (analogous to the purpose of normalization). 500 We refer to this quantity as the "relative variability"; 501 3. Another complexity is that, similar to the monthly normalized mean deviations, 502 the relative variability can be very noisy, and directly utilizing these quantities as 503 model weights can result in noisy and unstable output. Therefore, we apply the 504 penalized regression splines to the relative variability, in order to obtain a more 505 consistent representation of uncertainty associated with each data record. 506 We use the result from the above procedures to represent the uncertainty estimates, so 507 every aggregated monthly mean deviation from each data source will have an associated 508 uncertainty estimate. We use the inverse of the squared uncertainty as the weight in the 509 model fitting process (Aitken, 1936), thus a data source with a higher sampling frequency 510 and/or lower variability receives a higher weight.

511
The model fitting procedure involves the balance of minimization of (weighted) model where · is the Euclidean norm, W is the weight matrix, D 2 F is the roughness penalty 517 term (see supplementary material for details), and f (h, t) and g s (h, t) are determined 518 through basis representations:

538
The GAMM methodology described above is now used to calculate the regional-scale ozone 539 vertical distributions, anomalies and trends above Western Europe (Section 4.1) and western 540 North America (Section 4.2), with a summary of the 2020 regional anomalies above these 541 two regions presented in Section 4.3.

571
The final fused product provides our best estimate of ozone variability above a region 572 of western Europe that covers approximately 400,000 km 2 , and therefore some of the interannual 573 variability will be impacted by variations in sampling frequency across the region. The  Figure S3 shows the diagnostics of statistical model fitting of this final fused product.

603
The residuals from the fitted model are plotted as a histogram and as a function  3. Figure S6 shows the relative variability that we used to weight the data, with the 615 high frequency IAGOS data receiving the greatest weighting, and the low frequency     Table   667 3). Overall, mid-tropospheric ozone above Western Europe increased at the rate of 0.65 668 ± 0.19 ppbv decade −1 , from 1994 to 2019, equal to a total increase of 1.6 ± 0.5 ppbv, 669 or approximately 3%. 670 Note that the results presented in Table 3 are based on our methodology that accounts  . 686 However, our updated analysis for this region spans a much longer period (1994-2020),

687
which is expected to be more robust due to an additional decade of available data.

702
It is worth mentioning that even though the western North America fused product has 703 the same number of in situ data sources as the European fused product, the overall sample 704 size (∼9,900 profiles) is much lower than that in Europe (∼45,700 profiles) (see Table   705 1). This difference in sample size might be the main reason why the fused product in Europe 706 shows more profound and detailed interannual structures, while the variability above western 707 North America is still rather indistinct between neighboring years. and their associated estimation uncertainty, can be consistently and systematically quantified.

775
An important implication from our finding is that regional trend assessments based on 776 a single data source may be less reliable due to uncertainties associated with limited data,