A long-term (2001–2022) examination of surface ozone concentrations in Tucson, Arizona

Abstract

Ground-level ozone (O₃) pollution in semi-arid regions like Tucson, Arizona, presents unique challenges due to the interplay of anthropogenic emissions, biogenic volatile organic compounds (BVOCs), meteorological conditions, and regional transport. Tucson is the second-largest city in Arizona and has received comparatively less attention than the most populated city of Phoenix despite experiencing elevated O₃ levels amid rapid population growth. This study provides a comprehensive 22 year analysis (2001–2022) of O₃ trends in Tucson using a combination of ground-based monitoring data, satellite observations, NEI emissions inventories, land cover classification and meteorological datasets. The findings reveal no statistically significant long-term trend in O₃ levels at northwest (NW), urban core, and south/southeast (S/SE) monitoring sites despite regulatory actions to reduce precursor levels. However, spatial differences persist with one S/SE site (Saguaro National Park) consistently exhibiting the highest O₃ concentrations and an urban core site (Rose Elementary) usually exhibiting the lowest values across all seasons. Satellite and surface-based data reveal a decline in NO₂ across the study period, in contrast to HCHO levels that show little long-term change, with a brief increase in 2020 likely linked to regional fire activity and higher temperatures, particularly in June. Consequently, FNR values (formaldehyde-to-NO₂ ratio) increased after 2005–2009, indicating a regional shift influenced by reductions in NO_x emissions, especially during fall/winter and spring. This shift helps explain the weakening of the weekend effect (i.e., higher weekend levels versus weekdays) over time and the emergence of the weekday effect earlier in the summer (June) in contrast to the late 1990s. Generalized additive model meteorology normalization suggests that 79% of the O₃ variability is attributed to interannual weather variability. FNR started to decline post-2020, suggesting changes in O₃ responsiveness to further NO₂ reductions, particularly in cooler months. These dynamics, along with recent fall/winter O₃ increases, highlight the complex, chemical regime-dependent response of O₃ to precursor changes. This study recommends improved VOC characterization to inform future air quality strategies in the region.

---

Environmental significance

Ground-level ozone pollution remains a pressing air quality challenge in semi-arid cities like Tucson, Arizona, where rapid urban growth intersects with complex emissions and meteorological dynamics. Despite regulatory efforts to reduce precursors, ozone levels show no long-term decline, with persistent spatial differences and recent seasonal changes in ozone behavior. This study reveals a chemical regime shift marked by changing formaldehyde-to-NO₂ ratios (FNR), weakening the traditional weekend effect, and altering ozone sensitivity, especially in cooler months. These findings underscore the critical need for improved, speciated VOC characterization to better constrain VOC-NO_x sensitivity. Understanding these evolving dynamics is essential for developing effective, region-specific air quality strategies in the face of climate change and ongoing urban expansion.

1 Introduction

Reductions in ambient O₃ levels for the U.S. in recent decades have been pursued through control of its key precursors, NO_x and volatile organic compounds (VOCs).^1,2 Enhanced surface level O₃ persists as a problem though, particularly in the semi-arid southwestern United States, where cities face unique challenges due to a combination of factors including their native climate, geographical features, wildfires, and rapid urbanization.^3,4 High temperatures and intense sunlight in Arizona cities such as Phoenix and Tucson promote O₃ formation, particularly in the summer. Phoenix and Tucson, with populations of 5 070 110 and 1 063 162,⁵ respectively, both face air quality challenges linked to these conditions. While much O₃ research has focused on the larger city of Phoenix,^6–9 less emphasis has been placed on Tucson. An updated view of Tucson's O₃ characteristics is needed as this city is also very populated and faces similar challenges as well due to regional factors, a growing population, and extreme heat.

In Tucson, foundational work by Diem and Comrie¹⁰ revealed the pronounced role of biogenic VOC (BVOC) emissions stemming from local vegetation within Tucson's arid landscape. Diem¹¹ observed that during April–September 1995–1998, monsoon-driven BVOC emissions (e.g., isoprene and monoterpenes) increased sufficiently to reverse Tucson's chemical regime from VOC-limited to NO_x-limited, a pattern reinforced by temperature-enhanced BVOC production seen in other cities.¹² By mapping NO_x and VOC emissions across Tucson, Diem and Comrie¹³ identified pronounced pollution gradients that helped explain why O₃ levels differed so considerably between the city's urban core and its rural sites. Abraham and Comrie¹⁴ developed a real-time O₃ mapping approach that blends regression and interpolation techniques to stitch together scattered monitoring data, directly linking emissions patterns to dynamic O₃ variability. Together, these studies and others^15,16 revealed how Tucson's O₃ levels hinged not just on emissions, but also on the complexities of arid meteorology including temperature and mixing height variability. Simon et al.² noted that relying on outdated studies risks having incorrect information about current air quality characteristics since O₃ non-attainment areas that were previously VOC-limited (e.g., Phoenix) years ago could have transitioned to being more in a NO_x-limited regime in summer months when O₃ exceedances are most likely. It is essential to understand whether an area is VOC- or NO_x-limited because inappropriate NO_x regulations could increase O₃ inadvertently.¹⁷

A relevant recent study examining rural and urban sites in Arizona between 2015–2021 showed that the urban core of Phoenix experiences VOC-limited regime conditions during winter and fall and shifts towards a transitional and NO_x-limited regime during spring and summer.⁷ It remains unclear how the chemical regimes change in Tucson. Expanding the lens to the broader Southwest, studies by Cooper et al.^18,19 and Mousavinezhad et al.²⁰ show that O₃ levels throughout the southwest U.S. are increasing due to emissions combined with regional transport. Similarly, Parrish et al.²¹ highlight that background O₃ levels, including contributions from wildfires, now dominate over anthropogenic emissions in this region.

This study presents a refreshed view of Tucson's O₃ characteristics while also considering a historical perspective by examining 22 years of data (2001 to 2022) beginning with near the end (∼2003) of when past studies examined O₃ in Tucson.⁴ Unlike prior research, which primarily focused on the summer months, this work incorporates additional monitoring sites and examines data spanning all seasons. This broader seasonal perspective allows for a more complete understanding of O₃ trends and variability across different times of the year. This study is guided by the following questions: (i) how has ground-level O₃ in Tucson changed from 2001 to 2022?; (ii) how do those changes vary based on location in Tucson and for different seasons?; and (iii) how do O₃ chemical regimes vary seasonally and over time, and how are these regimes related to O₃ trends?

2 Methods

2.1 Study area

Tucson, situated in Pima County in southern Arizona (32.253° N, 111.113° W), is the second-largest city in the state with a population that increased by 27% from 828 905 in 1999 to 1 063 162 in 2023 based on Pima County population data.⁵ Tucson sits at an elevation of roughly 700 m above sea level and is approximately 180 km southeast of the largest city, Phoenix. Tucson is uniquely situated within a basin encircled by mountain ranges including the Santa Catalina, Santa Rita, Rincon, Tucson, and Tortolita Mountains, which lie approximately 10 to over 40 km from the city center with peaks ranging from ∼1400 m to nearly 2800 m. As a hot and semi-arid city, Tucson is prone to naturally occurring windblown dust, as well as emissions from vehicular traffic, suburban development, and industrial emissions.^22,23 Tucson experiences limited precipitation, with the most occurring during the summer monsoon season,¹¹ as also shown in Fig. S1. Despite the dry conditions, the Tucson metropolitan region has considerable leaf biomass^10,11 as shown in Fig. S2. Its “urban forest” blends native and non-native species, including mesquite, palo verde, and eucalyptus. The vegetation in desert areas is predominantly composed of shrubs, including bursage and creosote bush,²⁴ as well as trees like palo verde, acacia, and mesquite.^11,25,26

