We applied a Monte Carlo simulation approach to characterize the likely, daily exposure incidence to hazardous occupational heat from 2010 to 2019 for outdoor workers across all census tracts in the contiguous US. To achieve this, we integrated information on (i) the types of jobs people do and where, (ii) the typical work characteristics of these jobs that influence susceptibility to heat exposure, (iii) the typical heat exposures where people are working and (iv) whether these workers meet the American Conference of Governmental Industrial Hygienists (ACGIH) criteria for likely hazardous heat exposure.

We leveraged (i) employment counts by major occupational groups from the US American Community Survey, (ii) work characteristics from the Occupational Information Network (O*NET), and (iii) county-level data on daily-average wet bulb globe temperature (WBGT) for the contiguous US only, respectively. This approach estimated daily exposure incidence for every census tract in the contiguous US, which we aggregated both spatially (i.e., county, region, nationwide) and temporally (i.e., monthly, yearly). We characterized exposure estimates at the census tract level, rather than county, because our goal was to construct the most spatially granular estimates that are possible with publicly available data sources. A visual overview of this approach is displayed in Fig. 2.

Example of census tract-level estimates are displayed for three large US cities: (a) Los Angeles, California; (b) Houston, Texas; and (c) New York City, New York. Note: WBGT wet bulb globe temperature, BLS OEWS Bureau of Labor Statistics Occupational Employment and Wage Statistics, DoL O*NET Department of Labor Occupational Information Network. US ACS United States American Community Survey.

Occupations are commonly classified in the US following the US Standard Occupational Classification (SOC) structure. At the most granular level in 2010, non-military/civilian occupations are classified in 820 detailed SOC codes (e.g., security guards), which are grouped in 22 major SOC codes (e.g., protective services). Data in this study were integrated from multiple, disparate sources, each with varying levels of spatial, temporal, and occupational resolution. To our knowledge, the most spatially granular employment count data at the detailed SOC code level is the US Bureau of Labor Statistics (BLS) Occupational Wage and Employment Statistics (OEWS). This data source provided 2019 metropolitan/nonmetropolitan area employment counts of detailed Standard Occupational Classification (SOC) codes [14]. To link this data source with the county-level WBGT data, we assumed that the employment counts by detailed SOC codes were the same for the counties within a metro/nonmetro area. Overall employment counts were not altered at the county level because our main interest in this data source is the relative occupational composition of detailed SOC groups in a county. Assuming detailed SOC employment counts did not vary day-to-day and year-to-year from 2010 to 2019, we assigned each county-level, detailed SOC employment count a county-level, daily, mean WBGT estimate using previously published estimates (Supplementary Fig. B1) [15], which were available for the contiguous US only. We made this assumption because OEWS estimates are aggregated over time and are not intended to be assessed temporally.

We then assigned heat-related work characteristics to each detailed SOC. These were outdoor work, long/irregular work schedules, metabolic rate, and work/rest allocation, all derived from the 2019 O*NET (version 24.1) [16]. While indoor workers can also be at increased risk of heat stress, we only considered environmental heat exposure among outdoor workers in this analysis. Identifying those who could work long or irregular work schedules helped us identify potential occupations that might work on weekends, which is important to not bias our exposure estimates low by only considering weekdays. Metabolic rate and work/rest allocation were included because these two characteristics influence the WBGT threshold that places workers at risk of heat stress. Clothing is another important factor. For example, double-layered woven clothing may increase a worker’s internal temperature, while loose clothing likely will not. We did not consider clothing adjustment values because this information is not available in the O*NET. More details are provided in Appendix A of the Supplementary Materials regarding how each characteristic was derived. Long/irregular work schedules, metabolic rate, and work/rest allocation were each independent characteristics applied to a detailed SOC group; outdoor work was a probabilistic factor (i.e., the proportion of days in a year that a detailed SOC group works outdoors, \(p\)). Thus, for each day, we conducted Monte Carlo simulation with 200 iterations to assign a binary outcome of whether the outdoor workers may have been working outdoors on that day (i.e., \( \sim (200,p)\)). We chose 200 iterations because increasing iterations did not substantially improve estimation convergence (Appendix A, Figure A3).

