Air pollution and meteorology as risk factors for COVID-19 death in a cohort from Southern California

Background Recent evidence links ambient air pollution to COVID-19 incidence, severity, and death, but few studies have analyzed individual-level mortality data with high quality exposure models. Methods We sought to assess whether higher air pollution exposures led to greater risk of death during or after hospitalization in confirmed COVID-19 cases among patients who were members of the Kaiser Permanente Southern California (KPSC) healthcare system (N=21,415 between 06-01-2020 and 01-31-2022 of whom 99.85 % were unvaccinated during the study period). We used 1 km resolution chemical transport models to estimate ambient concentrations of several common air pollutants, including ozone, nitrogen dioxide, and fine particle matter (PM2.5). We also derived estimates of pollutant exposures from ultra-fine particulate matter (PM0.1), PM chemical species, and PM sources. We employed Cox proportional hazards models to assess associations between air pollution exposures and death from COVID-19 among hospitalized patients. Findings We found significant associations between COVID-19 death and several air pollution exposures, including: PM2.5 mass, PM0.1 mass, PM2.5 nitrates, PM2.5 elemental carbon, PM2.5 on-road diesel, and PM2.5 on-road gasoline. Based on the interquartile (IQR) exposure increment, effect sizes ranged from hazard ratios (HR) = 1.12 for PM2.5 mass and PM2.5 nitrate to HR ∼ 1.06–1.07 for other species or source markers. Humidity and temperature in the month of diagnosis were also significant negative predictors of COVID-19 death and negative modifiers of the air pollution effects. Interpretation Air pollution exposures and meteorology were associated the risk of COVID-19 death in a cohort of patients from Southern California. These findings have implications for prevention of death from COVID-19 and for future pandemics.


Introduction
The COVID-19 pandemic represents one of the largest threats to population health in more than a century. Currently, more than 624 million people worldwide have been diagnosed with COVID-19, resulting in more than 6.5 million deaths (World Health Organization, 2022). Researchers have extensively investigated the etiology of COVID-19, yet considerable uncertainties remain on how risk factors may influence COVID-19 incidence, severity, and death. Recent evidence from the North America, Asia, and Europe implicates air pollution as a risk factor for COVID-19 incidence, prognosis, and death (Brandt et al., 2020;Li et al., 2020;Wu et al., 2020;Zhang et al., 2020;Zhu et al., 2020;Lippi et al., 2019;Coker et al., 2020;Travaglio et al., 2021;Yao et al., 2021;Huang et al., 2021;Chen et al., 2021;Berg et al., 2021;Zhou et al., 2021).
In reviewing the growing literature on air pollution exposure and COVID-19 outcomes, we found only four other mortality studies have used individual level data with some levels of control for potential confounders (Chen et al., 2021;Bozack et al., 2021;Elliott et al., 2021;Nobile et al., 2022;Chen et al., 2022). These studies focused on the early phases of the pandemic, which may have led to lower statistical power due to a relatively small number of deaths. While some of the mortality studies used high-quality exposure estimates, none assessed source contributions or ultrafine particle concentrations. In addition, none of these studies examined interactions between air pollution and meteorological variables such as temperature and humidity. Here we expand the evidence base with a large sample of individual data, a longer study period, exposure models capable of assessing particle species and sources, and meteorological variables. In this context, we addressed two research objectives. Firstly, we assessed whether higher air pollution exposures led to greater risk of death in confirmed COVID-19 cases among patients who were members of the Kaiser Permanente Southern California (KPSC) healthcare system. Secondly, we investigated whether meteorology variables influenced the risk of COVID-19 death or modified associations between air pollution and COVID-19 death.