Land cover classification data from 2001 and 2022 reveal a clear expansion of developed land in and around Tucson (Fig. S2a and b), with corresponding changes in vegetative cover over time. The land cover change map (Fig. S2c) highlights areas of increased development (shown in red), suggesting that urban growth has reshaped the landscape in ways that may influence emissions. Since 1999, eight long-term O₃ monitoring sites have been operational throughout the metropolitan area (Fig. 1). These include sites in northwest (NW) Tucson (Tangerine and Coachline), urban core sites (Children's Park, Craycroft, Rose Elementary), and to the south/southeast (S/SE) of the inner city (Saguaro National Park, Fairgrounds, Green Valley).

Fig. 1 Map showing the eight O₃ monitoring sites classified as northwest (NW) (blue oval), urban core (purple oval), and south/southeast (S/SE) (red oval). Next to each site is a O₃ pollution rose diagram based on 13 years (2010–2022) of data over Tucson and representative of the full 24 hours of a day. Black lines are arrows connecting each monitoring site to its corresponding pollution rose diagram, and thin gray lines represent basemap roads for orientation only.

2.2 Datasets

2.2.1 U.S. EPA Air Quality System (AQS). Hourly and maximum daily 8 hour average (MDA8) O₃ were obtained from the U.S. EPA Air Quality System (https://aqs.epa.gov/aqsweb/airdata/download_files.html), covering the period from 1 January 2001 to 31 December 2022 for the eight stations shown in Fig. 1. Hourly NO₂ data were obtained for Craycroft and Children's Park sites from 2010 to 2022, as these were the only locations with consistent data availability. Data values reported as zero or below were excluded, as in past work.²⁷ Values reported below the method detection level (MDL) were substituted with one-half of the MDL for hourly data.²⁸ The MDL represents the lowest concentration reliably distinguished from zero with 99% confidence.²⁹

2.2.2 U.S. EPA emissions inventory data (NEI). Triennial emissions data from the U.S. EPA National Emissions Inventory (NEI; https://www.epa.gov/air-emissions-inventories/get-air-emissions-data-0) are used to assess long-term trends in pollutant emissions that may help explain observed changes in surface O₃ concentrations. County-level totals for NO_x, VOCs, and carbon monoxide (CO) were extracted for Pima County for the years 2008, 2011, 2014, 2017, and 2020. Emissions were categorized by source sector, including on-road, non-road, point, fire, and biogenic categories.

2.2.3 Meteorological data. Temperature, wind speed, and wind direction data were obtained from the U.S. EPA AQS (https://aqs.epa.gov/aqsweb/airdata/download_files.html) for the same eight locations in Fig. 1, covering the period from 1 January 2001 to 31 December 2022. PBLH data were downloaded at 3 hour time resolution from the Modern Era-Retrospective Analysis for Research and Applications, version 2 (MERRA-2) at a spatial resolution of 0.5° × 0.625° (ref. 30) covering the period from 2010 to 2022. PBLH is considered as it affects surface concentrations and local pollutant dispersion.

Precipitation data are obtained from the Precipitation Estimation from Remotely Sensed Information using Artificial Neural Networks (PERSIANN) system,³¹ which spans latitudes from 60° S to 60° N with a resolution of 0.25° × 0.25°. To examine precipitation trends, daily precipitation data within the spatial bounds of the Tucson metropolitan area (111.12716° W, 110.737116° W, 31.87952° N, 32.425261° N) were obtained for the period 2001 to 2022.

For an in-depth examination of meteorological trends and their impact on pollutant distribution, see SI Fig. S1 and the related discussion in Section S1. Briefly to summarize key points, the summer months are marked by elevated temperatures and peak precipitation during the monsoon season (July–August), in contrast to the cooler, drier winters. PBLH most follows the monthly air temperature trends.

2.2.4 Satellite data. This study uses data from the O₃ Monitoring Instrument (OMI) for the historical period (2005–2020) and Tropospheric Monitoring Instrument (TROPOMI) for recent years (2020–2022). NO₂ and HCHO concentrations (2005–2020) from the OMI instrument were used to examine the formaldehyde to NO₂ ratio (FNR) trends. OMI quantifies backscattered solar radiation to analyze atmospheric composition, with data at a spatial resolution of 0.25° × 0.25°. The L3 NO₂ vertical column density (VCD) global gridded data used here include modifications such as application of refined air mass factors.³² The HCHO vertical column L2 data developed under the QA4ECV3 project^33,34 were regridded to a resolution of 0.25° × 0.25° using area-weighted averaging to align its spatial resolution with other datasets.

Additionally, we utilized TROPOMI Level 2 satellite data for the most recent three years in the study period: 2020–2022. The Sentinel-5P satellite's TROPOMI instrument measures ultraviolet, visible, near-infrared, and shortwave infrared wavelengths with a spatial resolution of about 5.5 km by 3.5 km.^35,36 HCHO values were filtered based on data quality assessments to a quality assurance value exceeding 0.50 and NO₂ values over 0.75 (ref. 33 and 37) after gridding the data to a resolution of 0.07° × 0.07° following Mirrezaei et al.,³⁸ enabling detailed analysis and facilitating site-specific evaluations at individual monitoring locations. Chemical regimes are sometimes classified using documented FNR thresholds (<3.2 for VOC-limited, 3.2–4.1 for transition, and >4.2 for NO_x-limited) from;¹² however, we caution that those definitions are most applicable to high-NO_x environments (i.e., megacities) and thus not as relevant for Tucson.

2.2.5 Land cover and land use data. To evaluate the impact of urbanization on O₃, we used land cover data from the National Land Cover Database (NLCD), obtained from the Multi-Resolution Land Characteristics Consortium viewer (https://www.mrlc.gov/viewer/). Classified land cover maps are examined for 2001 and 2022 at 30 m spatial resolution to quantify changes in land use across the Tucson metropolitan region. A difference map between 2022 and 2001 was generated to identify areas where shrubland, grassland, or barren surfaces were converted into developed classes (low, medium, or high intensity). These spatial changes were then compared with O₃ monitoring site categories and used to contextualize observed shifts in O₃ levels and weekly patterns.

For another quantitative perspective on land type, we utilized monthly 1 km Moderate Resolution Imaging Spectroradiometer (MODIS) Normalized Difference Vegetation Index (NDVI) data obtained from NASA's Application for Extracting and Exploring Analysis Ready Samples (AppEEARS; https://appeears.earthdatacloud.nasa.gov/task/area). NDVI values were used for correlation analysis with air-quality variables.

2.3 Calculations

MDA8 O₃ values were utilized for daily, monthly, and seasonal assessments. Data from 2001 to 2022 were utilized to evaluate long-term trends in O₃. Satellite-derived NO₂, VOC, and FNR ratio data were analyzed for the period 2005–2022, while surface NO₂ measurements from two monitoring sites (2006–2022) were used to corroborate the satellite NO₂ trends. Time series analysis was conducted before categorizing the data into four intervals (2001–2009, 2010–2019, 2020, and 2021–2022), with 2020 separated to account for the COVID-19 pandemic. The Mann–Kendall trend test was applied to detect monotonic trends and determine statistical significance. This non-parametric test is optimal for environmental data as it does not assume a normal distribution and is resilient to outliers. Statistically significant differences were determined using p-values, with a threshold of p < 0.05 indicating a significant trend. These statistical analyses aid in ascertaining whether detected discrepancies stem from authentic environmental changes or random data variations.

We report two complementary O₃ metrics: (i) seasonal medians of daily MDA8 O₃ to describe the central tendency of ambient conditions, and (ii) ozone design values (ODV; the 3 year average of the annual 4th highest MDA8 O₃) to quantify the regulatory upper end, relevant to NAAQS attainment. Part of this study focuses on day-of-week differences in O₃, for which mean O₃ concentrations for the four time periods were calculated for weekdays (Monday–Friday) and weekends (Saturday–Sunday), allowing for an examination of the weekly cycle. In contrast to some studies, such as those by Koplitz et al.³⁹ and Simon et al.,² which focused on specific days such as Tuesday to Thursday for weekdays and Sunday for weekends, our study encompasses all days of the week. This decision was informed by comparing both approaches, which indicated no significant differences in outcomes between selecting specific days and including all days. Seasonal analysis is conducted using the following months: winter (DJF), spring (MAM), dry summer (J), monsoon summer (JA), and fall (SON). Note that we separate the summer into two periods as the monsoon period has important effects on air pollutants such as O₃ (ref. 7) and particulate matter due partly to enhanced BVOC emissions.^11,40–42 For simplicity of presentation of results, we also combine winter and fall seasons due to their similar results as in other work for the study region.⁷

To assess temporal changes in NO₂, HCHO, and FNR over time, we computed both climatological medians for each season and also seasonal anomalies. The former represents the median of each variable's value for a particular season over either the range of 2005–2020 (for OMI) or 2020–2022 (for TROPOMI). The seasonal anomaly is conducted separately for OMI and TROPOMI periods and is the difference between a season's median value for a single or subset of years and the overall climatological median. These anomalies are visualized as heatmaps and help interpret shifts in O₃ chemical regimes over time.

To quantify the relationship between land cover, emissions, and O₃, we calculated Pearson and Spearman correlations, using monthly de-seasonalized anomalies of median NDVI, MDA8 O₃, along with median NO₂ and HCHO columns. De-seasonalized anomalies were derived by subtracting the long-term monthly climatology for each variable using data between 2001 and 2022. Statistical significance was assessed with three p value thresholds: p <0.05, <0.01, and <0.001.

2.4 Generalized additive models (GAM) for de-weathering trend analysis

To quantify the influence of meteorological factors on MDA8 O₃ and isolating emission factors, we applied a meteorological normalization procedure relying on Generalized Additive Models (GAMs).^43,44 GAMs are a flexible extension of generalized linear models that allow both linear and non-linear terms, making them suitable for modelling complex, nonlinear interactions between meteorology and air quality.^45–48

We followed approaches used in previous de-weathering studies^45–50 We combined daily MDA8 O₃ data from eight Tucson monitoring sites into a single domain dataset using the median value. These aggregated data values were then modelled as the dependent variable, with input predictors including daily values of temperature (°C), relative humidity (%), precipitation (mm), PBLH (m), wind speed (m s⁻¹), wind direction (°), and also day of year (DOY), day of week (DOW), and month. Wind speed and direction were transformed into orthogonal u- and v-components to represent east–west and north–south flow, respectively. The GAM was fitted for Tucson as a domain using smooth functions s(x) for continuous variables and factor terms f(x) for categorical variables:

O₃ = s(T) + s(RH) + s(u_wind) + s(v_wind) + s(precipitation) + s(PBLH) + s(DOY) + f(DOW) + f(month) + ε (1)

where ε is the residual error term.

Model fitting was implemented using the pyGAM library in Python, selecting smoothing parameters via grid search to optimize the generalized cross-validation score. Meteorology-normalized O₃ was estimated by replacing daily meteorological predictor values with their site-specific day-of-year climatological means while retaining the observed temporal and categorical predictors (DOY, DOW, month). This generated a weather-normalized prediction for each day that reflects conditions under average meteorology. The difference between the actual GAM prediction and the normalized prediction represents the estimated meteorological effect on O₃ for that day. This approach enables separation of emission-driven and meteorology-driven variability in O₃, allowing for more robust trend detection and attribution.

3. Results

3.1 Annual O₃ design values for cumulative dataset

To provide context for the subsequent analysis, brief discussion begins with ODVs from 2001 to 2022 across monitoring sites in the Tucson region (Table 1) due to their relevance for policy. Values exhibit some similar characteristics including fairly steady values from 2001 till 2008, before experiencing a decline until about 2017 likely due to regulatory activities reducing precursor levels, especially NO_x. After 2017, ODVs began to stabilize and increase slightly. An exception to these trends was Saguaro National Park (SNP), which experienced increasing ODVs initially from 2001 through 2005 unlike other sites. We note that ODV emphasizes the upper end as it is sensitive to the frequency and magnitude of high O₃ episodes even when seasonal medians remain relatively stable.

Table 1 Ozone design values (ppb) for each location by year from 2001 to 2022

Year	Tangerine	Coachline	Children's Park	Rose Elementary	Craycroft	Saguaro National Park	Fairgrounds	Green Valley
2001	71	NA	72	NA	72	70	69	NA
2002	72	NA	73	NA	73	73	71	NA
2003	72	NA	73	NA	72	74	69	NA
2004	72	NA	72	NA	72	76	69	NA
2005	72	65	73	65	72	77	69	67
2006	72	68	72	66	71	76	68	68
2007	73	67	73	68	70	76	71	68
2008	72	68	71	67	68	74	70	66
2009	69	65	69	65	66	71	69	65
2010	69	64	67	65	65	70	68	65
2011	69	63	67	65	66	70	69	67
2012	69	65	66	67	68	71	70	69
2013	68	66	68	68	67	73	71	70
2014	67	65	67	66	64	71	68	68
2015	65	63	67	65	63	70	66	64
2016	65	61	64	64	63	68	64	63
2017	66	64	65	66	64	69	67	64
2018	67	66	67	66	65	71	68	66
2019	68	67	68	64	67	70	68	64
2020	68	67	69	63	68	69	67	64
2021	68	66	68	63	68	68	67	64
2022	69	67	70	65	70	69	68	67

---

Starting in 2003, SNP consistently recorded the highest values in the region, approaching or exceeding the current 8 hour ozone NAAQS threshold of 70 ppb. This highlights SNP as a critical area of concern. The data also suggest a narrowing O₃ gradient over time between urban core and the NW and S/SE sites. For example, while peripheral sites like Green Valley and Coachline initially reported lower ODVs than urban core sites such as Craycroft and Children's Park, this difference has diminished, pointing to more regionally uniform O₃ levels in recent years. Overall, these trends underscore both the progress made in reducing O₃ pollution and the ongoing challenges marked by the ODV rebound in recent years that warrants increased attention to factors limiting continuing reductions in surface-level O₃.

3.2 Interannual seasonal trends for cumulative dataset

This section discusses seasonal interannual trends in MDA8 O₃ concentrations from 2001 to 2022 across Tucson's NW, urban core, and S/SE sites. Results for individual sites are in Table 2 and a graphical view of the three site category results are in Fig. 2. Fig. S3 additionally shows the results for each site in Table 2 in the form of spatial maps. Assisted by prevailing wind patterns from the urban core, the SNP site among the S/SE sites consistently exhibits the highest O₃ levels with Green Valley exhibiting the lowest coincident with prevailing winds from the west and southwest where there are no major urban emissions. SNP shows the highest levels of any site considered, with Fairgrounds also usually near the highest levels after SNP. Overall, the S/SE site mean values tend to be highest, and this is most pronounced in fall/winter seasons where the next higher levels are in NW sites and then urban core locations.

Table 2 MDA8 median O₃ (ppb) concentrations for each location, by season and across four binned periods. Mn summer refers to monsoon summer months of July–August and dry summer is June

Binned years	Seasons	Tangerine	Coachline	Children's Park	Rose Elementary	Craycroft	Saguaro National Park	Fairgrounds	Green Valley
2001–2009	Winter/fall	40.74	36.02	35.40	34.83	34.99	42.18	40.98	39.63
	Spring	53.90	52.08	52.28	49.51	51.18	54.80	52.64	51.97
	Dry summer	52.95	49.85	53.08	47.50	52.95	55.84	51.37	49.78
	Mn summer	53.06	48.61	52.35	47.67	51.28	55.95	52.00	48.50
2010–2019	Winter/fall	40.12	36.72	36.78	36.05	36.80	42.10	40.16	40.13
	Spring	53.01	51.11	51.81	51.22	51.02	55.33	53.21	53.20
	Dry summer	51.51	48.33	50.86	49.09	50.63	55.17	52.26	49.86
	Mn summer	50.59	47.47	49.52	47.38	50.01	52.84	50.76	48.00
2020	Winter/fall	42.68	39.75	41.31	39.92	42.23	43.85	42.96	43.15
	Spring	49.69	49.38	49.85	46.08	49.69	51.31	51.15	49.69
	Dry summer	48.60	48.80	49.60	45.75	52.00	54.50	52.25	50.75
	Mn summer	51.00	50.50	53.00	49.12	54.50	53.12	52.25	49.88
2021–2022	Winter/fall	41.85	39.94	40.58	39.07	41.65	43.40	43.17	41.28
	Spring	51.44	51.88	52.58	47.85	52.92	53.61	52.35	51.70
	Dry summer	52.62	51.33	54.11	49.86	53.00	54.50	54.29	50.87
	Mn summer	52.82	51.11	55.00	52.27	54.88	54.06	54.41	50.78

---

Fig. 2 Time series of seasonal MDA8 O₃ concentrations for NW, urban core, and S/SE sites from 2001–2022. The solid line indicates median values and shaded areas represent interquartile ranges (25th to 75th percentile).

The fall/winter median MDA8 O₃ levels are on average about ∼2.19 and ∼3.61 ppb higher for S/SE sites relative to NW and urban core sites, respectively. This pattern is consistent with observations from other metropolitan areas such as New York,⁵¹ and also Atlanta, Chicago, Dallas–Fort Worth, and Phoenix, where downwind sites recorded mean peak O₃ concentrations 5–20 ppb higher than nonurban upwind locations during comparable periods.⁵² The elevated O₃ concentrations at Tucson's S/SE sites are consistent with how afternoon O₃ levels are higher in downwind areas of the Phoenix area and correlated with morning time NO₂ photolysis rates at upwind areas.⁹ Thus, daytime transport likely contributes to the spatial trend of O₃ in Tucson leading to peak values usually at SNP and Fairgrounds.

Other notable features in Table 2 are that Rose Elementary exhibits usually the lowest O₃ levels even though it is in the urban core. A potential explanation is that its O₃ pollution rose (Fig. 1) is different from the other urban core sites in that it has more influence from westerly winds rather than northwesterly winds that presumably transport more of the emissions from the urban core and upwind areas in the NW cluster and even areas farther upwind like Phoenix. Among the NW sites, Tangerine usually exhibited higher levels than Coachline even though they are in close proximity and one possible reason is that Coachline is positioned by a major source of vehicular missions (Interstate 10) with prevailing daytime winds being directed from Interstate 10 to the northwest towards this site. Because the difference in median levels between the two sites is most substantial in winter/fall, this may point to the influence of increased NO_x emissions from vehicles leading to O₃ reductions assuming this area is in a VOC-limited regime during those seasons as suggested by Greenslade et al.⁷ and explored further in a subsequent section. However, this is speculative and requires more investigation.

The Mann–Kendall trend test results indicate that there is generally no statistically significant variation (p > 0.05) over the study period regardless of site and season, suggesting that the observed fluctuations are largely driven by natural variability rather than a persistent long-term trend. This result does not rule out anthropogenic influences; rather, competing drivers, including precursor emission reductions, background O₃, and meteorological variability, likely offset one another, producing little net change. Exceptions include a significant positive trend observed in the urban core during the fall/winter season (tau = 0.44, p < 0.001). Specifically, Children's Park (tau = 0.57, p < 0.001) and Rose Elementary (tau = 0.39, p < 0.02) exhibited significant increasing trends.

3.3 Day of week patterns

One way to understand O₃ concentrations in populated cities is to examine day-of-week (DOW) variations of O₃ and its precursors. Such methods capitalize on differences in NO_x emissions between weekdays and weekends due to fluctuating vehicle traffic among other sources.^53,54 DOW O₃ patterns are well-documented in the literature for cities such as Los Angeles,^55,56 Phoenix,^6,7 and New York City,^2,57 among several urban and rural areas.^52,58 Building on the foundational work of Diem,¹¹ who identified distinct weekday–weekend O₃ patterns in Tucson during the 1990s, this study provides an updated analysis of seasonal O₃ dynamics in the region. We examine the seasonally-resolved weekly cycle of MDA8 O₃ for the four binned periods (2001–2009, 2010–2019, 2020, 2021–2022) across NW, urban core, and S/SE sites. Fig. 3 shows results for individual days of the week whereas Fig. 4 shows spatially how the difference between weekend and weekday MDA8 O₃ levels have changed over time for each site and season.

Fig. 3 Seasonal weekly cycle of MDA8 O₃ concentrations for NW, urban core, and S/SE sites for four binned periods (2001–2009, 2010–2019, 2020, and 2021–2022). The solid line indicates median values, and dash lines represent interquartile ranges (25th to 75th percentile).

Fig. 4 Seasonal variation of weekend–weekday differences of O₃ among distinct site categories—NW, urban core, and S/SE—over four-time intervals (2001–2009, 2010–2019, 2020, 2021–2022). The period 2001–2009 serves as the baseline and is not compared to any earlier interval, which is why it shows no relative change. For panels in the farthest three columns, each triangle indicates the percent change in O₃ relative to the immediately preceding interval (i.e., 2010–2019 is compared to 2001–2009, 2020 is compared to 2010–2019, and 2021–2022 is compared to 2020). Upward-facing triangles indicate positive percentages in O₃ concentrations, while downward-facing triangles indicate negative percentages relative to the respective baseline period. Blue shading indicates a weekday effect (weekday > weekend); red shading indicates a weekend effect (weekend > weekday).

During the fall/winter seasons, all three site categories exhibited a clear weekend peak in O₃ concentrations, typically reaching the highest levels on Sundays and the lowest during weekdays as is consistent with Greenslade et al.⁷ for urban Arizona sites. This pattern is characteristic of a VOC-limited atmospheric regime in an urban area whereby reduced NO_x emissions on weekend days due to reduced anthropogenic activity like driving can result in an increase in O₃. However, in 2020, this pattern changed, with the weekend–weekday difference dropping to the point that O₃ peaked on weekdays with the urban core category showing the clearest ramp up in O₃ towards Wednesday and a gradual decline towards Friday and the weekend days. In the subsequent period (2021–2022), the pattern reversed again with weekend–weekday differences increasing with a general reduction in O₃ on weekdays and increased levels around Friday–Sunday (peak usually on Sunday). Fig. 4 shows how urban core and NW sites have a more pronounced weekend effect (i.e., higher weekend–weekday difference) outside of 2020. Across sites, DOW differences were significant (p < 0.05) at all sites in 2001–2009, mixed in 2010–2019 (not significant at Children's Park and Coachline), and not significant at any site in 2020 or 2021–2022. Spring initially exhibited a weekend effect during the earliest binned period (2001–2009), consistent with the Diem¹¹ observation of a persistent weekend effect in April in downtown Tucson, which gradually expanded to more sites across Tucson in May and June, based on data from 1995 to 1998. However, in subsequent periods (2010–2019 and 2020), the spring pattern shifted to a weekday effect before eventually reverting to a weekend effect in 2021–2022. In 2020, coinciding with the onset of COVID-19 lockdowns in March, O₃ concentrations decreased markedly, coincident with reduced anthropogenic activity and lowered regional scale O₃ levels as noted by Greenslade et al.⁷ For spring, most sites showed significant (p < 0.05) weekend–weekday differences in 2001–2009, except Green Valley, Fairgrounds, and Saguaro National Park, whereas several sites showed no significant differences in 2010–2019, except Children's Park, Tangerine, and Coachline. In 2020, only Children's Park and Tangerine showed a significant weekend–weekday difference. In 2021–2022, all sites were significant except Coachline.

Dry summer (June) consistently showed a weekday effect throughout the study period, contrasting with Diem (2000) who reported a weekend effect during the same month. For June, Diem (2000) noted that only one site in their analysis showed a significant weekend effect, and that was SNP, which stands out in our findings as well for consistently lacking a weekend effect during the dry summer period (Fig. S4). A notable result is that the dry summer season of 2020 shows a distinctly different weekly profile with a pronounced ramp-up in O₃ from Sunday to Wednesday–Thursday for all three site categories. This pattern is further corroborated by weekly surface NO₂ measurements in Children's Park and Craycroft (Fig. S5), which exhibit a similar weekly progression, reinforcing the link between weekday-emission patterns and O₃ formation during that period. This seems to have lingered into 2021–2022 albeit slightly less pronounced. In this season, no sites showed statistically significant weekend–weekday differences between 2001-2009 and 2021–2022. In 2010–2019, most sites showed significant weekend–weekday differences except Craycroft and Children's Park. In 2020, most sites showed significant weekend–weekday differences except Saguaro National Park and Green Valley.

Monsoon summers predominantly displayed a weekday effect, aligning with findings from Diem¹¹ for July and August. This pattern is also consistent with findings by Buysse et al.,⁵⁹ who documented a weekday effect in Sequoia National Park (California) during spring and summer under high and moderate temperatures. The observed weekday effect in our summer results aligns with the decreasing trend of the weekend effect reported in other regions such as for Los Angeles after 2000,^12,60 particularly between 2011 and 2015. This shift from a traditional weekend peak to a weekday peak underscores changing emission patterns and atmospheric dynamics, reflecting a broader trend beyond Tucson. For July–August, no sites showed significant weekend–weekday differences in 2001–2009 or 2020. In 2010–2019, only Fairgrounds showed significant weekend–weekday differences, and in 2021–2022, most sites exhibited significant (p < 0.05) differences except Children's Park, Green Valley, and Saguaro National Park.

3.4 Precursors and chemical regimes

We next examine satellite-based OMI (Fig. 5) and TROPOMI data (Fig. 6) for HCHO, NO₂, and FNR (HCHO/NO₂) to better interpret results from Sections 3.1–3.3. We do not have routine surface VOC observations for Tucson over this period; consequently, we interpret satellite HCHO columns as a proxy for reactive VOCs and use NEI data for context about emissions types. Fig. 5 and 6 show climatological seasonal medians and also anomalies in each season (described in Section 2.3). Beginning with the OMI period (2005–2020), NO₂ anomalies are systematically higher in the earliest years (2005–2009) and also generally highest in fall/winter (Fig. 5). Consequently, NO₂ anomalies are negative between 2010–2019 and 2020 with the largest reduction in fall/winter. There is a statistically significant decreasing trend in NO₂ across all seasons between 2005 and 2020 with Kendall tau values for fall/winter, spring, dry summer, and monsoon summer being −0.87, −0.83, −0.50, and −0.75 (p < 0.05), respectively, consistent with decreasing NO₂ emissions due from improved emissions control and regulatory policies.¹ Fig. S6 further supports these results based on surface measurements of two urban core sites (Children's Park and Craycroft) with a notable decline in NO₂ from 2010 to 2020. However, levels in 2021 and 2022 appear slightly elevated compared to 2020, suggesting a partial rebound in NO₂. TROPOMI data (Fig. 6a) are in line with surface observations, showing mostly positive NO₂ anomalies during 2021 and 2022 relative to 2020, albeit with seasonal variability. Tau values for the period 2010–2022 based on surface measurements at Children's Park were −0.74 for both fall/winter and spring, −0.56 for dry summer, and −0.59 for monsoon summer, and Craycroft exhibited even stronger trends, with tau values of −0.82, −0.82, −0.79, and −0.69 for those same seasons in order (all p < 0.01). These trends are corroborated by NEI emissions data in Fig. S7, which show similar reductions in NO_x across major source sectors. Together, satellite, surface, and NEI emissions data reinforce the broader pattern of sustained NO₂ reductions that likely offset much of the urbanization-related pressure on O₃ levels and possibly help explain the lack of long-term change in seasonal medians of O₃.

Fig. 5 Top gray row of each panel shows the seasonally-resolved climatological median value for Ozone Monitoring Instrument (OMI) data for (a) NO₂, (b) HCHO, and (c) FNR for 2005–2020 over Tucson. The bottom three rows of each panel represent color-coded seasonal anomalies (described in Section 2.3) of OMI data.

Fig. 6 Same as Fig. 5 except for the period between 2020 and 2022 using data from the Tropospheric Monitoring Instrument (TROPOMI) for (a) NO₂, (b) HCHO, and (c) FNR.

HCHO anomalies did not show clear systematic changes between time intervals or within specific seasons of given years. This is supported by the lack of statistically significant trends across seasons, with Kendall's tau values for fall/winter, spring, dry summer, and monsoon summer being 0.27, 0.25, 0.02, and 0.33, respectively, for OMI. However, it should be noted that in 2020 (Fig. 5b) there were generally higher and positive HCHO anomalies relative to previous years, with the most pronounced anomalies for fall/winter and spring. TROPOMI HCHO anomalies (Fig. 6b) show that 2020 levels were elevated relative to 2021 and 2022, particularly during the spring and summer months. One possible contributing factor to the higher 2020 HCHO levels is regional fire activity since other studies^9,61 reported hotspots related to fire activities in northeastern areas of Phoenix and Tucson during 2020. Smoke transport from these fires may have contributed to elevated amounts of VOCs capable of enhancing HCHO levels.

NEI VOC estimates are bottom-up totals for all emitted species, including long-lived compounds that contribute little to formaldehyde, whereas satellite HCHO columns reflect the oxidation of short-lived, highly reactive VOCs (e.g., isoprene). As a secondary product, HCHO is also influenced by meteorology, transport, and NO_x levels. Furthermore, NEI values are reported every three years as annual totals, while satellite observations are daily and capture episodic events (e.g., drought, wildfire smoke) that inventories do not. These differences explain why trends in the satellite and NEI data may not align. The FNR anomalies from OMI (Fig. 5c) show a general increase from 2001–2009 to 2020, particularly during fall/winter and spring. The long-term increase in FNR from the 2000s to the 2020s is consistent with the results of Simon et al.² in their analysis of 51 U.S. non-attainment areas between 2002 and 2019. TROPOMI results from 2020 to 2022 (Fig. 6c) show that 2020 had elevated FNR relative to 2021 and 2022, as reflected by predominantly negative anomalies for 2021 and 2022 across all seasons. This pattern suggests that 2020 was an anomalous year, potentially influenced by pandemic-related emission reductions and enhanced VOC contributions from regional fire activity.

The FNR results explain how most sites exhibited a clearer weekend effect in 2005–2009, especially in fall/winter, with a reversal after that to a weekday effect until 2021 when some sites reverted again to a weekend effect (Fig. 3 and 4). This is because as FNR increases and moves further away from VOC-limited conditions, O₃ no longer increases when NO_x levels decrease, which was what drove higher weekend O₃ values previously due to less anthropogenic NO_x emissions on weekend days. Also, while MDA8 O₃ has not exhibited any significant long-term trend for any season across the entire study period, it is noteworthy that the levels have been mostly steadily rising from 2019 to 2022 (Fig. 2). Fig. 4's intriguing result that in spring the sites mostly all reverted back to a weekend effect in 2021–2022 can possibly now be explained by the TROPOMI results showing that FNR values dropped in those two years compared to previous years with the exception of 2005–2009 when it was low also consistent with only other period showing the weekend effect for spring. The post-2020 FNR decline represents a decrease from 2020, consistent with a positive anomaly in NO₂ and negative HCHO anomaly, rather than a return to pre-2020 levels or reversal of the long-term shift.

4. Discussion

4.1. Land type influence

Over the study period (2001–2022), changes in Tucson's land cover and emissions corroborate key findings regarding MDA8 O₃ trends and spatial variability. The land cover classification maps (Fig. S2) representing the period between 2001 and 2022, illustrate the extent of urban expansion with large areas of herbaceous and shrubland converted to medium and high intensity development. These spatial changes are consistent with the rise in population within the Tucson metropolitan statistical area, with an increase of 27% from 828 905 in 1999 to over 1.06 million in 2023,⁵ which led to increases in anthropogenic activity, including vehicular traffic and energy use. This aligns with the relatively elevated O₃ levels recorded at S/SE sites such as Saguaro National Park and Fairgrounds, where downwind transport from expanding urban zones likely contributed to increased precursor availability.

In addition, NDVI data can be used in a quantitative analysis of the relationship between land cover, emissions, and O₃. De-seasonalized monthly NDVI anomalies between 2001 and 2022 exhibited a substantial negative correlation (Tables S1 and S2) with O₃ anomalies (Pearson: r = −0.22, p < 0.01; Spearman: r = −0.28, p < 0.001), suggesting that increased vegetative cover is associated with reduced O₃ levels. Similar negative correlations between NDVI and O₃ have been documented in other studies. Dong et al.⁶² identified a weaker negative correlation (r = −0.08) between NDVI and O₃ in a large urban greenspace study in China, associating enhanced vegetation with decreased O₃ exposure in children. Similarly, Liu et al.⁶³ noted that elevated NDVI along urban roads correlated with reduced concentrations of pollutants such as NO₂, PM_2.5, and PM₁₀, suggesting that vegetation may mitigate air pollution via deposition. The relationships between NDVI with NO₂ and HCHO were weak and statistically insignificant. NDVI exhibited a positive correlation with relative humidity and a negative correlation with temperature, suggesting that climatic variability influences the observed relationship. Collectively, these findings support the hypothesis that urban growth and related vegetation changes might affect O₃ variability, although, the weak correlation strength suggests a complicated interaction among emissions, chemistry, and meteorology. Furthermore, the usage of correlations is limited in that it cannot explain causality and thus more research into land use change impacts on O₃ trends is warranted.

4.2. De-weathering analysis

The differing behaviors of ODVs and seasonal medians capture different parts of the O₃ long-term variability. While O₃ medians exhibit little statistically significant change (Table 2), ODVs showed a general decline until 2017 followed by an increase (Table 1). This difference (stable O₃ medians and changing ODVs) suggests that ODVs are sensitive to high-end events (e.g., the fourth-highest daily MDA8 O₃ each year), making them sensitive to meteorological variability (e.g., warmer and drier season, stagnation, wildfires) and emission changes, whereas medians, which are potentially more affected by background O₃, regional transport, and evolving VOC-NO_x chemical regime behavior, showed minimal long-term change.

The lack of a pronounced long-term trend in O₃ (2001–2022) is not only attributed to changes in emissions and chemical regime shifts, but interannual meteorological variability also likely plays a major role. Meteorology normalized analysis using GAM for data between 2010 and 2022 shows that most of the relatively small increasing trend in Tucson's O_3, though not significant, is largely controlled by year-to-year meteorological variability. The observed trend of 0.21 ppb per year is reduced to 0.04 ppb per year after normalization, suggesting that 79% of the increase (0.17 ppb per year) is meteorologically driven. The mean meteorological offset over this period is −0.41 ppb relative to the fixed climatology. This result reconciles the stability of seasonal medians of O₃ with ODVs and highlights the importance of interpreting long-term O₃ trends in the context of both emissions and meteorology.

5. Conclusion

This study comprehensively analyzes seasonal O₃ trends in Tucson, Arizona, over a 22 year period, revealing important spatial and temporal variations reflecting evolving emissions patterns, urban expansion, and population growth. Spring and summer exhibit the highest O₃ levels, particularly during the dry and monsoon summers, when photochemical activity and biogenic VOC emissions peak. Although no statistically significant long-term trend in O₃ levels was detected, differences across seasons and sites highlight important spatiotemporal factors. Two of the three S/SE locations (especially Saguaro National Park) consistently recorded the highest seasonal O₃ concentrations.

Despite sustained reductions in NO₂ emissions, continued land development and changes in HCHO point to a more complex, nonlinear relationship between emissions and O₃. O₃ levels have remained relatively stable, suggesting that shifting chemical regimes and seasonal meteorology continue to promote high O₃ levels in the region. Meteorology-normalized results further suggest that much of the interannual variability is weather-driven, with only a small change attributable to non-meteorological factors. The weekend effect has weakened over time as weekday O₃ levels have increased relative to weekend days, including having a weekday effect evident earlier in the summer (June) as compared to conditions in the late 1990s.¹¹

The integration of land cover change, meteorology, and emissions data reinforces that O₃ variability cannot be understood without accounting for the broader urban footprint and its evolving emissions landscape. This underscores the need to pair air quality monitoring with land use planning, population trends, and high-resolution satellite data. Our analysis is limited by the lack of widespread surface-level VOC measurements across Tucson's monitoring sites. Expanding surface monitoring for the combination of O₃, NO₂, and VOCs at multiple sites would improve the accuracy of chemical regime classification (e.g., VOC-limited, NO_x-limited, transition) and help refine our understanding of what drives O₃ levels. Such surface in situ measurements are necessary as there are limitations in using remote sensing data from the surface and space to determine the chemical regime at the surface. Given the complexity of these interactions, adaptive air quality management strategies that consider both anthropogenic and natural emission sources will be essential for effectively mitigating O₃ pollution in the region.

Author contributions

Taiwo Ajayi: writing – review & editing, writing – original draft, methodology, investigation, formal analysis, conceptualization. Mohammad Amin Mirrezaei: writing – review & editing, investigation, data curation. Ellis Robinson: review & editing, writing – original draft, methodology, investigation, conceptualization. Avelino F. Arellano: writing – review & editing, validation, supervision, project administration, methodology, investigation. Armin Sorooshian: writing – review & editing, writing – original draft, visualization, validation, supervision, resources, project administration, methodology, investigation, formal analysis.

Conflicts of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability

Daily MDA O₃ and NO₂ data and meteorological data from the U.S. Environmental Protection Agency Air Quality System are available at https://aqs.epa.gov/aqsweb/airdata/download_files.html (last access: 05 January 2025). Emission inventories data from the U.S. Environmental Protection Agency National Emission Inventory are available at https://www.epa.gov/air-emissions-inventories/get-air-emissions-data-0 (last access: 20 May 2025). Planetary boundary layer height (PBLH) data are available from the Modern Era-Retrospective Analysis for Research and Applications, version 2 (MERRA-2) at https://disc.gsfc.nasa.gov/datasets?project=MERRA-2/ (last access: 31 March 2025). Precipitation data are available from the Precipitation Estimation from Remotely Sensed Information using Artificial Neural Networks (PERSIANN) at https://chrsdata.eng.uci.edu/ (last access: 31 March 2025). OMI HCHO and NO₂ L2 data from the NASA Goddard Earth Sciences Data and Information Services Center (GES DISC) can be found at https://disc.gsfc.nasa.gov/ (last access: 18 March 2025). Population totals from the U.S. Census Bureau are available at https://data.census.gov/ (last access: 16 January 2025). TROPOMI L2 data are available from NO₂: https://disc.gsfc.nasa.gov/datasets/S5P_L2__NO2____1/summary (last access: 16 January 2025) and HCHO: https://disc.gsfc.nasa.gov/datasets/S5P_L2__HCHO___HiR_1/summary (last access: 16 January 2025). Land cover data from the National Land Cover Database (NLCD) are available at https://www.mrlc.gov/viewer/ (last access: 20 May 2025). Normalized Difference Vegetation Index (NDVI) data from Moderate Resolution Imaging Spectroradiometer (MODIS) using NASA's Application for Extracting and Exploring Analysis Ready Samples (AppEEARS) are available at https://appeears.earthdatacloud.nasa.gov/task/area.

Supplementary information: figures illustrating meteorological variables, land cover changes, spatial and temporal variability in O₃ and NO₂, emission trends of O₃ precursors, and tables of NDVI–air quality correlations. See DOI: https://doi.org/10.1039/d5ea00072f.

Acknowledgements

Funding was provided by the Arizona Board of Regents (ABOR) Regent's Grant from the Technology and Research Initiative Fund (TRIF).

References

1. H. Simon, A. Reff, B. Wells, J. Xing and N. Frank, Environ. Sci. Technol., 2015, 49, 186–195.
2. H. Simon, C. Hogrefe, A. Whitehill, K. M. Foley, J. Liljegren, N. Possiel, B. Wells, B. H. Henderson, L. C. Valin, G. Tonnesen, K. W. Appel and S. Koplitz, Atmos. Chem. Phys., 2024, 24, 1855–1871.
3. G. Prabhakar, E. A. Betterton, W. Conant and B. M. Herman, DOI: DOI:10.1175/JAMC-D-14-0001.1.
4. A. Sorooshian, A. F. Arellano, M. P. Fraser, P. Herckes, G. Betito, E. A. Betterton, R. A. Braun, Y. Guo, M. A. Mirrezaei and C. Roychoudhury, ACS EST Air, 2024, 1, 62–79.
5. Census Bureau Data, https://data.census.gov/, accessed 16 January 2025.
6. C. M. Atkinson-Palombo, J. A. Miller and R. C. Balling, Atmos. Environ., 2006, 40, 7644–7658.
7. M. Greenslade, Y. Guo, G. Betito, M. A. Mirrezaei, C. Roychoudhury, A. F. Arellano and A. Sorooshian, Atmos. Environ., 2024, 334, 120703.
8. J. M. Heuss, D. F. Kahlbaum and G. T. Wolff, J. Air Waste Manage. Assoc., 2003, 53, 772–788.
9. J. A. Miech, P. Herckes, M. P. Fraser, A. F. Arellano, M. A. Mirrezaei and Y. Guo, Atmosphere, 2024, 15, 555.
10. J. E. Diem and A. C. Comrie, Proceedings of the 1998 91st Annual Meeting & Exposition of the Air & Waste Management Association, 1998, p. 6.
11. J. E. Diem, Atmos. Environ., 2000, 34, 3445–3451.
12. X. Jin, A. Fiore, K. F. Boersma, I. D. Smedt and L. Valin, Environ. Sci. Technol., 2020, 54, 6518–6529.
13. J. E. Diem and A. C. Comrie, J. Environ. Manage., 2001, 63, 425–447.
14. J. S. Abraham and A. C. Comrie, J. Air Waste Manage. Assoc., 2004, 54, 914–925.
15. E. K. Wise and A. C. Comrie, J. Air Waste Manage. Assoc., 2005, 55, 1208–1216.
16. E. K. Wise, Int. J. Climatol., 2009, 29, 87–97.
17. H. Akimoto and H. Tanimoto, Atmos. Environ., 2022, 277, 119033.
18. O. R. Cooper, R.-S. Gao, D. Tarasick, T. Leblanc and C. Sweeney, J. Geophys. Res. Atmos., 2012, 117(D22) DOI:10.1029/2012JD018261.
19. O. R. Cooper, D. D. Parrish, J. Ziemke, N. V. Balashov, M. Cupeiro, I. E. Galbally, S. Gilge, L. Horowitz, N. R. Jensen, J.-F. Lamarque, V. Naik, S. J. Oltmans, J. Schwab, D. T. Shindell, A. M. Thompson, V. Thouret, Y. Wang and R. M. Zbinden, Elementa: Science of the Anthropocene, 2014, 2, 000029.
20. S. Mousavinezhad, M. Ghahremanloo, Y. Choi, A. Pouyaei, N. Khorshidian and B. Sadeghi, Atmos. Environ., 2023, 300, 119693.
21. D. D. Parrish, I. C. Faloona and R. G. Derwent, Atmos. Chem. Phys., 2025, 25, 263–289.
22. N. Lothrop, N. Lopez-Galvez, R. A. Canales, M. K. O'Rourke, S. Guerra and P. Beamer, Int. J. Environ. Res. Public Health, 2022, 19, 3173.
23. A. Sorooshian, A. Wonaschütz, E. G. Jarjour, B. I. Hashimoto, B. A. Schichtel and E. A. Betterton, J. Geophys. Res.:Atmos., 2011, 116, 16.
24. K. Jardine, L. Abrell, S. A. Kurc, T. Huxman, J. Ortega and A. Guenther, Atmos. Chem. Phys., 2010, 10, 12191–12206.
25. F. Dubertret, F.-M. Le Tourneau, M. L. Villarreal and L. M. Norman, Remote Sens., 2022, 14, 2127.
26. C. S. A. Wallace, M. L. Villarreal and L. M. Norman, Development of a High-Resolution Binational Vegetation Map of the Santa Cruz River Riparian Corridor and Surrounding Watershed, Southern Arizona and Northern Sonora, Mexico, U.S. Geological Survey, 2011.
27. B. Chu, S. Zhang, J. Liu, Q. Ma and H. He, J. Environ. Sci., 2021, 99, 346–353.
28. Y. Zhang, X. Wang, D. R. Blake, L. Li, Z. Zhang, S. Wang, H. Guo, F. S. C. Lee, B. Gao, L. Chan, D. Wu and F. S. Rowland, J. Geophys. Res. Atmos., 2012, 117(D15) DOI:10.1029/2011JD017356.
29. O. US EPA, Procedures for Detection and Quantitation – Documents, https://www.epa.gov/cwa-methods/procedures-detection-and-quantitation-documents, accessed 16 January 2025.
30. R. Gelaro, W. McCarty, M. J. Suárez, R. Todling, A. Molod, L. Takacs, C. A. Randles, A. Darmenov, M. G. Bosilovich, R. Reichle, K. Wargan, L. Coy, R. Cullather, C. Draper, S. Akella, V. Buchard, A. Conaty, A. M. da Silva, W. Gu, G.-K. Kim, R. Koster, R. Lucchesi, D. Merkova, J. E. Nielsen, G. Partyka, S. Pawson, W. Putman, M. Rienecker, S. D. Schubert, M. Sienkiewicz and B. Zhao, J. Clim., 2017, 30, 5419–5454.
31. P. Nguyen, E. J. Shearer, H. Tran, M. Ombadi, N. Hayatbini, T. Palacios, P. Huynh, D. Braithwaite, G. Updegraff, K. Hsu, B. Kuligowski, W. S. Logan and S. Sorooshian, Sci. Data, 2019, 6, 180296.
32. L. N. Lamsal, N. A. Krotkov, A. Vasilkov, S. Marchenko, W. Qin, E.-S. Yang, Z. Fasnacht, J. Joiner, S. Choi, D. Haffner, W. H. Swartz, B. Fisher and E. Bucsela, Atmos. Meas. Tech., 2021, 14, 455–479.
33. I. De Smedt, N. Theys, H. Yu, T. Danckaert, C. Lerot, S. Compernolle, M. Van Roozendael, A. Richter, A. Hilboll, E. Peters, M. Pedergnana, D. Loyola, S. Beirle, T. Wagner, H. Eskes, J. van Geffen, K. F. Boersma and P. Veefkind, Atmos. Meas. Tech., 2018, 11, 2395–2426.
34. A. Lorente, K. Folkert Boersma, H. Yu, S. Dörner, A. Hilboll, A. Richter, M. Liu, L. N. Lamsal, M. Barkley, I. De Smedt, M. Van Roozendael, Y. Wang, T. Wagner, S. Beirle, J.-T. Lin, N. Krotkov, P. Stammes, P. Wang, H. J. Eskes and M. Krol, Atmos. Meas. Tech., 2017, 10, 759–782.
35. J. van Geffen, K. F. Boersma, H. Eskes, M. Sneep, M. ter Linden, M. Zara and J. P. Veefkind, Atmos. Meas. Tech., 2020, 13, 1315–1335.
36. A. Ludewig, Q. Kleipool, R. Bartstra, R. Landzaat, J. Leloux, E. Loots, P. Meijering, E. van der Plas, N. Rozemeijer, F. Vonk and P. Veefkind, Atmos. Meas. Tech., 2020, 13, 3561–3580.
37. Sentinel-5P Nitrogen Dioxide Level 2 Product Readme File, https://copernicus.eu/user-guides/sentinel-5p-tropomi/document-library/-/asset_publisher/w9Mnd6VPjXlc/content/sentinel-5p-nitrogen-dioxide-level-2-product-readme-file, accessed 18 October 2024.
38. M. A. Mirrezaei, A. Arellano, Y. Guo, C. Roychoudhury and A. Sorooshian, Environ. Res. Commun., 2024, 6, 051009.
39. S. Koplitz, H. Simon, B. Henderson, J. Liljegren, G. Tonnesen, A. Whitehill and B. Wells, ACS Environ. Au, 2022, 2, 206–222.
40. J.-S. Youn, Z. Wang, A. Wonaschütz, A. Arellano, E. A. Betterton and A. Sorooshian, Geophys. Res. Lett., 2013, 40, 3468–3472.
41. E. Crosbie, J.-S. Youn, B. Balch, A. Wonaschütz, T. Shingler, Z. Wang, W. C. Conant, E. A. Betterton and A. Sorooshian, Atmos. Chem. Phys., 2015, 15, 6943–6958.
42. A. Raman, A. F. Arellano and A. Sorooshian, Atmosphere, 2016, 7, 24.
43. T. Hastie and R. Tibshirani, Stat. Sci., 1986, 1, 297–310.
44. S. N. Wood, Generalized Additive Models: an Introduction with R, Second Edition, Chapman and Hall/CRC, New York, 2nd edn, 2017.
45. Z. Gao, C. E. Ivey, C. L. Blanchard, K. Do, S.-M. Lee and A. G. Russell, Environ. Pollut., 2022, 307, 119503.
46. H. Lee and D. A. Jaffe, J. Air Waste Manage. Assoc., 2024, 74, 116–130.
47. Y. Ma, B. Ma, H. Jiao, Y. Zhang, J. Xin and Z. Yu, Atmos. Environ., 2020, 224, 117342.
48. S. Munir, G. Coskuner, M. S. Jassim, Y. A. Aina, A. Ali and M. Mayfield, Atmosphere, 2021, 12, 504.
49. A. Chauhan, G. P. Sai and C.-Y. Hsu, Heliyon, 2025, 11, e41762.
50. H. Lee and D. A. Jaffe, Environ. Sci. Technol., 2024, 58, 14764–14774.
51. H. Luo and C.-H. Lu, J. Geophys. Res. Atmos., 2024, 129, e2023JD039206.
52. C. L. Blanchard, S. Tanenbaum and D. R. Lawson, J. Air Waste Manage. Assoc., 2008, 58, 1598–1615.
53. L. C. Marr and R. A. Harley, Environ. Sci. Technol., 2002, 36, 4099–4106.
54. L. C. Marr and R. A. Harley, Atmos. Environ., 2002, 36, 2327–2335.
55. L. R. Chinkin, D. L. Coe, T. H. Funk, H. R. Hafner, P. T. Roberts, P. A. Ryan and D. R. Lawson, J. Air Waste Manage. Assoc., 2003, 53, 829–843.
56. E. M. Fujita, D. E. Campbell, B. Zielinska, J. C. Sagebiel, J. L. Bowen, W. S. Goliff, W. R. Stockwell and D. R. Lawson, J. Air Waste Manag. Assoc., 2003, 53, 844–863.
57. S. Singh and I. G. Kavouras, Atmosphere, 2022, 13, 114.
58. D. Jaffe and J. Ray, Atmos. Environ., 2007, 41, 5452–5463.
59. C. E. Buysse, J. A. Munyan, C. A. Bailey, A. Kotsakis, J. A. Sagona, A. Esperanza and S. E. Pusede, Atmos. Chem. Phys., 2018, 18, 17061–17076.
60. S. Baidar, R. M. Hardesty, S.-W. Kim, A. O. Langford, H. Oetjen, C. J. Senff, M. Trainer and R. Volkamer, Geophys. Res. Lett., 2015, 42, 9457–9464.
61. Y. Guo, C. Roychoudhury, M. A. Mirrezaei, R. Kumar, A. Sorooshian and A. F. Arellano, Geosci. Model Dev., 2024, 17, 4331–4353.
62. P.-P. Dong, G. Chen, Y. Wen, J. Heinrich, E. Fuertes, T. Zhao, N. S. Idrose, L.-Z. Lin, Z.-H. Gui, X. Qin, H. Tu and G.-H. Dong, Sustain. Cities Soc., 2025, 124, 106329.
63. C. Liu, A. Dai, Q. Sheng and Z. Zhu, Atmos. Pollut. Res., 2024, 15, 102003.

---