If the simulation indicated that a detailed SOC group within a given county was working outdoors for a particular day, we then determined if the workers could be potentially exposed to hazardous heat using the ACGIH screening criteria for unacclimatized workers [17], which was determined as a function of a detailed SOC group’s metabolic rate, work/rest allocation, and the daily mean WBGT level (specific thresholds are displayed in Appendix A, Table A1). This assumes that all workers are unacclimatized to heat, which overestimates the number of individuals actually at risk of heat strain. However, the number of workers truly acclimatized to heat is unknown, so assuming all workers are unacclimatized provides a more realistic estimation of those potentially exposed to hazardous heat and is more inclusive of workers at increase susceptibility of heat strain. As a sensitivity analysis, we also re-ran our simulations assuming all workers are acclimatized (see Appendix B, Figs. B6, B8).

Thus, for each iteration-day, we estimated the number of workers exposed to hazardous levels of heat by detailed SOC group for each county in the contiguous US. The central estimate was the median across all 200 iterations. The unlikely, best-case scenario was the 2.5th percentile, while the unlikely, worst-case scenario was the 97.5th percentile. We aggregated these detailed SOC estimated counts of workers exposed to the major SOC level. For each major SOC group, we thus estimated the proportion exposed (i.e., number of workers exposed divided by total number of workers) every county-day. Since our primary goal was to estimate occupational heat exposure at the most spatially granular level possible (i.e., census tract), we applied these daily-county-major SOC proportions to the census tract level using annual employment counts by major SOC group for each census tract using poststratification. Thus, for each day from 2010–2019, we estimated the proportion of workers potentially exposed to hazardous levels of occupational heat for every census tract in the contiguous US.

For each year from 2010 to 2019, we extracted census tract-level estimates from 5-year, US American Community Survey data of the percentage of (i) racial and ethnic minorities (i.e., individuals who do not identify as non-Hispanic White); (ii) low-income individuals (i.e., <200% of the federal poverty limit); (iii) individuals with low educational attainment (i.e., without a high school diploma), and foreign-born populations. We further broke down the percentage of racial and ethnic minorities to four separate groups: non-Hispanic Black, non-Hispanic Asian, Hispanic, and all other minority race and ethnicity groups. To control for the likely non-working population of a census tract in our regression modeling, we extracted census tract-level estimates of the percentage of (i) individuals <16 years old, (ii) individuals >64 years old, and (iii) individuals who are unemployed, in the armed forces, or not in the labor force.

All analyses were conducted in R version 4.2.1. We summarized and visualized our estimations in a variety of ways. Treating each worker-day of exposure as an incident event, we aggregated the daily exposure proportions to monthly and yearly exposure rates per 100 workers.

$$=\frac__}_}}_\times }\times 100$$

for i days, j census tracts. An example interpretation of a rate of “10 exposures per 100 workers in 2010” could be “there were, on average, 10 exposures to occupational heat per day for every 100 workers in 2010”. We also assessed the number of worker-days of exposure across various temporal levels. These estimates could be aggregated to the county, state, region, and nationwide level (which could be further broken down by major SOC group) to show disproportionate burdens. Yearly estimates by census tract, county, state, and nationwide are publicly available at https://doi.org/10.5281/zenodo.15199026.

To test our hypothesis that marginalized populations of workers are disproportionately burdened by hazardous heat exposure at work, we ran four separate Poisson generalized additive models (GAM) using the yearly census tract-level percentage of racial and ethnic minority groups, low-income individuals, individuals with low educational attainment, and foreign-born individuals as the main, independent covariate of each model (using a smoothing spline with 4 degrees of freedom) and the estimated number of yearly worker-days exposed to occupational heat from 2010 to 2019 as the main outcome for each model. We included an offset term for the total number of residents in a census tract. We adjusted our models for the yearly, census tract-level percentage of individuals <16 years old, individuals >64 years old, and individuals who are unemployed, in the armed forces, or not in the labor force to control for individuals unlikely to be working. Models included a categorical variable for year to account for temporal autocorrelation and the census tract centroid’s latitude/longitude (tensor product with 40 degrees of freedom) to account for spatial autocorrelation.

We ran an additional regression model using the aforementioned four racial and ethnic subgroups (i.e., Hispanic, non-Hispanic Black, non-Hispanic Asian, and all other racial and ethnic minority subgroups) as the main independent covariates to assess whether specific minority groups may be highly exposed to occupational heat.