Diverging Ozone Trends Above Western North America: Boundary Layer Decreases Versus Free Tropospheric Increases

This study has produced an improved percentile and seasonal (median) trend estimate of free tropospheric ozone above western North America (WNA), through a data fusion of ozonesonde, lidar, commercial aircraft, and field campaign measurements. Our method combines heterogeneous data sets according to the consensus data characteristics and inherent uncertainty in order to produce our best fused product. In response to different data collection environments (in situ or ground‐based), we investigate the ozone variability based on a wide range of percentiles, which is preferable for trend detection due to tropospheric ozone's high degree of heteroscedasticity (i.e., inconsistent trends and variability between different ozone percentiles). We then compare the ozone trends and variability above the California sub‐domain to the full WNA region for better understanding of the correlations between different regional scales. In California, the 1995–2021 percentile (from the 5th to 95th) and seasonal trends are clearly positive in terms of high signal‐to‐noise ratios. The magnitude of the trends is generally weaker over WNA compared to California, but reliable positive trends can still be found between the 10th and 70th percentiles, as well as winter and summer, whereas autumn shows a negative trend over the same period. In addition, dozens of rural surface sites across the region are selected to represent the boundary layer variability. In contrast to increasing free tropospheric ozone, we find overall strong negative surface trends since 1995, with the greatest divergence found in summer. Throughout the analysis implications of the COVID‐19 economic downturn on ozone variability are discussed.

report (AR6; Gulev et al., 2021;Szopa et al., 2021) assessed output from an ensemble of chemistry climate models and concluded that the tropospheric ozone burden increased by 45% from 1850 to the present day, mainly due to the burning of fossil fuels. This modeled increase is consistent with limited observations (surface and free tropospheric) from the mid-20th century through the early 21st century (Tarasick et al., 2019). The IPCC AR6 report also concluded with high confidence that observed free tropospheric ozone has increased at northern mid-latitudes and across the tropics since the mid-1990s. A combined product from NASA's OMI (Ozone Monitoring Instrument) and MLS (Microwave Limb Sounder) satellite instruments shows positive trends of the tropospheric column ozone across the northern tropics and northern mid-latitudes over -2021(Dunn et al., 2022Ziemke et al., 2022), indicating that the broad increases of tropospheric ozone are continuing to the present day.
Meanwhile, the COVID-19 (SARS-CoV-2) pandemic emerged in late 2019 and triggered a worldwide economic downturn in 2020, which reduced emissions of ozone precursors. Studies have shown a weak but consistent ozone drop in the mid-troposphere across the Northern Hemisphere in response to the emissions reduction in 2020 (Bouarar et al., 2021;Chang et al., 2022;Clark et al., 2021;Miyazaki et al., 2021;Steinbrecht et al., 2021). Specifically, Steinbrecht et al. (2021) and Ziemke et al. (2022) observed approximately 7% or 4 ppbv ozone reductions at multiple ozone monitoring sites within the 1-8 km column during April-August 2020. By integrating ozonesonde, lidar and commercial aircraft data, Chang et al. (2022) quantified an overall free tropospheric ozone reduction of 3.6 [±1.8] ppbv above Europe and 2.8 [±1.9] ppbv above western North America (WNA) in 2020 with respect to 1994-2019. The diminished ozone levels in 2020 slightly weakened the positive free tropospheric ozone trends observed for the period 1994-2019. In terms of model studies, a 5%-15% ozone reduction was estimated for spring and summer of 2020 above the extratropical Northern Hemisphere (Bouarar et al., 2021). Miyazaki et al. (2021) estimated an ozone decrease of 5 ppbv at the surface and 3 ppbv in the mid-troposphere (500 hPa) above northern mid-latitudes in spring of 2020. Compared to the same period in March-June 2019, Miyazaki et al. (2021) also show that free tropospheric ozone (700 hPa) reductions can be observed from the CrIS (Cross-Track Infrared Sounder) satellite instrument, with one of the most remarkable reductions (∼10 ppb) found in June above the western USA. By including an additional year of data in 2021, this study will examine if the COVID-19 economic downturn has had a lasting impact on the long-term free tropospheric ozone trends.
At the surface, several studies found that ozone increased broadly in many urban areas during the lockdown period (Cooper et al., 2021;Gkatzelis et al., 2021;Keller et al., 2021;Sokhi et al., 2021), likely because reduced NO emissions tend to limit the ozone destruction in highly polluted urban areas (Sillman, 1999). Jaffe et al. (2022) analyzed surface data from 22 urban monitoring sites in the U.S. and found that COVID emission reductions significantly reduced peak ozone in the Eastern U.S., but in the western U.S. peak ozone increased in 2020 due to wildfire emissions. While those surface studies focused on urban areas, this study will investigate the COVID-19 impact on surface ozone trends and variability at rural sites across WNA.
This study describes a refined methodology for exploring trends across the full range of the ozone distribution (e.g., from the 5th percentile to the 95th percentile). Conventionally, heteroscedasticity in regression analysis refers to data variability that is not constant and varies over time. If data are homogeneous, all percentile trends should be similar. Given that previous studies generally found diverse tropospheric ozone percentile trends (Cooper et al., 2010;Gaudel et al., 2020), heteroscedasticity is an expected characteristic of ozone trend analysis. Previous investigations of ozone percentile trends have typically calculated the sample percentile at each time interval and then applied the mean or median-based trend detection technique to the percentile time series (Cooper et al., 2012;Lin et al., 2017;Michael et al., 2019;Strode et al., 2015). However, that approach implicitly assumes that the data samples are abundant enough to derive a representative sample percentile for each given month. Therefore, such an approach to study the trends in extreme events is usually limited to surface observations measured at a regular and high frequency (e.g., hourly data), instead of sparsely sampled free tropospheric ozone profiles. An alternative approach to evaluate the percentile trends is through a well-developed technique known as quantile regression (Koenker & Hallock, 2001). This method is designed to detect heterogeneous distributional changes, and has been applied to the study of surface ozone extreme events (Baur et al., 2004;Munir et al., 2012;Porter et al., 2015), and to explore extreme percentile variability from commercial aircraft data .
Whereas these two approaches discussed above have been used to successfully evaluate ozone percentile trends, the rationale behind them is slightly different. For example, if we aim to derive the 95th percentile change of an ozone time series based on daily values, the difference between these two approaches to derive the percentile trends can be understood as, respectively: 1. A 2-step method: the 95th sample percentile is first calculated from daily values for each individual month, and then the trend is fit to the resulting monthly time series (i.e., focusing on targeted sample percentiles). 2. A single-step method: quantile regression is directly applied to the full daily time series for deriving the 95th percentile trend, as it is designed for investigating distributional change (i.e., based on all available data samples).
Even though these two approaches often arrive at similar conclusions, quantile regression provides a more stable estimation when data sample sizes are deemed to be less sufficient, since all data points contribute to the estimation, rather than merely a portion of data (Krock et al., 2022). Nevertheless, the discussion so far only considers ozone time series collected from a single source (e.g., a single monitoring station). If an assessment of trends is based on multiple data sources, the complications associated with varying data characteristics should be taken into account in advance, such as spatial variability (Chang et al., 2017;Wood et al., 2017) or various measurement platforms with differing sampling frequencies .
This paper improves upon previous estimates of the free tropospheric ozone trends above WNA Cooper et al., 2010;Lin et al., 2015); through a more detailed investigation of seasonal and percentile ozone trends based on a data fusion of multiple sources of data platforms . The purpose of data fusion is to integrate multiple heterogeneous data sets to produce more consistent and accurate information. We highlight the benefit of an integration of all available data sources, which can increase our ability to detect trends in extreme percentiles for both the large domain of WNA and the much smaller sub-region of California. Further evidence of increasing free tropospheric ozone is provided by investigating ozone seasonal and percentile variability at the Mt. Bachelor Observatory, where the nighttime data are considered to be representative of the lower free troposphere, and a strong ozone increase has been previously reported (Gratz et al., 2015). An additional analysis of surface trends at rural sites across the western United States is included to explore the differences in ozone trends between the free troposphere and the boundary layer.
Section 2 outlines the sources of free tropospheric ozone profiles and rural surface measurements used in this study, and the statistical methodology for data fusion and percentile trend estimation. The results of trend analyses for the free tropospheric ozone and the rural surface ozone are provided in Sections 3 and 4, respectively. A summary of this study is provided in Section 5.