KPSC cohort and health data
KPSC is a large integrated health care system with a racially/ethnically and socioeconomically diverse membership of 4.7 million people, living across nine southern California counties. KPSC's membership approximately represents the underlying population of the second largest urban region in the United States; further details of the KPSC membership are described elsewhere (Koebnick et al., 2012). KPSC's Electronic Health Record (EHR) is an integrated data system that captures all aspects of patient care, including diagnoses, inpatient and outpatient encounters, pharmacy encounters, and laboratory tests.
Clinical care changed rapidly during the first months of the pandemic. We therefore began our observation period on 06/01/2020 when new standards of COVID care, such as lying patients in the prone position, had become more common. We identified patients with KPSC COVID-19 molecular diagnostic tests and diagnoses (ICD-10 codes: J12.89, J20.8, J22, J80, B34.2, B97.29, U07.1) from 06/01/2020 to 01/ 30/2021. We include both diagnoses and COVID-19 tests because patients may have been tested outside of KPSC and received a diagnosis at KPSC without being re-tested.
The study cohort is comprised of patients who were 18 years or older at the time of diagnoses or positive COVID-19 test. We limited our sample to members who had at least 1 year of membership before their COVID-19 diagnoses/test to reliably assess co-morbidities. We defined COVID-19 hospitalizations as hospitalizations occurring within 21 days of COVID-19 diagnoses or positive test (N = 316,224) (Nau et al., 2021). We used hospitalized patients from the cohort rather from all those who tested positive because testing could have occurred after possible contact with an infected person or after the onset of severe illness at the point of hospital admission. This would result in uncertainty about the time window at which the test could have occurred among different study patients that would introduce substantial errors in our follow up times, which would lead to biased results in the statistical models. Restricting to those hospitalized eliminated this potential problem, as no uncertainty existed in the time of hospitalization. After applying eligibility and exclusion criteria, the analytic cohort consisted of 21,415. Deaths were included up to 90 days after the initial hospitalization (see Online Data Supplement [ODS] for further details on death ascertainment). We excluded patients who lost membership during our 90-observation window and who were hospitalized for childbirth. This study was approved by the Kaiser Permanente Institutional Review Board.
The KPSC EHR provides information on patient age and sex. Member race/ethnicity categories have been created using a validated algorithm that uses multiple data sources (Nau et al., 2021).
Body mass index (BMI: kg/m 2 ) has been found to be an important risk factor for COVID-19 mortality (Tartof et al., 2020). The most recent BMI available in the EHR was used to adjust for this potential confounder (Tartof et al., 2020). We cleaned BMI data using validated algorithms to delete biologically implausible values.
Five broad comorbidity categories that have been used in prior COVID-19 research were created to identify co-morbidities that may increase a person's risk of severe COVID-19 outcomes (Nau et al., 2021;Quan et al., 2005). We use Elixhauser disease categories to create COVID-19 relevant disease categories (see ODS for further details).
Smoking status and the Exercise Vital Sign (EVS) data are collected during each KPSC in-person outpatient health care encounter. Smoking status was coded (ever or never) based on the information provided during the last encounter before the COVID-19 test/diagnoses reaching back up to four years. The EVS queries on usual exercise is coded in the EHR in minutes/week of moderate to vigorous exercise. All EVS information for the past four years was identified for every patient. The median value of minutes of exercise per week was calculated and used in our analysis (Young et al., 2018;Zhou et al., 2021).
We identified patients who were enrolled at KPSC via MediCal to identify patients with very low income. In sum, four individual-level confounders were considered: smoking status, BMI, Medicaid (low income), and EVS.
We queried vaccination status and found only 33 members of our cohort were vaccinated prior to hospitalization; thus about 99.85 % of the cohort was unvaccinated during the study period.
We also followed common practice in analyses of EHR data and added predictors of community-level SES to help proxy individual SES and adjust for community level effects of social determinants of health (Krieger, 1992;Diez-Roux et al., 2001;Geronimus and Bound, 1998). Community-level predictors at the census block-group level were drawn from the American US Census Bureau, 2018; (Messer et al., 2006). They include a validated neighborhood deprivation index (NDI), a measure of crowding (the proportion of households with more than one occupant per home), and the proportion of workers aged 16 and older who commute to work via public transportation (Messer et al., 2006).
GridMET are high-spatial resolution (~4-km) surface meteorological data covering the contiguous U.S. We acquired GridMET daily maximum temperature and daily maximum relative humidity for our entire study period through Google Earth Engine (https://developers.google.com/ea rth-engine/datasets/catalog/IDAHO_EPSCOR_GRIDMET?hl = en). We aggregated the GridMET data to monthly means for the home address of every study participant up to the month where they were hospitalized with COVID-19.

Exposure Assessment: Chemical transport model
Exposure simulations were carried out across California using the source-oriented UC Davis-California Institute of Technology (UCD-CIT) 3D reactive chemical transport model (CTM)  .The UCD/CIT model predicts the evolution of gas and particle phase pollutants in the atmosphere in the presence of emissions, transport, deposition, chemical reaction, and phase change. The pressing timeline for the current study during an ongoing public health crisis necessitated leveraging past efforts that prepared and validated CTM inputs. We previously reported CTM exposure fields with 4 km resolution over California for the years 2000-2016 . The most recent year (i.e., 2016) of this time window was selected as the starting point to characterize chronic exposure in the current study. Meteorology and emissions inputs for the year 2016 were downscaled to improve spatial resolution to 1 km. Bias in the raw CTM output fields was removed using a constrained regression model based on source apportionment tags and the difference between predicted and measured concentrations. See ODS for further details on CTM methods.
CTM predictions include a wide range of pollutants. For our study area, we estimated PM 2.5 total mass, PM 2.5 nitrates, PM 2.5 organic carbon (OC), PM 2.5 elemental carbon (EC), ultra-fine particle mass or PM 0.1 for particles with diameters of 100 nm or less, nitrogen dioxide (NO 2 ), and ozone (O 3 ). We also extracted source tracers for on-road diesel, onroad gasoline, and biomass burning. These exposure fields were assigned to the geocoded home address of the cohort members. Although the exposure fields were restricted to 2016, we accounted for population mobility by assigning exposures to each address for any member of the cohort who moved in the past five years. We then did time-weighted averaging of the exposures to account for mobility effects for those who had moved in the preceding 5 years.

Statistical models
We used Cox proportional hazards models with adjustment for potential individual and neighborhood confounders. All models were stratified at baseline for age, sex, and race-ethnicity. Age was included in 5-year intervals. We controlled for potential non-independence at the census-tract level with a sandwich estimator, which allowed for robust variance estimation. All analyses were run in the R package version 4. 0.4 (2021-02-15).
The Cox model estimates the instantaneous hazard of dying during the follow up as: where, h ij (t): hazard function for the ith subject in jth census tract neighborhood; h 0s (t): the baseline hazard function for stratum s (i.e., age, race and sex); P ij : air pollution exposure metric of interest (e.g., PM 2.5 ) standardized to the interquartile range for individual i in census tract j; X ij : individual risk factors (i.e., smoking, exercise, BMI, poverty) for individual i in census tract j; Z ij : neighborhood risk factors (i.e., deprivation index, proportion taking public transit, crowding) for individual i in census tract j; and.
W itj : weather conditions (i.e., maximum temperature and humidity) for individual i at the t th month of admission in census tract j.
Equation (1) above represents the general form of the model. Confounders, however, were selected for each pollutant with the following procedure: We ran unadjusted models stratified by age, race/ethnicity, and sex for each pollutant exposures. We tested every possible confounder (BMI, smoking, etc) one at a time with each pollution estimate. We included any confounder that changed the unadjusted pollution coefficient by at least 10 %. We subsequently ran the adjusted models for all pollution exposures that included variables meeting the 10 % criterion. Exposures were standardized for comparison across pollutants by dividing each by their respective interquartile ranges (IQR).
For pollutants that had statistically significant effects at conventional levels after adjustment (i.e., p < 0.05), we then conducted stratified analyses on variables that could modify the association between COVID-19 death and air pollution, including race/ethnicity, sex, age, and number of chronic diseases categories.
We also tested for interaction by running models with a multiplicative term with one pollutant and one meteorological variable. When significant interactions were present based on the p-value of the interaction term, we stratified the HR estimates for the pollutant by tertile of the meteorological variable.
We examined two-pollutant models (i.e., O 3 and NO 2 , NO 2 and PM 2.5 mass, and O 3 and PM 2.5 mass). We also explored the concentration-response (CR) functions for each pollutant that had a significant individual effect in a fully adjusted model. The CR functions were estimated via the pspline function in the gam library in R.
We also investigated the potential contribution of different SARS-COV2 variants by performing sensitivity analyses that were restricted to periods when the Delta variant was dominant. The California Department of Public Health has done retrospective genomic analyses on specimens from all stages of the pandemic (https://covid19.ca.gov/ variants/). In the early part of our study, five different variants were circulating. For much of the study period, the Delta variant was dominant. In the last month or so of our study, the Omicron variant became dominant. It is likely, however, that many of the hospitalizations that would have occurred in the last weeks to month of our study would have resulted from Delta due to latency in the infection time and the time required for a person to become ill enough to be hospitalized.
Our sensitivity analyses focused on the period of 06/19/2020 to 01/ 03/2021. The first date corresponds to the initial timepoint at which the Delta variant accounted for more than 50 % of the cases. We then identified the point at which Delta lost dominance (i.e., greater than 50 %) as 12/19/2020. We added a two week buffer to this end date on the assumption that many of the hospitalizations and subsequent deaths from Delta would have taken at least two weeks to occur. This yields a conservative estimate for the end of the sensitivity analysis to be 01/03/ 2021. We reran our analyses for PM 2.5 during this restricted time period so that results could be compared to the main analysis.

Role of the funder
The California Air Resources Board (CARB) funded most of this research and oversaw peer review of an unpublished final report documenting the methods and results, which was required by the terms of the funded contract. CARB staff also offered comments intended to improve the clarity of presentation in the final report. The Health Effects Institute also funded some of the study, but had no active role in the research. Neither funder had any role in the decision to publish this manuscript. Table 1 displays the descriptive statistics for the cohort of the 21,415 KPSC patients who were hospitalized with COVID-19, of whom 4,815 died. Cohort characteristics of were as follows: median age 64 (IQR: 52, 75), 58 % male, 56 % Hispanic origin, 23 % white, 11 % Asian/Pacific Islanders, 8.6 % Black, and 1.6 % were of other or unknown race/ ethnicity. Some 37 % were ever smokers, and 13 % had health insurance through MediCal, a government health program for low-income people. Fig. 1 shows the hospitalizations over the entire time period, which included a major surge in admissions in November 2020 to the end of our study period.

Descriptive statistics
Most hospitalized patients were overweight or obese, with 29 % meeting criteria for overweight, 27 % for obesity class 1, and 30 % for obesity class 2 or higher. Some 41 % had a history of cardiovascular disease, 59 % had hypertension, 22 % had chronic obstructive pulmonary disease, 45 % had diabetes, and 65 % had another chronic condition.
Patients who died within 90 days of their first hospitalization were older than those who did not (median age 74 vs 61 years), more likely to be male (63 % vs 56 %), and more likely to be ever smokers (46 % vs 34 %). Patients who died had more comorbidities (median Elixhauser index 5.0 vs 2.0) and greater prevalence of chronic diseases. Of those who died, 63 % had prevalent cardiovascular disease (vs 35 % in survivors), 76 % had hypertension (vs 54 %), 55 % had diabetes (vs 42 %), and 26 % had chronic obstructive pulmonary disease (COPD) (vs 20 %).
Descriptive statistics for the pollutants are shown in Table 2. Table 3 shows many of the pollutants had moderate to high correlations with one another (e.g., PM 2.5 and PM 2.5 nitrate r ~ 0.9). Ozone was the least correlated with the other pollutants and, as expected, had negative associations with NO 2 (r ~ 0.25) and some of the particle species or source tracers. Fig. 2 shows spatial distribution of several pollutants across Southern California, including: PM 2.5 mass, PM 2.5 nitrates, PM 2.5 EC, PM 0.1 as well as on-road gasoline and diesel tracers. Substantial differences exist in the spatial patterns among several pollutants. For example, on-road gasoline displayed variation consistent with highways that carry large volumes of traffic, while PM 2.5 and PM 2.5 nitrate had more smoothly-varying exposures across the region, likely due to a large contributions to the mass from secondary formation in the atmosphere. All pollutants had relatively higher concentrations in the inland areas of San Bernardino and Riverside.

Results from adjusted models
The confounders selected for each pollutant in our adjusted models are shown Table 3 of the ODS. Fig. 3 Table 5 show the main results on associations between air pollution and COVID-19 death. After confounding adjustment, we found several air pollutants were related to COVID-19 death among hospitalized patients including: PM 2.5 mass (HR = 1.12, 95 % CI 1.06, 1.17); PM 2.5 nitrates (HR = 1.12, 95 % CI 1.07, 1.17); PM 2.5 EC (HR = 1.07, 95 % CI 1.03, 1.12); PM 0.1 mass (HR = 1.06, 95 % CI 1.02, 1.10); PM 2.5 on-road diesel (HR = 1.06, 95 % CI 1.03, 1.10); and PM 2.5 on-road gasoline (HR = 1.07, 95 % CI 1.02, 1.13). Effects of PM 2.5 mass were partly confounded by NO 2 in the two pollutant models, but remained significantly elevated (Fig. 3). For the Delta variant-dominant period, results were similar to those from the main model for PM 2.5 mass (HR = 1.13, 95 % CI 1.07, 1.20). Effects for gaseous species were sensitive to co-pollutant adjustment. NO 2 had a significant association with the risk of death (HR = 1.10, 95 % CI 1.04, 1.16), while ozone had positive but insignificant effects (HR = 1.02, 95 % CI 0.96, 1.08). Because the inverse spatial pattern that can lead to positive confounding (Quiros et al., 2013), we also ran copollutant models with ozone and NO 2 included. In these models, NO 2 remained significantly elevated, but ozone remained null. PM 2.5 confounded the NO 2 effect to null when both were included in the same model (Fig. 3).

Stratification analyses
All of the subgroup analyses were insignificant based on the Q statistic shown at the bottom of each table, meaning these variables had no significant impact on the air pollution concentration-response association with COVID-19 death (see Tables 6-9 for stratification analyses in the ODS).

Interaction models with meteorological variables
After determining that temperature and humidity significantly modified the effects of air pollution on COVID-19 death, we stratified by tertile for these variables to visualize the effect modification for PM 2.5 (see Fig. 4). See ODS Figure 13 for other pollutants. For most of the pollutants, elevated risks appear only in the lower two tertiles of temperature. Effect modification was particularly pronounced for humidity, with most pollutants showing a graded decline in effects as humidity went up. PM 0.1 mass and PM 2.5 nitrates followed a slightly different trend with the largest effect in the middle tertile. Overall, most effects were present only in the lowest two tertiles of humidity.

Concentration-Response analysis
Concentration-response curves are shown in the Fig. 5. For most of the pollutants, we observed fairly linear curves when sufficient data was available to support the spline derivation. Some pollutants such as EC and on-road diesel displayed a supra-linear response with a steeper response curve at the low exposure levels. This supra-linear function has been observed in many air pollution-mortality studies (Burnett et al., 2018). Humidity displayed a clear linear negative association with risk of COVID-19 death. Temperature had a U-shape risks appear to be higher at lower temperatures, although where there was sufficient data to support the spline derivation, the inverse curve appeared linear.

Discussion and conclusion
Here we evaluated whether chronic exposure to air pollution and meteorology at the time of diagnosis affected the risk of death in patients with a COVID-19-related hospitalization. We found significant associations between the risk of COVID-19 death following hospitalization and PM 2.5 mass, PM 0.1 mass, and several of the particle species or source tracers. Effects for PM 2.5 mass were reduced when NO 2 was included in the model, but remained significantly elevated, while NO 2 was confounded to null in the two-pollutant model.
Meteorology has been associated with COVID-19 transmission in some studies (Zoran et al., 2022), and recent studies show that meteorology likely affected COVID-19 death rates in Europe (Kifer et al., 2021). These researchers proposed that humidity can interfere with viral defenses of nasal mucosa tissues and with the sputum deeper in the airway, which can lead to more severe infection and subsequently contribute to a poor prognosis after the virus establishes itself in the respiratory tract, particularly in the nose (Kifer et al., 2021;Weaver et al., 2022). Temperature and humidity can also affect the size of the viral droplets and its persistence in ambient air, but the extent to which this would affect severity is unknown (Bourdrel et al., 2021). Significant negative effects were present for both temperature and humidity in our study. We also found significant effect modification of the air pollution associations with lower temperature and humidity being associated generally with larger air pollution effects. If the viral defenses are influenced by meteorology, both direct effects of humidity and temperature and the effect modification of the pollution effect have biological plausibility.
In comparing our results to other mortality studies, (Chen et al., 2021) investigated the association of air pollution on COVID-19 severity and mortality using data from KPSC members with a CALINE dispersion model, which estimated traffic exposures as NOx (non-freeway and freeway). The odds of intensive care admission were 1.11 (95 % CI: 1.04, 1.19) and death were 1.10 (95 % CI: 1.03, 1.18) for each SD increase in non-freeway NOx. Several other freeway exposures, however, had protective effects (Chen et al., 2021). Including regional PM 2.5 and NO 2 as confounders attenuated the effects by 19-26 %, and this adjustment caused the freeway NOx to become significantly protective for mortality (HR = 0.93, 95 % CI: 0.87-1.0). Possibly, exposure measurement error may have been present due to the inability of the CALINE dispersion model to deal with complex traffic, terrain, and meteorological conditions, all which exist in Southern California Dhyani et al., 2013). Our findings may have also differed due to the longer follow up in the present study (with about 4.5 times as many deaths).
A follow up to this Chen et al. (2021) study using the same health data, but relying on inverse-distance averaging to interpolate from government monitors, found significant chronic effects associated with PM 2.5 exposure and sub-chronic effects from NO 2 (Chen et al., 2022). This study, however, also had a high probability of exposure measurement error given the likely level of spatial variation in these pollutants and the sparse data support available from the government monitors of which there are relatively few covering thousands of square kilometers.
Another study from the UK relied on the Biobank data and used an agnostic exposomic statistical approach to evaluate many factors for risk of COVID-19 incidence and death (Elliott et al., 2021). Although mild associations were present in univariate models with PM 2.5 , these were eliminated in multivariate models, leading the authors to conclude there was little evidence of an independent association between air pollution and COVID-19 death. This study, however, had relatively few deaths and may have lacked power to detect subtle effects from air pollution.
A study using data from hospitalized patients in New York City reported an association between PM 2.5 and risk of mortality (risk ratio, 1.11, 95 % CI: 1.02-1.21) per 1 μg/m 3 increase) (Bozack et al., 2021). Evaluated across the reported IQR of 0.7 μg/m 3 the rate ratio would be ~ 1.08. Neither black carbon nor NO 2 had a significant association with COVID-19 death. This study also found Hispanic ethnicity significantly modified the air pollution risk for COVID-19 death, which differs from our finding of no significant subgroup interaction. This study lacked individual information on some potential risk factors for COVID-19 death, including obesity and smoking. Consequently, residual confounding cannot be ruled out.
In a large administrative cohort from Rome Italy, significant associations with COVID-19 mortality were found with both NO 2 and PM 2.5 (Nobile et al., 2022). Associations found in the Rome study were somewhat smaller than what we have found here, but the confidence intervals for the two studies overlap The range of exposure in Rome was much smaller than what we observed in Southern California, which may in part explain the smaller effects observed there.
On limitations, while we did control for several individual confounders such as smoking and obesity, the KPSC data did not include potentially important confounding variables such as occupational status. Nascent research suggests increased risk of mortality in some occupational groups in California, particularly in the farming, material moving, transportation, and construction sectors, all of which could have higher air pollution due to occupational exposures (Cummings et al., 2021). While numerous complexities exist in analyzing and interpreting COVID-19 mortality risk in different occupations (Cummings et al., 2021;Pearce et al., 2021), it is plausible that lack of occupational status could have biased our results.
In addition, a temporal mismatch existed between the exposure fields from 2016, which predated the study by some three years; however, overall spatial patterns of exposure are unlikely to change during this period. Some portions of our study, however, overlapped with the "lock down" period when traffic emissions were lower . Cohort members may therefore have experienced lower exposures than they would have if normal conditions had prevailed, which would not have been accounted for in our exposure or statistical models. The impact would have been to overestimate their exposures near-source traffic exposures, which may have biased some results toward null. The near-road pollutants such as PM 2.5 EC, PM 0.1 on-road diesel and onroad gasoline had risks that were smaller than regionally-varying PM 2.5 and PM 2.5 nitrate, which might have been due to the lock down effect not captured in our exposure model. Nevertheless, despite this limitation, several near-source pollutants still displayed significant associations with COVID-19 death. We were also unable to account for acute effects, which may have contributed to risks of COVID-19 death. Currently, we are extending the CTM exposure modeling to derive contemporaneous estimates of acute and chronic exposure.
Another concern with observational studies of COVID-19 and mortality rests in the different variants that emerged and gained dominance through the pandemic. If certain variants were more virulent than others  as appears likely (Adjei, 2022) and these emerge coincidentally at times when air pollution is high, then associations between air pollution and mortality could be confounded by the virulence of the variant. In this study, the Delta variant was dominant for the majority of the study period. We performed sensitivity tests on the PM 2.5 model by restricting the analysis period to times when Delta was dominant. The results from the restricted analysis were virtually the same as the results from the full analysis. Based on this similarity, we conclude that it was unlikely that our results are confounded by variants with different virulence.
We used time-to-event data after hospitalization to avoid having bias in our follow up times, which could have varied considerably if we began the study at the point of COVID-19 diagnosis. Although necessary for unbiased statistical inference, this restriction reduces the generalizability of the results to hospitalized individuals rather than the general population.
Other environmental variables have also been implicated in the spread and severity of COVID-19, including wind speed and ultraviolet radiation. Both variables were explicitly included in our CTM exposure model. Wind speed in particular has a major impact on ambient concentrations of several pollutants, and we were concerned that including wind speed as its own variable would induce collinearity into the model. In reviewing the literature on wind speed we also found that most of the influence of this meteorological parameter affected the spread of COVID-19, not the severity of symptoms or risk of death (Weaver et al., 2022;Bourdrel et al., 2021).
UVB potentially operates through vitamin D deficiency, which has been identified as a risk factor for more extreme COVID-19 outcomes . We visually explored the 1 km UVB fields used as inputs in the CTM modeling. UVB levels were higher inland and lower near the Pacific Coastline likely due to fog and cloud cover. Recent UVB exposure modeling, however, estimates that personal behavior and occupation are much more important predictors of UVB exposure than ambient levels alone, which often account for little of the explained variation in objectively measured UVB (Dadvand et al., 2011;Soueid et al., 2022). Thus, the ambient levels are unlikely to be reasonable proxies for exposure and subsequent deficiency.
We also queried our data base to identify patients who were vitamin D deficient and run stratification analyses to assess whether air pollution contributed to worse outcomes in these patients. Some 4,142 (19.64 %) of the hospitalized cohort had a vitamin D lab test within 1 year prior to COVID test date and out of that group, 1,524 (7.23 % of total cohort) were vitamin D deficient (25-HYDROXYVITAMIN D lab result < 30 ng/ mL) based on their most recent vitamin D lab prior to hospitalization. Because this is relatively small proportion of the cohort and likely represents an underestimate that may bias results, we were unable to stratify the analysis by vitamin D deficiency.
Virucidal activity also decreases in the presence of higher UVB radiation (Weaver et al., 2022;Bourdrel et al., 2021), but this would be more likely to affect the spread of the virus and not the severity. In addition, it is likely that most of the infections occurred from contact in the indoor environment where ambient UVB would likely have a minimal impact on the virucidal activity (Weaver et al., 2022). Future research is nevertheless needed to assess whether vitamin D deficiency modifies air pollution effects on COVID-19 severity.
Saturation of capacity of the ICU care is another possible factor affecting survival that may have acted as a confounder. Internal data and consultation with attending physicians indicated that despite this surge (see Fig. 1), at no time were the ICU units saturated beyond capacity. KPSC did not run out of ventilators or physical space for admitting seriously ill COVID-19 patients. An overflow facility that could have accepted KPSC patients was never used. Consequently, saturation of the ICU is unlikely to confound our results.
The observational nature of this study precludes causal interpretation. Based on our results, however, we can conclude that chronic exposure to air pollution and meteorology in the month of diagnosis in Southern California is associated the risk of death from COVID-19.
Better knowledge about environmental variables such as air pollution and meteorology could be used by communities and local governments to target neighborhoods with higher risks for COVID-19 death. Such information could also be brought into healthcare systems to assist clinicians with better estimating the likely severity of disease in patients from high air pollution areas. Minimizing the spread and reducing the severity of COVID-19 through non-pharmaceutical interventions (NPIs) such as masking and economic shutdowns remains problematic over the longer term (Imai et al., 2020) given the social and economic costs involved. In addition, modeling suggests that NPI measures have the potential to increase the severity of other respiratory viral outbreaks in the future (Baker et al., 2020). Pharmaceutical measures such as vaccines continue to have mixed results in part due to vaccine hesitancy in some high prevalence locations and population groups (Sallam, 2021). In contrast, air pollution is a modifiable environmental risk factor that could affect disease severity across the entire population. Reducing air pollution may thus provide a more sustainable means of reducing COVID-19 severity that would have substantial population benefits. It may also reduce the risks for catastrophic outcomes from future pandemics fueled by novel viruses, while also having beneficial effects on a wide array of other health endpoints.

Funding
California Air Resources Board and the California Environmental Protection Agency, agreement no. 19RD030 and the Health Effects Institute award no. #4979-RFA20-1B/21-2.

Declaration of Competing 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
The data that has been used is confidential.