Data
The free tropospheric ozone measurements used in this study include ozonesonde, lidar, commercial aircraft, and field campaign observations. We focus on the period 1995-2021 when the data density is great enough to calculate regional-scale trends. Since different data sets report ozone values at different resolutions, measurements from each profile are aggregated into a vertical resolution of 10 hPa over 700-300 hPa, so that all data sets are aligned to the same vertical coordinate. A map view of flight tracks and data densities from the IAGOS (In-Service Aircraft for a Global Observing System) program and AJAX (Alpha Jet Atmospheric eXperiment) project, and the locations of ozonesonde and lidar stations are shown in Figure 1.
The data sets are as follows: (a) ozonesondes launched at Boulder (Colorado, 1967(Colorado, -2021, Trinidad Head (THD, California, 1997-2021, Edmonton (Canada, 1970(Canada, -2021, Kelowna (Canada, 2003 and Port Hardy (Canada, 2018(Canada, -2021, with a sampling frequency of approximately once per week (Sterling et al., 2018;Tarasick et al., 2016); (b) The lidar operated at the Jet Propulsion Laboratory Table Mountain Facility (TMF, California) has typically measured 2-5 profiles per week since 1999 (McDermid et al., 2002), with an increased frequency in the most recent years (Chouza et al., 2019); (c) The IAGOS program has measured ozone from multiple commercial aircraft since 1994 (Blot et al., 2021;Nédélec et al., 2015;Petzold et al., 2015). Ozone profiles above WNA mainly occurred above the international airports at Vancouver, San Francisco, Portland, and Los Angeles; a limitation of the IAGOS data set is the occurrence of many intermittent data gaps over this region (see Figure 1 for the annual sample size); and (d) the NASA AJAX field campaign data set includes 199 flights during 2011-2018 (Iraci et al., 2021;Yates et al., 2016Yates et al., , 2017, and the SNAAX (Scientific Aviation NASA Ames Airborne 10.1029/2022JD038090 4 of 26 eXperiment) data set includes 10 flights during 2018-2019. While most of these data sets were previously used in the recent study by Chang et al. (2022) we have been able to increase the size of the overall data set with the inclusion of the AJAX and SNAAX profiles, and additional IAGOS profiles from the southwestern USA: Denver (2), Fresno (2), Las Vegas (90), Los Angeles (404), Phoenix (24), San Jose (2) and San Diego (120). These additional 853 profiles increased the size of the overall WNA data set (10,403 total profiles) by 9%, and allowed the creation of a new aircraft ozone time series for the southwestern USA (mainly California), when merged with the 437 IAGOS profiles above San Francisco (previously used by Chang et al., 2022 in combination with IAGOS profiles above Portland, Seattle, and Vancouver). In addition, this analysis includes one additional year (2021) of data from all sites, if available, compared to the Chang et al. (2022) data set, which ended in 2020.
This part of the analysis places the focus on free tropospheric ozone variability and trends at 700-300 hPa over the period 1994-2021 in two domains (surface data are not included): 1. California: including THD ozonesondes, TMF lidar, AJAX flights, and IAGOS flights collected from San Francisco, Los Angeles, San Diego, San Jose, Fresno, and Las Vegas (this airport is geographically close Figure 1. A map view of flight tracks and data densities from the aircraft measurements above western North America , and the locations of ozonesonde (circles) and lidar (square) stations. Also shown are numbers of annual profiles from different data platforms (colored by the range of sample sizes), where (S) represents the ozonesonde station, (L) represents lidar station, and (all) represents the regional sum.
10.1029/2022JD038090 5 of 26 to California). Because we treat the AJAX and IAGOS flights as a single data set (see later discussion), this region consists of a total of three data sets; 2. WNA: includes the THD ozonesondes and TMF lidar data sets described above, plus the ozonesonde data from Boulder, Edmonton, and Kelowna/Port Hardy (these two sites are combined as a single data set), as well as other IAGOS data over WNA (e.g., Vancouver, Portland, Phoenix, Denver, and Seattle, combined with aircraft data from California as a single data set). This region consists of a total of six data sets.
Note that the study region over WNA is different from Chang et al. (2022), which did not include the AJAX data and IAGOS data collected over the Southwest U.S.. To be explicit, Chang et al. (2022) included the same ozonesonde and lidar data used in this study, but only incorporated the IAGOS data from San Francisco, Vancouver, Portland, and Seattle.
In order to better understand the difference in ozone trends between the free troposphere and the surface, several rural surface stations are included in the analysis. We use the hourly surface ozone data measured at the Mt Bachelor Observatory (MBO, 44.0°N, 121.7°W, 2,731 m), Oregon, which has been operated by the University of Washington (UW) since 2004. NOAA Global Monitoring Laboratory (GML) has operated a separate ozone instrument at the site. It should be noted that the GML instrument appears to have a positive bias in heavy smoke plumes (Bernays et al., 2022); for this reason, we use the UW ozone data, which have been processed and calibrated consistently through the data record (L. . Trend detection is based on daily time series aggregated by nighttime data only (8:00 p.m. to 7:59 a.m. local time or 4:00 a.m. to 15:59 p.m. UTC) to avoid upslope flow which transports ozone-depleted air from the surrounding valleys up to the summit. At night the localized mountain-valley circulations are diminished and MBO is more likely to sample baseline air from the Pacific Ocean, or polluted USA air that has been transported long distances (such as from California).
The analysis of rural boundary layer trends across the western USA relies on 26 high-elevation rural monitoring sites (see Table S1 in Supporting Information S1), maintained by the EPA Clean Air Status and Trends Network (CAST-NET) and the National Park Service. The daily observations were downloaded from the EPA data website (see data availability section) as maximum daily 8 hr average values (MDA8). This metric focuses on the time of day when the boundary layer is well mixed and therefore more likely to be regionally representative, and avoids the nighttime hours when the shallow nocturnal boundary layer leads to localized ozone deposition (Cooper et al., 2012). While all of these sites are located between 1.0 and 3.2 km above sea level, none are mountaintop sites. In contrast, MBO at 2.8 km above sea level is a mountaintop site, and as described earlier, the nighttime observations are made in laminar flow that lies above the nocturnal boundary layer and can, at times, be representative of the lower free troposphere or of baseline air that has been transported to the summit from the Pacific Ocean (Fischer et al., 2011;Gratz et al., 2015). However, the site is known to be heavily influenced by forest fires in summer and air masses at the site may have been influenced by the U.S. boundary layer during the previous day or days. Therefore, while we expect MBO to sample a different set of conditions and air masses from the other 26 rural surface sites, and while it should be more strongly influenced by the lower free troposphere, some influence from the U.S. boundary layer is still likely to occur. The links for all the data sets are provided in the data availability section.

Data Fusion of Free Tropospheric Measurements From Multiple Platforms
Deriving the regional ozone trends in the free troposphere based on infrequent observations at sparsely distributed monitoring stations is challenging, due to multiple sources of uncertainty associated with limited data Cooper et al., 2010;Lin et al., 2015). In our previous work, we developed a statistical framework for combining various non-independent vertical profile records, such as data sets from commercial aircraft, lidar, and ozonesonde measurements . The methodology is designed to, (a) determine the contribution of each data source through its deseasonalized and normalized deviations (ND), according to the sampling frequency and inherent data uncertainty; and (b) borrow the strength from neighboring (vertically) correlated structure to identify systematic ozone variability between data sets, remove unstructured variations within each data set, and thus achieve an enhancement of signal-to-noise ratio .
The rationale of data fusion developed in Chang et al. (2022) can be summarized by the following statistical structure: various data sources = fused consensus process + discrepancy from individual dataset + random noise (see Chang et al., 2022 for technical details) This methodology is implemented through the framework of the generalized additive models (Hastie & Tibshirani, 1990;Wood, 2006), by assigning separate penalized regression spline functions to the consensus process and discrepancy term. The regression splines are chosen so that the systematic variations between different data sets can be constituted (Wood et al., 2017). Therefore, this methodology is designed to separate common regional variations from each individual data source.
For additional implementation details, it should be noted that the methodology in Chang et al. (2022) was only applied to monthly mean deviations above WNA, but the variation in extreme percentiles is of particular interest to the research community, due to their relevance for impacts on human health and vegetation (Fleming et al., 2018;Mills et al., 2018). However, the statistical properties of extreme percentiles can be very different from the sample average (Berrocal et al., 2014;Eastoe & Tawn, 2009;Smith, 1989), thus the same methodology is not directly applicable. More specifically, Chang et al. (2022) used an inverse of standard error of the sample mean to represent the uncertainty associated with sample frequency and data variability, but this statistic is not necessarily representative of the uncertainty of the extreme percentiles (especially ozone data, which is often non-normally distributed), and might not be appropriate for weighting extreme percentiles from different sources of data sets. Even though the uncertainty of extreme percentiles can be assessed through bootstrap-based approaches (Gilleland, 2020;Kyselỳ, 2008), the weighting scheme for extreme values is currently absent and has not been thoroughly investigated. Therefore, in this study, we treat each data set equally when fusing the percentiles from different data sets.
While the equal-weight strategy for different data sets may seem to be less sophisticated, it is a desirable feature in this study. Recall that six data sets are used for data fusion, including Edmonton (53.6°N), Kelowna (49.5°N) combined with Port Hardy (50.4°N), THD (41.1°N), Boulder (40.0°N), TMF (34.4°N), and aircraft data combined from the IAGOS, AJAX and SNAAX (their sampling latitudes vary over time). Therefore 5 out of 6 data sets are sampled at fixed latitudes (accounting for 83% of the weights in the data fusion, the impact of the remaining 17% aircraft data on regional percentile trends will be discussed later). This weighting scheme preserves sufficient influence from each site spread over different latitudes. For example, the number of TMF lidar profiles in 2021 is roughly equal to the sum of the ozonesonde profiles added from Edmonton, Port Hardy, Boulder, and THD, so the equal-weight strategy avoids an out-sized influence from TMF. It should be emphasized that even though we use the equal-weight scheme in the algorithm, this does not imply each data set has the same contribution to the fused product, since the fused product will be determined from the common systematic variations between different data sets.
Spatial heterogeneity is an unavoidable complication in modeling regional ozone trends. The most straightforward approach is to explicitly account for the spatial variability (e.g., Chang et al., 2017Chang et al., , 2021. However, since the spatial coverages of free tropospheric observations vary greatly between different years above WNA (see Figure 1), aircraft data and observations at a few fixed locations are generally insufficient to determine the full spatial patterns. In terms of spatial representativeness, Liu et al. (2009) found the free tropospheric ozone observations are correlated over distances of 500-1,000 km. On the other hand, Stauffer et al. (2022) recently examined the stability of global ozonesonde network and concluded that those data are suitable to study trends (albeit the accuracy depends on the sampling frequency in each station), but trends at the individual station would much likely reflect localized variability (and partial regional variability). Therefore it is important to consider all available data sources as a whole for better characterizing large-scale variability Cooper et al., 2010;Stauffer et al., 2022). Even though spatial heterogeneity is not fully accounted for in our current methodology, we evaluated the sensitivity of our data fusion method to spatial heterogeneity by also calculating the free-tropospheric ozone trends using a geographically weighted modeling approach. This second method is described in Section S1 in Supporting Information S1 and the results are discussed in Section 3.3, demonstrating that similar trends are produced by the data fusion and the geographically weighted methods.

Percentile Trend Estimations
As discussed above, conventional trend analyses of tropospheric ozone time series have been carried out using monthly, seasonal, or annual aggregated values (depending on the ozone metrics). The statistical purpose of such aggregations is to diminish heterogeneous data variability (e.g., hourly or daily values are expected to be even more variable and heterogeneous than monthly values), allowing traditional hypothesis tests to gain statistical power. On the other hand, in response to the need for understanding the trends of extreme events, previous attempts were often made by first deriving a single percentile per temporal scale, and then applying the same 7 of 26 trend model to the resulting ranked value time series. However, when we consider the above two issues together, the well-developed quantile regression method becomes a practical solution for both challenges. Quantile regression has been shown to be a well-suited approach for detecting heterogeneous distributional changes (Chang et al., 2021). The percentile trends can be quantified by directly applying the methodology to daily surface data or sparsely sampled ozone profiles, which is advantageous because no further data aggregation is needed, and we can keep as much information as possible for improved analysis of the extreme events .
The trend estimation is based on the following model, with adjustments for the potential correlation of the El Niño-Southern Oscillation (ENSO) and quasi-biennial oscillation (QBO) with free tropospheric ozone Neu et al., 2014;Ziemke et al., 2019): where y t is the ozone time series with a temporal index t, β 0 is an intercept, β 1 is the linear trend estimate, {β 2 , β 3 , β 4 } is a set of coefficients associated with ENSO and QBO at 30/50 hPa, respectively, and ɛ t is the random noise (assumed to be an AR(1) process). Note that the data are deseasonalized in advance, so the seasonal component is not included in Equation 1. Equation 1 describes how the statistical relationship is formulated, the next step is to provide context for evaluation of the trend and regression coefficients by median-and quantile-based methods. Median regression (also known as least absolute deviations, LAD) is well known for its robustness to outliers and non-normally distributed data, as compared to the least squares (mean-based) method. Quantile regression extends LAD by enabling quantification of other percentile changes (Koenker & Hallock, 2001). If we rewrite Equation 1 as ɛ t = y t − X t β, where t = 1, …, T, then the percentile estimate for the trend and each regression coefficient can be evaluated by minimizing the following weighted least absolute error criterion: where 0 < q < 1 represents the quantile to be estimated (e.g., q = 0.5 represents the median).
Although autocorrelation tends to cause an underestimation of uncertainty associated with the trend estimate (Weatherhead et al., 1998), the influence of daily autocorrelation on the long-term trends is expected to be rather weak. On the other hand, bootstrapping methods are known to produce a consistently greater estimation uncertainty than standard methods, when data are heteroscedastic and/or non-independent (Gonçalves & White, 2005;Politis et al., 1997). Therefore, in order to account for as much variability as possible, we adopt a bootstrapping approach to determine the standard error associated with the trend (i.e., uncertainty estimate) from quantile regression (Koenker, 1994). In this study, we use SNR (signal-to-noise ratio, defined as the ratio between trend and its uncertainty) value to assess the reliability of the trend estimate. An SNR value of 2 (or 3) represents the threshold of approximately the 95% (or 99%) confidence level. The higher the absolute SNR value, the higher the confidence of the trend estimate.
It should be emphasized that in this study quantile regression is only applied to surface measurements at each monitoring site. For the free tropospheric ozone observations from different platforms, we first account for the heterogeneity between different data sets. Therefore, the data fusion described in Section 2.2.1 is applied to the profile data from each data set in order to produce the fused ozone percentile (vertical) distribution over time, and then the mean trend estimates are derived from the resulting fused percentile product. To be consistent with the uncertainty from quantile regression, the uncertainty of the trend estimate for the fused percentile time series is determined by the bootstrapping method implemented by the R package boot (Canty, 2002;Kunsch, 1989).

Preliminary Analysis
To ensure the field campaign data can be incorporated into the trend assessment, we evaluated the comparability between different data sets to avoid suspicious jumps or offsets, which can bias the trends. Field campaign data are usually scheduled for certain years, which makes it difficult to derive reliable long-term monthly means to deseasonalize the data. Since our methodology fuses data through the deseasonalized and ND, we assume the AJAX data (in the free troposphere) can be deseasonalized by the seasonal cycle based on the IAGOS data. Figure  S1 in Supporting Information S1 shows the boxplots of annual ozone observations and their deseasonalized anomalies from different data sets in California over 2011-2020. We can see some discrepancies between the annual distributions of the AJAX (as well as IAGOS) and ozonesonde/lidar data (e.g., 2014 and 2017), mainly because the AJAX data are not uniformly sampled for all months, but after deseasonalization, the discrepancies are diminished. This result indicates that the AJAX data are compatible with other data sets and can be incorporated into the trend assessment. Despite the fact that the IAGOS data set has a large gap over 2006-2016 in California, incorporation of the AJAX data set partially fills the data gap and produces a more accurate trend estimate.
Trend analysis without taking seasonality into account (or deseasonalization) is not only prone to estimation bias (when data coverage is limited in some years), but also much more likely to produce an inflated trend uncertainty (regular patterns should not be considered to be a part of trend uncertainty). Figure S2 in Supporting Information S1 shows the median ozone time series from each data set, with a linear trend estimate and a Loess (locally estimated scatterplot smoothing, Cleveland, 1979) nonlinear fit to the deseasonalized anomalies. Even though different linear trends and uncertainties are found for different data sets (due to localized differences), the general variability is consistent among time series, except for the early period where the confidence intervals between IAGOS and ozonesonde data are separated for a short period of time (but this does not adversely affect the trend assessment, as indicated by a sensitivity test of trends derived from different starting years; see Section 3.3). Indeed, even though IAGOS can be considered to be a reference data set due to its known calibration history (Tarasick et al., 2019), and ozonesonde data are of high quality and stability (Stauffer et al., 2022), Tarasick et al. (2019) concluded that ozonesonde observations are moderately higher (5%-8%) compared to IAGOS data. This mean difference could produce a bias in the trend estimate, but deseasonalization of each data set removes this offset.
The merged annual ozone data set and its anomaly distributions based on all available data over California are shown in Figure 2. We can see how the ozone distribution varies across the years, and this plot also demonstrates the importance of deseasonalization for time series analysis, as it transforms non-normal distributions (presumably due to a non-uniform or sporadic sampling schemes) into normal-alike distributions (see the distribution in 2003 and 2004 as an example). However, a positive skewness can still be observed in some deseasonalized anomaly distributions; in this case, a difference between mean and median trends might be expected (i.e., heteroscedasticity), and an investigation of the extreme percentiles is desirable for better understanding of a greater extent of the ozone variability.

Mean Trends
Since average ozone mixing ratios and their variations typically increase with altitude, we transform the ozone profile data into ND by deseasonalizing and normalizing ozone measurements at each 10 hPa pressure surface, so that the data are more comparable across different pressure levels and different data sets. This consideration is critical because statistical models perform better when data have similar levels of variability. The data fusion methodology is then applied to the deviations from each data set. Figure 3 shows the fused ozone distributions in the free troposphere, based on the ND and a transformation back to the units of ppbv, over California and WNA respectively. The plot of ND enables us to enhance the delicate ozone variability across different pressure surfaces; we can observe that the general patterns are well correlated between California and WNA in terms of the distribution of interannual positive and negative anomalies.
To avoid the impact of incomplete data in 1994 above California and anomalously low ozone in 2020 Steinbrecht et al., 2021), the trend estimates are reported for three periods 1995-2019, 1995-2020, and 1995-2021, and based on the column average over 700-300 hPa of the fused monthly anomalies (e.g., the last two panels of Figure 3). The positive trends above California are found to be slightly stronger than WNA, but the interannual variability is similar in general between the two regions. The California data are a subset (50%) of the WNA data set, and cover a much smaller region. The fact that similar interannual variability and SNR (signal-to-noise ratios) can be derived from the California sub-domain indicates that the observed ozone increase is a robust feature of the WNA domain (see Table 1).
To highlight the fact that intermittent field campaign data available in certain years can still add value to regional trend detection, Figure S4 in Supporting Information S1 shows the fused ozone distributions above California, with and without the AJAX data. With additional data included in the analysis, we can see a much more detailed variability in 2011-2018 with an enhanced SNR value in trend detection.
10.1029/2022JD038090 9 of 26 To compare with our previous result , we provide the trend estimates with reference to 1994 in Table S2 in Supporting Information S1. Our previous regional mean trend estimate of free tropospheric ozone over WNA was 0.4 [± 0.2] ppbv/decade (1994-2019), but as mentioned in Section 2.1, our previous result did not include the AJAX data or the IAGOS data over the Southwest (i.e., data collected from Las Vegas, Los Angeles, and San Diego). The regional mean trend above WNA in this study is much stronger (1.0 [± 0.3] ppbv/decade) over the same period. We conclude that the discrepancy with our previous trend estimate is due to a greater ozone enhancement in the second half of the study period revealed by additional data above the Southwest.
The following analysis places the focus on the comparison of a specific portion of ozone variability (i.e., percentile and seasonal trends), but we should keep in mind that the overall trends remain fairly consistent between California and WNA.

Percentile Trends
The transport and chemical processes that drive heterogeneous trends are typically different for high and low ozone in the free troposphere. The high ozone values are more likely affected by stratosphere-troposphere exchange processes (Langford, 1999), long-range transport of air pollution plumes (Nowak et al., 2004), or wildfires ( 10.1029/2022JD038090 10 of 26 Davies et al., 1998;Grant et al., 2000). The data fusion methodology is applied to the 5th, 50th, and 95th percentiles, and the resulting fused distributions above California are provided in Figure 4: the 95th percentile trend is positive and stronger than the 50th percentile, and then followed by the 5th percentile. To better compare the variability from different percentiles, the curtain plots based on units of ppbv are scaled to their relative percentile values (instead of zero centered anomalies as in Figure 3). Noticeable common interannual variabilities include.
• The enhanced ozone episode over 1997-1998, previously attributed to enhanced stratosphere-troposphere exchange following a strong El Niño event (Koumoutsaris et al., 2008;Langford, 1999), but also explained by extreme biomass burning in southeast Asia (Fiore et al., 2022); • An enhancement in 2015 in the 50th and 5th percentiles. An attribution study of this free tropospheric enhancement has not yet been conducted, however, ozone enhancements at the surface and in the boundary layer of the western USA have been linked to a high-temperature anomaly and increased biogenic emissions ); • Strong ozone anomalies in 2018 in the 95th, 50th and 5th percentiles. An attribution study of this free tropospheric enhancement has not yet been conducted. A previous study found strong correlation between enhanced ozone and fires for 2017, 2018, and 2020 in urban areas over WNA , and we suggest that a future attribution study could explore the possibility that wildfire smoke impacted free tropospheric ozone.
The corresponding fused distribution over WNA (which includes the California subdomain) is shown in Figure 5: the magnitude of the trend in all three percentiles is weaker with respect to California, and the 5th percentile trend is nearly flat over 1995-2021, but the interannual variability shows good agreement with the California subdomain. Figure S5 in Supporting Information S1 shows the changes of trends presented by 10 hPa resolution vertical profiles for both regions based on the ending year of 2019, 2020, and 2021, respectively. The vertical structures of the percentile trends do not reveal any obvious differences due to the choice of ending year. Therefore, we summarize the overall free tropospheric ozone trends (700-300 hPa) in all three percentiles in the first part of Note. SNR (signal-to-noise ratio) is the ratio between trend and its standard error. Table 1. We can observe a coherent drop of trends for both regions when 2020 data are included, except for the 95th percentile above California (which is presumably associated with extensive wildfires in 2020 as discussed above). When the study period is extended to 2021, we can observe increasing positive trends in the 95th percentiles, no change in the median trends, and much weaker positive trends in the 5th percentile.
To investigate the heterogeneity of regional percentile trends, we further extend the analysis by applying the same statistical methodology to every 5-percentile interval between the 5th and 95th percentiles. The resulting distributions of percentile trends are shown in Figure 6. The corresponding curtain plots of ozone distributions by every 10th percentile are shown in Figure S6 in Supporting Information S1 (above California) and S7 (above WNA). We summarize the findings as follows: • At first glance, we can see that neither California nor WNA shows a strictly monotonic or homogeneous structure of percentile trend distributions, which might indicate that the 5th, 50th, and 95th percentiles may not be sufficient to fully represent the heterogeneous ozone variability. For example, the strongest positive trends over WNA are found within the 10th-20th percentile range, but the 5th percentile trends are nearly flat over 1995-2021. • Prior to the COVID-19 economic downturn (1995-2019), nearly all percentiles above California and WNA were reliably positive. Inclusion of the pandemic years  shows that the positive trends in the lower percentiles (5th-25th percentiles) above California were greatly diminished. • Above California, all trends between the 95th and 5th percentiles remain reliably positive (SNR > 2, i.e., the lower bound of the 2-σ interval greater than zero) after the 2020 COVID-19 impact. Above WNA, despite the reduction in the lower percentile trends, the general pattern of percentile trend distribution has not changed (e.g., trends in the 10th-70th percentile range remain reliably positive.  Figure 7 summarizes the overall quantified annual percentile anomalies in the free troposphere, and the merged monthly median time series (with observed fluctuations), above California and WNA, respectively. Note that the magnitudes of the 95th and 5th percentiles are average differences from the medians. The annual median anomalies are strongly correlated (r = 0.84) between California and WNA, followed by the 5th (r = 0.72) and 95th  (r = 0.66) percentiles. Note that the uncertainty for the 95th percentile is constantly greater than the 50th and 5th percentiles, this is not unexpected due to a positive skewness and a heavy tail of the ozone distribution.
As pointed out previously, ozonesonde and lidar data have a fixed weight in the data fusion (accounting for 83% of the regional contribution), therefore a 17% weight is contributed by aircraft data sampled at varying latitudes. Figure S8 in Supporting Information S1 shows the annual data density made by the sampling latitudes from the IAGOS and AJAX data. Since 2016 there are fewer data sampled from high latitudes (over >45°N), which could potentially create a bias by being heavily weighted toward low latitude data. Nevertheless, the downward tendency of the 5th percentile anomalies over 2015-2020 can be clearly observed above California (the top panel of Figure 7), but not observed over a greater domain above WNA, this indicates the regional trends over WNA should not be heavily impacted by imbalanced sampling latitudes from aircraft data in recent years.
To examine the reliability of our methodology for deriving fused percentile trends, in Section S1 in Supporting Information S1 we adopt a geographically weighted modeling approach to derive the overall trends by an adjustment to the spatial heterogeneity, and in Section S2 in Supporting Information S1 we apply different approaches for estimating the regional ozone percentiles (Fasiolo et al., 2020;Hyndman & Fan, 1996). Both results show great similarity to the discussion above, confirming the effectiveness of our approach. Results from a further sensitivity analysis are presented in Figure S12 in Supporting Information S1 to investigate the percentile trend estimates and uncertainties based on different starting years (1995, 1998, 2000, and 2005) and ending years (2019 and 2021). The year 1998 is selected as it represents a well-known enhanced ozone episode that was observed across the Northern Hemisphere, and has been recently attributed to extreme biomass burning (Fiore et al., 2022). In general, the big picture remains similar, except that, (a) the lower percentile trends above California are shifted from positive over 1995-2021 toward negative over 2005-2021 (as discussed above, lower percentile trends in California are heavily affected by the pandemic years); and (b) if the estimates are based on the starting year in 1998 with unusually high ozone, its impact can be observed mainly on the reduction of trends at the higher percentiles above WNA, but the sensitivity is weak at the lower percentiles. Nevertheless, with the exception of lower percentile trends above California, the trends and uncertainties are comparable regardless of the starting year (1995, 2000, or 2005), and we conclude that our result is robust.

Seasonal Trends
Figure 8 provides seasonal ozone distributions and trends (1995-2021) based on the fused medians from different data sets above California and WNA (the corresponding curtain plots for ND are provided in Figure S13 in Supporting Information S1). In California, all seasons show positive trends (very high confidence, SNR > 4), with the strongest trend observed in JJA (June-July-August), followed by MAM (March-April-May), DJF (December-January-February), and SON (September-October-November). In WNA, we find strong positive trends for DJF and JJA (high confidence, SNR > 2), weak positive trends for MAM, and strong negative trends for SON. In both regions, the most confident (i.e., the highest SNR value) positive trends occur in DJF. The magnitude and SNR of the trends above WNA are weaker than California.
A notable feature in Figure S13 in Supporting Information S1 and Figure 8 is that the interannual variability is more dynamic in MAM and DJF. To explore the reason behind this phenomenon, additional analyses are provided as follows: • We investigate the stratospheric influence on tropospheric ozone by showing the percentage (from all available profiles) of ozone values exceeding a threshold at 100,125,150,175,200, and 225 ppbv above WNA, respectively in Figure S14 in Supporting Information S1; such high ozone values in the free troposphere are explained by stratospheric intrusions. The frequency of observing ozone exceedance across 700-500 hPa seems to be more common in recent years than the late 1990s. We made similar plots focusing on MAM and JJA in Figure S15 in Supporting Information S1: For a threshold of 150 ppbv and in the 350-300 hPa layer, 42 months (50%) in JJA and 10 months (12%) in MAM have zero exceedances over 1995-2021. For the months with non-zero exceedance rates, an averaged exceedance rate in MAM is also higher than JJA. The greater frequency of stratospheric intrusions in spring compared to summer is consistent with many previous studies (Cooper et al., 2014;Stohl et al., 2003). Given that these high ozone values occur in very few profiles in any given month, any impact on ozone trends would be limited to the higher percentiles. • We extracted meteorological output at San Diego, THD, and Vancouver from the NASA Global Modeling and Assimilation Office's MERRA-2 (Modern-Era Retrospective analysis for Research and Applications, Version 2) reanalysis (Gelaro et al., 2017), including seasonal geopotential height and specific humidity at 500 hPa ( Figure S16 in Supporting Information S1) and monthly mean and standard deviation of potential vorticity (PV) at 300 hPa ( Figure S17 in Supporting Information S1). The variability of PV at all three locations is highest in DJF and MAM, indicating a greater frequency of upper-level long wave troughs at this time of year, which are commonly associated with stratospheric intrusions and therefore high ozone variability, consistent with the extreme ozone observations described above. Geopotential height and specific humidity are greatest during summer at all three locations, indicating a higher frequency of upper-level ridges and moist tropical air extending into mid-latitudes. These conditions are consistent with a northward migration of long-wave troughs to high latitudes, minimizing the occurrence of stratospheric intrusions. These dynamical indicators point to conditions in summer that are relatively quiescent with increased tropical characteristics, which limit ozone variability.
The lower half of Table 1 reports the seasonal trend estimates for both regions. The seasonal positive trends over 1995-2019 diminish after the 2020 data are included, but overall, the impact of the 2020 COVID-19 economic downturn on tropospheric ozone trends is slight in SON above California and more pronounced in other seasons. When the 2021 data are included, the diminished trends persist in DJF and MAM, but there was a rebound of the positive trend in JJA (consistent with the findings from three NASA satellite products (Ziemke et al., 2022)).
Above WNA Cooper et al. (2010) found the median trend of free tropospheric ozone in springtime to be 0.

Result: Rural Surface Ozone Trends
The previous section focused on a data fusion of the free tropospheric ozone profiles from several different data sets, with trend estimates derived from the fused product. This section will focus on surface ozone observations and trends.

Surface Trends at the Mt Bachelor High Elevation Site
As described above, heterogeneous ozone trends and variability at different percentiles were found in the free troposphere in association with the different data platforms, which had to be taken into account when applying the data fusion technique. However, if the time series data are from a single source (e.g., the same location and same measurement technique), quantile regression is a straightforward and preferable approach for trend detection of ozone observations. Even though the applications of quantile regression have gained popularity in the trend detection of atmospheric composition and pollutants (Baur et al., 2004;Gaudel et al., 2020;Munir et al., 2012;Sousa et al., 2009;Xu & Lin, 2018), the practical adjustments for atmospheric time series analysis, such as the potential impact of inconsistent seasonality at different percentiles on the percentile trend estimate, have not been investigated. Therefore, we provide a comprehensive technical evaluation of quantile regression specifications in Supporting Information S1 as follows: • Section S3 Supporting Information S1 demonstrates that deseasonalization is a critical consideration for accurate percentile trend detections. As a demonstration of quantile regression, Figure 9 shows the (deseasonalized) percentile time series at MBO with their mean trend estimates, as well as daily time series with percentile trend estimates derived by quantile regression. We can see that these two approaches produce a consistent median trend estimate, but some discrepancies are found at the 5th and 95hr percentile trends, likely due to a short record being considered here (18 yr). These discrepancies are expected to be reconciled when longterm data are available.
• A series of sensitivity tests is presented in Section S4 in Supporting Information S1 to demonstrate that (moderate) intermittent data gaps are expected to have a small impact on the percentile trend estimates.
Since the above inconsistency in the extreme percentile trends is very likely due to a short data record, and since quantile regression is tailored to the purpose of detecting distributional changes (Reich, 2012), we therefore adopt quantile regression to produce more robust percentile estimates. Section S5 in Supporting Information S1 provides the details on how an attribution of data variability can be performed at different percentiles.
To investigate the impact of the 2020 COVID-19 economic downturn on ozone trends at MBO, Table 2 provides the trend estimate, 2-σ uncertainty and SNR value for the 5th, 50th, and 95 percentiles based on three periods, 2004-2019, 2004-2020, and 2004-2021, respectively. The positive trends over 2004-2019 are clearly weakened when 2020 data are included, with a 50% reduction in mean and a 37% reduction in median. We also see that the weakened trends rebounded in 2021, especially for the 95th percentile. Figure 10 shows the percentile trends over the above three periods, from the 5th to 95th percentile by every 5 percentile interval. In addition, it is well known that wildfires have an impact on ozone variability at MBO in summer (Baylon et al., 2015;Laing et al., 2016), thus we also show the same analysis, but exclude the wildfire season (June-September). The 2020 COVID-19 impact on percentile trends can be clearly observed in both scenarios: a stronger reduction is found in the lower percentiles and a weaker change is found in the higher percentiles (consistent with the free tropospheric results above WNA and California). Another noteworthy feature is that the positive trends in the 95th and 90th percentiles became much weaker over 2004-2021 when the data in the wildfire season were removed, indicating that the rebounding 95th percentile trend is very likely driven by ozone production associated with the extensive fires in 2021.

Surface Trends at Rural Sites Across the Western USA
Twenty-six high-elevation, rural monitoring sites across the western United States were selected to represent long-term, regional-scale ozone variability in the boundary layer; see Table S1 in Supporting Information S1 for site details (station name, state, longitude, latitude, and elevation). Some of these sites were previously analyzed by Cooper et al. (2012), who found increasing median ozone trends in MAM, and a range of trends in JJA for the period 1990-2010. Note. SNR (signal-to-noise ratio) is the ratio between the trend and its standard error. Percentile trends for each site are shown in Figure S24 in Supporting Information S1 over 1995-2021 (any sites with measurements beginning after 2005 are excluded from this particular plot, so a total of 21 stations are shown). The trend values and their associated uncertainties are based on the 2000-2019 reference period and quantile regression of MDA8 ozone. The results show that 12 (57%), 17 (81%), and 20 (95%) sites have reliably strong decreasing trends in the 5th, 50th, and 95th percentiles, respectively, so that an overall decrease of regional ozone is not unexpected from this consistent pattern of negative trends across WNA.
To better summarize the regional variations across the domain, we applied a geographically weighted regression model to all available time series data, based on a Gaussian spatial process model fitted through the generalized additive models (Chang et al., 2017(Chang et al., , 2021. The fitted result can be considered to be the overall common trend after an adjustment of the irregularity of the spatial distribution of monitoring stations. The resulting regional 5th, 50th, and 95th percentiles displayed by individual years from 1995 to 2021 are shown in Figure 11 (their continuous time series plots are shown in Figure S25 in Supporting Information S1). A clear ozone reduction was observed in MAM and JJA of 2020 in all of these percentiles, except for a strong enhancement in the 95th percentile that coincides with the wildfires in August 2020 (Filonchyk et al., 2022). A clear rebound of ozone in 2021 to 2000-2019 climatological mean can also be observed (except for MAM where the ozone level is still lower than the mean over 2000-2019). Over 1995-2021 the 95th percentiles have the strongest negative trends, followed by the 50th and 5th percentiles ( Figure S25 in Supporting Information S1).
The percentile and seasonal trend estimates for the regional rural surface sites are provided in  . When comparing the seasonal trends before and after 2020, weak changes in JJA/SON may reflect the transition from emissions reductions during the lockdown to ozone enhancement due to the wildfires in the late summer of 2020 (such a reflection within 2020 can also be observed in Figure 11).

Summary of Differences in Trends Between the Free Troposphere and Rural Surface Measurements
A highly consistent pattern in Table 1 is that all of the percentile trends show a substantial difference between the free troposphere and rural surface when considering all months in the year, that is, opposite in magnitudes and no overlap of the 2−σ uncertainties. For seasonal trends, the largest difference in trends is found in JJA, followed by MAM, and similar trends are only found in SON. This finding adds to the growing body of evidence that surface trends are frequently disconnected from the general increases observed in the free troposphere , as assessed by IPCC (Gulev et al., 2021), and also reproduced by a recent chemistry-climate model simulation (Fiore et al., 2022). In order to summarize the overall difference in (median) trends between the free troposphere and the rural surface over WNA, we compare the trends at the rural surface sites to the free tropospheric trends above California and WNA, based on individual months for the periods 1995-2021 and 2004-2021 (with the MBO record included) in Figure 12. From the top panel (1995-2021) we can clearly observe opposite trends in summertime and inseparable trends in wintertime between the free troposphere and rural surface, which clearly indicates that different processes must be affecting the decreasing ozone in the WNA boundary layer compared to the increasing ozone in the free troposphere. When we limit the study period to a shorter time frame (2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020)(2021), similar seasonal patterns can still be observed. Consistent positive trends can be observed at MBO, with much stronger enhancements in August-October presumably due to the wildfire impacts (Baylon et al., 2015;Jaffe et al., 2022;Laing et al., 2016). A better agreement of trends between MBO and the free troposphere above WNA/California is observed in April-July.

Discussion
This study provides an extensive analysis of tropospheric ozone trends and variability over WNA, using (a) free tropospheric ozone profiles from ozonesondes, TMF lidar, IAGOS commercial aircraft, and the NASA AJAX field project; and (b) surface measurements from MBO, and the National Park Service and EPA CASTNET monitoring networks. Although we present the results starting with the free troposphere and then proceed to the rural surface sites across WNA, the trends are summarized for three different scenarios: 1. A surface time series measured at a fixed location: Since surface station data are usually measured at a high frequency (i.e., hourly values are available, albeit subject to data gaps), applying quantile regression to daily aggregated values is appropriate for studying the extreme events. 2. Multiple surface time series collected from a regional monitoring network: For each predetermined percentile, the spatial heterogeneity (from this specific percentile variability) and the irregular distribution of monitoring locations need to be taken into account to produce the common regional time series, and then the overall trends can be derived. 3. Multiple time series made by sparsely sampled vertical profiles from different platforms: For each predetermined percentile, various data sets are combined (via a data fusion methodology) by taking into account the vertical correlation structures and data characteristics (from this specific percentile variability). The fused result is produced as the ozone vertical distribution over time, and the overall trends are derived from the column averages of the fused product.
All the percentile trends examined in this study show heterogeneity (i.e., consistent trends are not found between different percentiles and seasons). This result suggests that only using the mean and a few percentiles (e.g., 5th and 95th) may be insufficient to fully represent the heterogeneous ozone variability, given the fact that the distribution of percentile trends is not monotonic between the 5th and 95th percentiles above WNA. As a candidate for median trend estimators, quantile regression has the same advantage as the classical Sen-Theil estimator, such as robustness to outliers and no assumption on error distribution. As a part of regression-based methods, quantile regression has more flexible model structures over the classical non-parametric methods, such as incorporation of covariates for trend attribution. We show that percentile trend estimates based on quantile regression are resistant to moderate data gaps (even with a short data record), and consistent with the result estimated through individual ranked values (if data samples are sufficient). Therefore, quantile regression is considered to be the preferred approach, especially when data samples are sparse.
The trend analysis of the fused free tropospheric observations and rural surface network in this study can be summarized as follows: 1. Despite being impacted by the 2020 COVID-19 economic downturn, median ozone trends remained positive in the free troposphere over 1995-2021. The interannual variabilities of median ozone above California and WNA are well correlated, but the magnitude of trends is much stronger above California (1.4 [±0.3] ppbv/ decade) than WNA (0.5 [±0.2] ppbv/decade). 2. In the free troposphere above California, nearly all percentile and seasonal trends are stronger and more confident than above WNA. The 2020 COVID-19 impact can be clearly observed from weakening trends in the lower and middle percentiles above California. After the 2021 data are included, we found that the positive 5th percentile was further diminished, with a rebounding positive trend in the 95th percentile, and no change in the median trend. The percentile and seasonal trends remain consistently positive (high confidence, SNR > 2) since 1995 through the end of 2021. 3. Weak positive mean and median trends can be detected reliably in the free troposphere above WNA over 1995-2021. While positive trends were also observed in the 5th and 95th percentiles, the confidence level is much lower; positive trends can also be found consistently (high confidence) between the 10th-70th percentile range. Even though weakened positive trends in the lower percentiles are found after the 2020 COVID-19 pandemic, the overall distribution of percentile trends remains similar. 4. After the 2020 COVID-19 pandemic, continuously weakened positive trends are found in DJF and MAM, and rebounding positive trends are found in JJA through 2021 in the free troposphere above WNA (as well as California The trends in all the percentiles show rebounding ozone from 2020 to 2021. 6. Twenty-six high-elevation, rural sites were selected from the EPA CASTNET and National Park Service air quality monitoring networks to represent boundary layer variability across the western U.S. Despite sharp ozone reductions in MAM and JJA in 2020, the regional mean, median, 5th and 95th percentile trends changed little and remained strongly negative through the end of 2021 (very high confidence, |SNR| > 3). The 95th percentile trends (−1.4 [±0.4] ppbv/decade) were roughly twice as strong as the trend of the 5th percentile (−0.6 [±0.3] ppbv/decade).

7.
Finally, yet importantly, overall opposite percentile and seasonal trends (no overlap of the 2-σ ranges) are found between the free troposphere and the rural surface sites across WNA, except for SON which shows consistent negative trends. By focusing on individual monthly results, we are able to further confirm that the greatest difference occurred in summertime, which suggests that different mechanisms are responsible for increasing ozone in the free troposphere and decreasing ozone in the rural regions.
From 2020 to 2021, even though we found rebounding trends in JJA and continuously diminished trends in DJF/MAM, the long-term mean and median trends show no noticeable change in the free troposphere above WNA. The only structural change we found after the pandemic period is the reduction of trends in the lower percentiles in the California free troposphere and MBO, which do not show evidence of rebound in 2021.
In terms of surface ozone at rural sites, albeit with very weak signals, a similar seasonal variability can be detected (rebounding in JJA and declining in DJF/MAM), but the overall mean and median trends remain unchanged.
Based on previous work the most likely explanation for the opposing trends in the free troposphere and at the surface is a sharp contrast in ozone precursor emissions trends in the U.S. boundary layer compared to upwind regions (Fiore et al., 2022). Several regional-scale studies based on observations and chemical transport models have shown that surface ozone is decreasing across most of the United States due to reductions of domestic ozone precursors, primarily NOx (Cooper et al., 2012;Hogrefe et al., 2011;Jaffe et al., 2018;Luo et al., 2020;Simon et al., 2015;Strode et al., 2015). On the hemispheric scale, a strong shift of ozone precursor emissions from high latitudes to low latitudes has increased ozone production in the tropics and subtropics of the Northern Hemisphere and the export of these precursors and resulting ozone to mid-latitudes has contributed to the observed increase of ozone in the free troposphere (Christiansen et al., 2022;Fiore et al., 2022;Gaudel et al., 2020;Y. Zhang et al., 2016Y. Zhang et al., , 2021, with additional contributions from the increase in commercial aircraft emissions (Wang et al., 2022), and a small contribution from the global increase of methane (Lin et al., 2017). Our present analysis, which focuses on observed ozone trends, is the first of a series of papers to explore ozone trends and attribution in the free-and lower-troposphere above WNA. We have work in progress to apply numerous satellite products and global atmospheric chemistry models to better characterize the source-attribution and main physiochemical processes driving the observed free-and lower-tropospheric trends.

Conflict of Interest
The authors declare no conflicts of interest relevant to this study.