Potential regional air quality impacts of cannabis cultivation facilities in Denver, Colorado

The legal commercialization of cannabis for recreational and medical use has effectively created a new and almost unregulated cultivation industry. In 2018, within the Denver County limits, there were more than 600 registered cannabis cultivation facilities (CCFs) for recreational and medical use, mostly housed in commercial warehouses. Measurements have found concentrations of highly reactive terpenes from the headspace above cannabis plants that, when released in the atmosphere, could impact air quality. Here we developed the first emission inventory for cannabis emissions of terpenes. The range of possible emissions from these facilities was 66–657 t yr⁻¹ of terpenes across the state of Colorado; half of the emissions are from Denver County. Our estimates are based on the best available information and highlight the critical data gaps needed to reduce uncertainties. These realizations of inventories were then used with a regulatory air quality model, developed by the state of Colorado to predict regional ozone impacts. It was found that most of the predicted changes occur in the vicinity of CCFs concentrated in Denver. An increase of 362 t yr⁻¹ in terpene emissions in Denver County resulted in increases of up to 0.34 ppb in hourly ozone concentrations during the morning and 0.67 ppb at night. Model predictions indicate that in Denver County every 1000 t yr⁻¹ increase in terpenes results in 1 ppb increase in daytime hourly ozone concentrations and a maximum daily 8 h average (MDA8) increase of 0.3 ppb. The emission inventories developed here are highly uncertain, but highlight the need for more detailed cannabis and CCF data to fully understand the possible impacts of this new industry on regional air quality.

1 Introduction

The rapid expansion of one of the United States' newest industries, the commercial production and sale of recreational cannabis, was recently likened to the millennial “dot-com” boom (Borchardt, 2017). With an increasing number of states passing bills to legalize recreational cannabis, the enterprise is set to rival all but the largest of current businesses. The cultivation, sale, and consumption of recreational cannabis annual sales revenues had reached USD 1.5 billion in the US state of Colorado by 2017 (CDOR, 2018b), exceeding revenues generated by grain farming in the state. The commercial cultivation and sale of cannabis is not subject to the same strict environmental monitoring and reporting procedures as other industries of similar size. While the relaxation of laws has provided certain medicinal and economic opportunities for the states involved, the potentially significant environmental impact on air quality due to the production of cannabis has largely been ignored.

Previous research on the wider impacts of cannabis production has been limited due to its federal status as an illegal or controlled substance (Crick et al., 2013; Eisenstein, 2015; Andreae et al., 2016; Stith and Vigil, 2016). As a result of this status, most studies have focused on the pharmacological and health effects of the psychoactive constituents of Cannabis spp. (Ashton, 2001; Borgelt et al., 2013; WHO, 2016), or the societal impacts associated with the illicit nature of the industry (IDCP, 1995; Sznitman and Zolotov, 2015; WHO, 2016). The few assessments to date on the environmental impacts of the production of Cannabis spp. have centered on the detrimental effects of outdoor cultivation on ecosystems and watersheds due to land clearance and high water demand (Bauer et al., 2015; Carah et al., 2015; Butsic and Brenner, 2016). Studies have also quantified the energy consumption of the industry and the resulting greenhouse gas emissions associated with indoor cultivation (Mills, 2012). Little attention has been paid to the possible biogenic volatile organic compounds (BVOCs) emitted from the growing of cannabis and its impact on indoor and outdoor air quality.

The only studies that have measured the composition of gaseous emissions from cannabis have been limited to headspace samples above the plants (Hood et al., 1973; Turner et al., 1980; Martyny et al., 2013). These studies have shown high concentrations of VOCs such as monoterpenes (C₁₀H₁₆), sesquiterpenes (C₁₅H₂₄), and cannabinoids. These studies also measured thiols, a sulfur-containing compound responsible for the characteristic odor of Cannabis spp. (Rice and Koziel, 2015a, b). The principle (trace) components are reported to be α- and β-pinene, β-myrcene, d-limonene, cis-ocimene, β-caryophyllene, β-farnesene, and α-humulene (Hood et al., 1973; Turner et al., 1980; Hillig, 2004; Fischedick et al., 2010; Martyny et al., 2013; Marchini et al., 2014; Rice and Koziel, 2015a). The precise mix of chemical species, however, was strongly dependent on strain and the growing conditions (Fischedick et al., 2010). It should be noted that the pharmacologically active ingredients, e.g., tetrahydrocannabinol (Δ⁹-THC), generally have low-volatility and therefore are rarely detected in the gas phase (Martyny et al., 2013). Measurements in (illicit) cannabis cultivation facilities (CCFs) in conjunction with law enforcement raids in Colorado in 2012 found VOC concentrations of terpenes to be 50–100 ppb within growing rooms (Martyny et al., 2013). In these cases, the CCF operation contained fewer than 100 plants, compared with the thousands of plants found in currently licensed premises (CDOR, 2018a). Further, the Spokane Regional Clean Air Agency (SRCAA) study in Washington state measured indoor VOCs in seven flowering rooms and two dry bud rooms across four different CCFs. The average terpene concentration was 361 ppb (27–1676 ppb) in those facilities (Southwellb et al., 2017). These indoor measurements indicate the presence of BVOCs, but only limited studies have actually determined the chemical profile of gases actually emitted by the growing plants. For comparison, summertime outdoor monoterpene concentrations in forested regions of Colorado are typically less than 4 ppb (Ortega et al., 2014).

Terpenoids, such as monoterpenes (C₁₀H₁₆) and sesquiterpenes (C₁₅H₂₄), are highly reactive compounds with atmospheric lifetimes ranging from seconds to hours (Fuentes et al., 2000; Seinfeld and Pandis, 2006). They are primarily biogenic in origin (Fuentes et al., 2000; Guenther et al., 2012) and their reactions alter the atmospheric oxidizing capacity, resulting in a range of low volatility products that can partition into the aerosol phase and, depending on the concentration of nitrogen oxides (NO_x), lead to the formation of ozone (Laothawornkitkul et al., 2009; Guenther et al., 2012). Both ozone and aerosols are climate-relevant components of the atmosphere as well as critical air pollutants (USEPA, 2016).

In Colorado, the commercial growing of Cannabis spp. is restricted to secure and locked premises, resulting in indoor operations in most counties (CDOR, 2018a). Since legalization, the number of cannabis cultivation facilities (CCFs) has risen to 1400 across the state of Colorado in 2018, including more than 233 registered recreational and 375 medical CCFs within the Denver city limits alone. In Denver, the CCFs are commonly housed in commercial warehouses and the majority of these are located near transport links such as train hubs and major interstate highways (CDOR, 2019; Mills, 2012). Denver and the Front Range area are currently classified as “moderate” non-attainment of the ozone standard (USEPA, 2017). Due to that status, a federally mandated State Implementation Plan (SIP) was developed and mutually agreed upon between the state of Colorado and the United States Environmental Protection Agency (EPA) (CDPHE, 2009). Under the terms of the SIP, Colorado Air Quality Control Commission (AQCC) developed regulatory models to predict reductions in ozone precursors (CDPHE, 2009). These studies have found that ozone concentrations in Denver are VOC-sensitive, meaning that an increase in VOC concentrations will increase ozone production (UNC-IE and ENVIRON, 2013). The location of CCFs in a VOC sensitive region in Denver suggests a potential emission source that may impact regional air quality (UNC-IE and ENVIRON, 2014). This work used the best available information to produce the first emission inventory of VOCs from CCFs in Colorado. Colorado's regulatory model was then used to determine the extent that these emissions could impact regional air quality.

2 Materials and methods

2.1 Emission rate calculation

Figure 1a shows the locations of the licensed 739 recreational and 733 medical CCFs in Colorado as of March 2018 (CDOR, 2018a). Equation (1) was first used to estimate an emission rate for each CCF, and then all CCFs were used to build a bottom-up BVOC emission inventory.

$$\begin{array}{}\text{(1)}& {\mathrm{ER}}_i=\sum _j{\mathrm{EC}}_{ij}\times {\mathrm{DPW}}_{ij}\times {\mathrm{PC}}_{ij},\end{array}$$

where ER_i (µg h⁻¹) is the total emission rate for CCF i based on the sum of emission rates for all j cannabis strains; EC_ij is the emission capacity ($\mathrm{\mu }\mathrm{g}{\mathrm{dwg}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{1}}$; dwg is dry weight in grams) for cannabis strain j in facility i, DPW_ij is the dry plant weight in grams per plant for cannabis strain j, and PC is the plant count number for strain j in facility i.

Figure 1 (a) The locations of medical (red) and retail (green) Cannabis cultivation facilities (CCFs) in Colorado as of 1 March 2018. The corresponding values are the number of CCFs found within each city. The base map layer of this figure was made by Esri (Esri et al., 2013). (b) The 36 km ×36 km resolution of Western Air Quality Study (WAQS) and nested inner 12 km ×12 km resolution domains and 4 km ×4 km resolution domain used by the Comprehensive Air Quality Model with extensions (CAMx). This map was made by ENVIRON and Alpine Geophysics (ENVIRON and Alpine Geophysics, LLC, 2017b).

Since state legalization only occurred in 2014, and given the current federal illicit status of Cannabis spp., there is a lack of available data for the three parameters used in Eq. (1). The following describes the assumptions made for a range of potential values of EC, DPW, and PC given the best information available.

2.1.1 Emission capacity (EC)

The only data of EC from a leaf enclosure measurement are of three strains, namely Critical Mass, Lemon Wheel, and Rockstar Kush, that were 45 d old (Wang et al., 2019). This study found that at this growth stage the EC for total monoterpenes varied among strains: 10 $\mathrm{\mu }\mathrm{g}{\mathrm{dwg}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{1}}$ for Critical Mass, 7 $\mathrm{\mu }\mathrm{g}{\mathrm{dwg}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{1}}$ for Lemon Wheel, and 6 $\mathrm{\mu }\mathrm{g}{\mathrm{dwg}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{1}}$ for Rockstar Kush. The Department of Revenue (DOR) in Colorado has classified Cannabis spp. in a CCF into four different growth stages: immature (0–24 d old), vegetative (25–79 d old), flowering (80–132 d old), and at harvest (132–140 d old) (Hartman et al., 2018a). Wang et al. (2019) only sampled during the vegetative stage, and to our knowledge emission rates of monoterpenes from buds or flowers do not exist. It is not known how much EC will change during these different growth stages, but the grey literature does report that CCFs actively select cultivars to maximize the amount of monoterpenes found in the bud tissues.

The Spokane Regional Clean Air Agency (SRCAA), in collaboration with Washington State University (Southwellb et al., 2017; Wen et al., 2017), measured monoterpenes in flowering rooms of CCFs in Washington State. They found concentrations of monoterpenes in the growing room with 80 d old plants (1660 ppb) to be >10 times higher than the 48 d old plants (150 ppb). CCFs in Colorado house a wide variety of strains at both vegetative and flowering stages of growth, suggesting that the emission rate of monoterpenes from CCFs is higher than that measured from foliage by Wang et al. (2019). Currently, no database exists that can provide the number of plants by strain and growth stage. For the base case, it was assumed that each CCF grew only one strain and that all plants were at the vegetative growth stage, resulting in a single and constant EC for each CCF, taken to be 10 $\mathrm{\mu }\mathrm{g}{\mathrm{dwg}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{1}}$ of total monoterpenes based on the reported EC from the Critical Mass cultivar (Wang et al., 2019). Given the uncertainty in EC and the variety of possible plant stages and cultivars, the EC used in simulation 1_EC was multiplied by a factor of 5 and 10 in simulations 2_EC and 3_EC as a sensitivity analysis.

2.1.2 Dry plant weight (DPW)

No published studies report the DPW of a Cannabis spp. plant. Both the states of Colorado (METRC, 2018) and Washington (LCB, 2017; Topshelfdata, 2017) track the mass of the commercially sold portion of the plant, the “dry bud”. The Colorado database, however, is not publicly accessible and was not available for this study. In Washington, using data from all types of facilities (outdoor and indoor) from August to October 2017, it was found that the average dry bud mass per plant was 210 g (0–586 g) shown in Fig. S1a in the Supplement. The Washington database also includes the “wet bud” weight defined as the mass of the bud after it was just harvested (Fig. S1b in the Supplement), but prior to the 7–10 d drying process. The total waste weight, or the remaining mass of the plant after the buds have been harvested, is also recorded. As shown in Eq. (2), the sum of these two masses should equal the total mass of the wet plant.

$$\begin{array}{}\text{(2)}& {\displaystyle }{M}_{\text{wet plant}}=M_{\text{wet buds}}+M_{\text{waste}},\end{array}$$

where M_{wet plant} is the mass of the entire wet plant (g), M_{wet bud} is the mass of the wet bud (g), and M_{wet waste} is the mass of the waste (g).

Data from August to October 2017 were used with Eq. (2) to estimate the wet plant weight resulting in an average of 3770 g (6–13 405 g) shown in Fig. S1c. The large range in mass is due to the different growing conditions found in CCFs, and the type of strain being grown. The ratio of the wet and dry bud mass data from Washington was used as a surrogate to determine the percentage of water found in the total plant material as shown in Eq. (3).

$$\begin{array}{}\text{(3)}& {\displaystyle }{R}_{\mathrm{D}/\mathrm{W}}=M_{\text{dry bud}}/M_{\text{wet bud}},\end{array}$$

where R_D∕W is the ratio of the masses of the dry to wet bud, and M_{dry bud} (g) is the mass of the harvested buds after 7–10 d of drying (Fig. S1d).

It was assumed that the same factor could be applied to the total wet plant weight to estimate the DPW as shown in Eq. (4).

$$\begin{array}{}\text{(4)}& {\displaystyle }\mathrm{DPW}=M_{\text{wet plant}}\times R_{\mathrm{D}/\mathrm{W}}\end{array}$$

The average of DPW was 754 g (1–2260 g). For the development of these emission inventories, a base value of 750 g was assumed for DPW based on the average calculated from the Washington database. As a sensitivity test, a DPW of 1500 g representing the mean plus 1 standard deviation range was chosen. Finally, a DPW of 2500 g, the maximum yield recorded by Washington State Liquor and Cannabis Board, was taken as the upper statistical boundary as shown in Fig. S1e. As the total plant count and reported yields are a factor of 3 and 4 higher, respectively, in Colorado than Washington State (LCB, 2017; Topshelfdata, 2017; Hartman et al., 2018a), we took this maximum on the assumption that Cannabis spp. cultivated in CCFs in Colorado in the summer season are grown under more optimal conditions than those grown in Washington State, resulting in considerably higher yields.

2.1.3 Plant count (PC)

Counts of all plants larger than 20.3 cm have been recorded by the Colorado DOR on a monthly basis since 2014. As of June 2018, there are a total of 1.06 million plants (Hartman et al., 2018a, b). We therefore used 1 million as the base number for the emission inventory. The DOR data only provide county-level information rather than actual number of plants per CCF. The plants were then distributed equally among the CCFs to calculate an average of 905 plants per facility in Denver County and 521 outside of the county.

Two sensitivity simulations were conducted based on the assumption that the cannabis industry in Colorado will continue to expand at similar rates in the future. From June 2016 to June 2018 the total number of plants recorded by DOR grew from 826 963 to 1 062 765, an annual average increase of 118 000. Assuming this rate of expansion remains constant, there would be 2 million plants in the state of Colorado by 2025 and this value was used in simulation 6_PC. It was assumed in simulation 7_PC that growth would accelerate in the future to the point at which each recreational and medical CCF would contain the maximum number of plants permitted under a Tier 1 license leading to a statewide total of nearly 4 million plants. The maximum number of plants that can be grown under each licensing tier is shown in Table S2 in the Supplement (CDOR, 2019). The average plant count per CCF for each PC sensitivity simulation is shown in Table S1.

2.2 Emission inventories for cannabis cultivation facilities (CCFs)

Given the large gaps in knowledge, this study will focus only on variabilities in EC, DPW, and PC and will hold other parameters constant. For example, to maximize growing conditions relative humidity, temperatures, CO₂ concentrations, and fertilizer usage are all optimized and vary widely by CCF. Further, this study did not consider other processes such as trimming, harvesting, and drying buds, which may also release BVOCs.

For this study, it was assumed that all CCFs operated in the same way at a temperature of 30 ^∘C and 1000 $\mathrm{\mu }\mathrm{mol}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$ of photosynthetically active radiation (PAR). In addition, it was assumed that all emissions from the plants inside a CCF enter the atmosphere. Ventilation to the atmosphere varies widely by the operation, and there are no current regulations or industry-wide practices that are being used to mitigate emissions.

In total, seven scenarios of emission inventories were created to explore sensitivities in EC, DPW, and PC as shown in Table 1. In scenarios 1–3, the PC was held to a total of 1 million and a 750 g DPW was assumed. The EC of 10 $\mathrm{\mu }\mathrm{g}{\mathrm{dwg}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{1}}$ as reported by Wang et al. (2019) was used in 1_EC, with a sensitivity that multiplied that rate by a factor of 5 (scenario 2_EC) and 10 (scenario 3_EC). The remaining scenarios in Table 1 kept the EC constant at 10 $\mathrm{\mu }\mathrm{g}{\mathrm{dwg}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{1}}$. Scenarios 4_DPW and 5_DPW explored the sensitivity of increasing DPW, and scenarios 6_PC and 7_PC increased the total plant count.

Table 1 Simulation scenarios and assumed values for emission capacity (EC) rate, dry plant weight (DPW), and the plant count (PC) for Colorado and Denver County. The base case (BC) scenario has no cannabis emissions.

2.3 Model description and analysis tools

2.3.1 Model protocols and evaluation

The Comprehensive Air Quality Model with Extensions, CAMx6.10 (ENVIRON, 2013; ENVIRON and Alpine Geophysics, LLC, 2017b), was used to predict ground-level ozone concentrations. The model and protocols used in this study are based on the Western Air Quality Study (WAQS) for 2011 (ENVIRON and Alpine Geophysics, LLC, 2017b; Adelman et al., 2016). The WAQS 2011b baseline model simulation period runs from 15 June to 15 September 2011, and is driven with meteorological data from WRF version 3.3 for the same time period and domain. The model was initialized using Three-State Air Quality Modeling Study standard boundary and initial conditions (ENVIRON and Alpine Geophysics, LLC, 2017b). The model domain is a two-way nested grid at 12 and 4 km grid cell resolutions (Fig. 1b). Anthropogenic emissions were derived from EPA National Emission Inventory (NEI) version 2011 NEIv2 with updates for point and area sources of oil and gas emissions in the western US. The biogenic emission inventory was based on the Model of Emissions of Gases and Aerosols from Nature version 2.1 (MEGANv2.1) (Guenther et al., 2012). All data and supporting documentation are publicly available via the Intermountain West Data Warehouse (IWDW) website (WAQS, 2017).

The revision 2 of the Carbon Bond 6 (CB6r2) (Ruiz and Yarwood, 2013) chemistry mechanism was used in all model runs. This groups all monoterpenes as a single compound species, TERP. Thus, the total monoterpene EC reported in Wang et al. (2019) was converted into the TERP species. TERP undergoes oxidation reactions with the nitrate radical (NO₃), the hydroxyl radical (OH), ozone (O₃), and singlet oxygen. It should be noted that the TERP category includes a wide variety of monoterpenes whose reaction rate constants may differ from TERP ($k_{\mathrm{298}}=\mathrm{6.77}\times {\mathrm{10}}^{-\mathrm{11}}$ molecules cm⁻³ s⁻¹). For example, the rate constant of β-myrcene with OH radical (Hites and Turner, 2009) is $\mathrm{3.35}\times {\mathrm{10}}^{-\mathrm{10}}$ molecules cm⁻³ s⁻¹ (k₂₉₈), which is 4 times higher than TERP and 5.6 times faster than α-pinene (Carter, 2010).

The details of the WAQS model setup protocol (ENVIRON and Alpine Geophysics, LLC, 2017b) and model performance (Adelman et al., 2016) can be found on the IWDW website. In summary, the model performance evaluation concluded that this simulation had met all performance goals for both maximum daily 1 h (MDA1) and maximum daily 8 h average (MDA8) ozone. In the performance review report, it was found that the WAQS model had a positive bias for ozone simulated in a 4 km ×4 km resolution domain, when compared with EPA Air Quality System (AQS) surface monitors (MDA1: 0.8 %, MDA8: 0.9 %). On days when ozone concentrations higher than 60 ppb were measured, the model had a negative bias of −6.2 % for MDA1 and −6.3 % for MDA8. The model evaluation result also noted that the model performance was best during the spring and summer months.

Table 2 The estimated BVOC and total VOC emission rates (t yr⁻¹) for the base case (BC) scenario. Also shown are the increases in VOC emissions for all scenarios shown in Table 1 for Colorado, Denver County, Colorado Springs, Pueblo, and Boulder. The numbers in parentheses are the percentage increases compared with the BC scenario.

2.3.2 Process analysis

CAMx runs used in this analysis had the process analysis (PA) option enabled (ENVIRON, 2013). The CAMx configuration used here produces two additional files needed for PA: the integrated reaction rate (IRR) and integrated process rate (IPR). These files include the rates of change in concentration of every species due to chemistry and transport for every grid cell and time step. Python-based Process Analysis (pyPA) and the Python Environment for Reaction Mechanisms/Mathematics (PERMM) (Henderson et al., 2010, 2011) were then applied to post-processing the CAMx PA output. PERMM was used to aggregate the chemical and physical process rates for selected model grid cells and layers, allowing for tracking of plumes within the planetary boundary layer (PBL).

3 Results

3.1 Emission inventory

The seven scenarios were used to estimate a range of emissions of BVOCs from CCFs for the entire state of Colorado. As shown in Table 2, the base case (BC) scenario estimates 731 442 t yr⁻¹ of all VOCs being emitted in Colorado, of which 47 % are BVOCs. The BC scenario does not include any emissions from the cannabis industry. Table 2 also shows the seven scenarios that did include CCF emissions ranked in order of their increases in statewide BVOC emissions. As expected the CCF BVOC emissions scaled linearly with each factor that was changed in Eq. (1). In scenario 3_EC, a 10-fold increase in the emission rate (100 $\mathrm{\mu }\mathrm{g}{\mathrm{dwg}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{1}}$) resulted in a 657 t yr⁻¹ increase. Similarly, scenario 2_EC assumes 50 $\mathrm{\mu }\mathrm{g}{\mathrm{dwg}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{1}}$ and produces 329 t yr⁻¹. Scenarios 4 and 5 showed the sensitivity of terpene emissions from CCFs to variation in DPW while holding PC constant and an EC of 10 $\mathrm{\mu }\mathrm{g}{\mathrm{dwg}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{1}}$. It was estimated that an additional 66 t yr⁻¹ of emissions is produced when a 750 g DPW is assumed. This doubles to 131 t yr⁻¹ with a DPW of 1500 g and reaches 219 t yr⁻¹ with a DPW of 2500 g. Comparing scenario 1_EC with scenarios 6 and 7 shows how the growth in PC will impact emissions of BVOCs. In Colorado, a doubling of the PC increases BVOC emissions by 131 t yr⁻¹ in scenario 6_PC and 261 t yr⁻¹ for the 4 million plants in scenario 7_PC. The largest increases in BVOC emissions were predicted in scenarios 3_EC and 2_EC showing that the total emission rate of BVOCs from CCFs was most sensitive to EC.

In March 2018, Denver County housed 41 % of CCFs and 55 % of all cannabis plants in Colorado (Hartman et al., 2018b). As a result, about 43 % of statewide CCF BVOC emissions occur there (Table 2). Current emission inventories of Denver County show negligible amounts of biogenic emissions accounting for only 0.1 % of the total statewide BVOC emissions. CCF emissions increased BVOC emission rates in Denver Country up to 136 % in scenario 3_EC. This changes the total VOC emission rate in Denver County by up to 3.5 %. Other cities in Colorado do not have as high of a concentration of CCFs, and thus the relative increases were smaller as shown in Table 2.

The introduction of additional cannabis BVOC emissions into model simulations increased the predicted TERP concentrations. Figure 2 shows the maximum increase in TERP concentrations for three scenarios for Denver County over the entire 90 d simulation period. Regardless of the scenario, the largest increases in TERP occurred near the largest concentrations of CCFs. The absolute maximum changes ranged from 0.5 to 5.0 ppb located at the Elyria–Swansea and Globeville neighborhoods in north-central Denver. Increases in TERP were also predicted to the north due to the dominant wind flows in that direction throughout the simulation period. Figure S2 shows the maximum increase in TERP concentrations for the 1_EC, 5_DPW, and 3_EC scenarios in the 4 km ×4 km domain for the entire 90 d simulation period. As expected substantially lower increases in TERP concentrations were predicted for other cities in Colorado: 0.26 ppb in Colorado Springs and 0.24 ppb in Pueblo. Figure 3 shows the hourly changes in TERP concentrations across the entire 4 km ×4 km domain. The largest increases for all scenarios occurred at night with a peak of 5 ppb at 04:00 local standard time (LST). Given that the hourly emissions of terpenes from CCFs were assumed constant for 24 h, these larger nighttime changes can be primarily ascribed to the lack of photochemistry and a shallow nocturnal PBL. These results suggest that the increases in TERP are highly correlated with locations of CCFs, accumulate at night, and have significant losses during the day.

Figure 2 The maximum increase in TERP concentrations (ppb) for Denver County and Front Range over the entire 90 d simulation for the (a) 1_EC, (b) 5_DPW, and (c) 3_EC scenarios. The black outlines Denver County and the grey lines are state and interstate highways.

Figure 3 The hourly changes in TERP concentrations across the entire 4 km ×4 km domain, over the 90 d simulation for the (a) 1_EC, (b) 5_DPW, and (c) 3_EC scenarios.

3.2 Regional ozone impacts

Predicted increases in hourly ozone concentrations in excess of 0.1 ppb only occurred when terpene emissions were in excess of 219 t yr⁻¹, with scenarios 4_DPW, 6_PC, and 1_EC having little impact on predicted ozone. Thus, this analysis will focus on two scenarios, 5_DPW and 3_EC, to explore potential regional ozone impacts in the present and future. Figure 4 shows the hourly changes in ozone concentrations across the entire 4 km ×4 km domain for these two scenarios. During the daytime, the increase in TERP emissions results in a peak ozone increase of 0.34 ppb at 09:00 LST for 3_EC with only minimal changes in 5_DPW. Figure 5 shows, for Denver County and the Front Range metropolitan area, the locations of the daytime (06:00–18:00 LST) maximum increases in hourly ozone concentrations for all 90 d when emissions were added for scenarios 5_DPW and 3_EC. Ozone increases for the entire 4 km ×4 km domain can be found in Fig. S3. The largest predicted ozone concentrations occurred in Denver County with impacts of 0.11 ppb in 5_DPW and 0.34 ppb in 3_EC as shown in Fig. 5. Both scenarios show that daytime increases in ozone were limited to Denver County and just to the northwest, west, and southwest of Denver County.

Figure 4 The predicted differences in hourly ozone concentrations (ppb) across the entire Colorado domain, over the 90 d simulation for the (a) 5_DPW and (b) 3_EC scenarios.

Figure 5 The predicted changes in hourly ozone concentrations for the Denver region from 06:00 to 18:000 LST for all 90 d of the simulation for the (a) 5_DPW and (b) 3_EC scenarios. The grey lines indicate major highways and the black line outlines Denver County.

There were also nighttime variations in ozone observed for the modeling domain. In scenarios 5_DPW and 3_EC, nighttime increases were more than double the increases predicted during the day. The largest changes in hourly ozone concentrations of 0.67 ppb occurred at 00:00 LST (i.e., midnight) for 3_EC. Figure 6 shows the location and magnitude of the maximum changes in hourly ozone concentrations during the night (18:00–06:00 LST) in 5_DPW and 3_EC. The extent of ozone increases at night is primarily to the north of Denver indicating a northern outflow. The maximum increase in hourly ozone for the whole of Colorado is shown in Fig. S3, with visibly little changes at night in other cities. These model results suggest that the additional emissions of TERP have immediate impacts on local ozone production chemistry during both the day and night, but little wider impact.

Figure 6 The predicted changes in hourly ozone concentrations for the Denver region from 18:00 to 06:00 LST for all 90 d of the simulation for the (a) 5_DPW and (b) 3_EC scenarios. Black regions within the map indicate ozone increase values greater than 0.5 ppb. The grey lines indicate major highways and the black line outlines Denver County.

A critical metric for the attainment of the National Ambient Air Quality Standards (NAAQS) ozone standard in Denver County is the maximum daily average 8 h ozone concentration (MDA8). Figure 7 shows the maximum difference in MDA8 for each grid cell centered on Denver County, across the entire 90 d simulation period for the 5_DPW and 3_EC scenarios. Maximum increases in MDA8 are 0.14 ppb for 3_EC (Fig. 7b) co-located with the maximum increases in TERP concentrations.

Figure 7 The predicted maximum increases in the maximum daily average 8 h (MDA8) ozone concentration (ppb) for the (a) 5_DPW and (b) 3_EC scenarios for the Denver region over the 90 d simulation period. The black indicates ozone increase values greater than 0.12 ppb.

3.2.1 Ozone impact at night

The maximum hourly ozone increase of 0.67 ppb for the 3_EC scenario occurred on Thursday, 28 July 2011, at 00:00 LST (i.e., midnight) near the largest concentration of CCFs (see Fig. 8). In subsequent hours the plume of ozone moved slowly to the east before being dispersed by the rise of the morning PBL at 06:00 LST.

Figure 8 For the 3_EC scenario on 28 July 2011, the largest hourly predicted ground-level ozone increases at (a) 27 July, 21:00 LST and for 28 July, at (b) 00:00 LST (i.e., midnight), (c) 03:00 LST, and (d) 06:00 LST.

To better understand why ozone increased at night, the PA model output was analyzed to quantify the chemical and physical processes producing ozone. Plume tracking was used so that only grid cells where the increase in ozone (i.e., the plume) occurred were included in our analysis, which ran from 27 July, 21:00 to 28 July, 06:00 LST. The number of vertical model layers included in the analysis also varied to incorporate the hourly evolution of the PBL. Figure S4 provides snapshots of the horizontal grid cells used and the vertical layers that were aggregated throughout the simulation time period. Figure S5 shows the changes in final ozone concentrations (compared to the base case) for the grid cells and vertical layers included in the analysis, as well as the physical and chemical process rates that account for these changes. Figure S5 shows that the process most responsible for increases in ozone concentrations was chemical production.

For the chosen vertical layers and grid cells Table 3a shows the total rate of the oxidation reactions with TERP across the entire period. Throughout this time, the additional TERP emissions lead to an increase in the number of oxidation reactions thereby generating more secondary VOC products and radical species. The chemical losses of TERP increased due to reactions with OH (from 0.01 to 0.1 ppb; +900 %), nitrate radical (NO₃) (from 0.39 to 1.58 ppb; +305 %), and O₃ (from 0.04 to 0.2 ppb; +400 %). Further analysis confirms that nighttime oxidation chemistry leading to changes in ozone concentration are driven by NO₃. In the 3_EC scenario, TERP emissions only increased the annual VOC emission in Denver County by 3.5 %, but this is sufficient to increase the VOC + NO₃ reaction rates by 125 %. These increases produce more peroxyl radicals (TRO₂=HO₂ + RO₂) driving further oxidation and further radical production. Table 3b also shows that the generation of OH radicals from reactions of TERP with O₃ increased by 267 %. Ultimately, these increases in initial TERP reactions with NO₃ and O₃ increase the NO-to-NO₂ conversions via the TRO₂ pathway by 44 %, reducing the availability of NO to react with O₃. Thus, the increased ozone concentration predicted at night is actually due to the 1 ppb (0.8 %) reduction in the loss of ozone to reactions with NO rather than an increase in actual production of ozone (Table 3c). The increased TERP emissions also increase production of NO_x termination products (NO_z) by 27 % with organic nitrate (NTR; representing ∼71 % of this NO_z product), increasing from 0.66 to 1.6 ppb (+142 %). This increase in NO_z production at night also results in lower NO concentrations and thus lower ozone titration.

Table 3 All data summed from 27 July, 21:00 LST to 28 July, 05:00 LST for grid cells and layers shown in Fig. S4. The base case (BC) scenario column shows the absolute predicted values and the subsequent columns show the predicted changes due to emissions from the 3_EC scenario. Percentages in parentheses are the changes in 3_EC relative to BC. Shown is the (a) total amount of VOC and TERP consumed due to oxidation (ppb), the (b) total amount of hydroxyl radical (OH) and total peroxyl radicals (TRO₂) that were generated and their sources (ppb), and the (c) total amount of nitrogen dioxide (NO₂) and NO_x termination products (NO_z) produced and their sources (ppb).

3.3 Ozone impact during the day

The maximum daytime hourly ozone increase of 0.34 ppb occurred at 09:00 on Monday, 18 July 2011, as shown in Fig. 9. On this day, the meteorological conditions favored the maximum possible production of ozone. This day featured “upslope flows” that are a common meteorological condition linked to ozone exceedance periods (Pfister et al., 2017). We thus chose to focus on 18 July to understand the daytime changes in chemistry that occur from increased BVOC emissions. As expected, the location of predicted ozone increases coincides with the location of the strongest terpene emissions in the domain as shown in Fig. 9a. For the daytime hours of 06:00–14:00 LST, the PA option was used to quantify changes in chemical processes for the grid cells and model layers shown in Fig. S6. For these grid cells and layers, Fig. S7 shows the changes in final ozone concentrations compared to the base case and the physical and chemical process rates that impact those concentrations. Table S3 sums the key chemical processes for these hours. The increases in CCF emissions resulted in a 100 % increase in OH reactions with TERP producing intermediate oxidation products and ultimately increasing OH production by 0.6 %. As a result of this oxidation chemistry, there was an increase of 0.9 % in NO-to-NO₂ conversion by the TRO₂ pathway, ultimately leading to a 0.1 % increase in ozone production.

Figure 9 For the 3_EC scenario on 18 July 2011 the largest hourly predicted ground-level ozone increases at (a) 09:00 LST, (b) 12:00 LST (i.e., noon), (c) 14:00 LST, and (d) 17:00 LST. The maximum of 0.34 ppb occurred at 09:00 LST.

3.3.1 Ozone impact sensitivity

The maximum modeled daytime hourly ozone increase due to additional CCF emissions occurred on 18 July. Using this day, multiple sensitivity simulations were performed, where CCF emissions from Denver County were incrementally increased up to 3800 t yr⁻¹. Figure 10 shows the increase in terpene emissions from Denver County versus the largest daily increase in hourly ozone concentrations. Figure 10a shows a linear relationship, indicative of a VOC-limited environment, where hourly ozone concentrations are predicted to increase by 1 ppb for every 1000 t yr⁻¹ increase in TERP emissions during the day and 0.85 ppb at night. Also shown is the sensitivity to the MDA8 ozone where there is a 0.30 ppb increase for every 1000 t yr⁻¹ of TERP emissions. According to projected emission inventories provided by the state of Colorado, the ozone non-attainment area was expected to see reductions of 26.4 % of NO_x and 24.6 % of VOC emissions by the year 2017 (ENVIRON and Alpine Geophysics, LLC, 2017a). Under these reduced anthropogenic emission scenarios, Fig. 10b shows how ozone would then respond to additional CCF TERP emissions. Figure 10b continues to show a linear relationship, where hourly ozone concentrations are predicted to increase by 1.5 ppb for every 1000 t yr⁻¹ increase in TERP emissions during the day and 1.8 ppb at night. In the future case, the MDA8 ozone increases by 0.38 ppb for every 1000 t yr⁻¹ of TERP emissions. Therefore, Denver will still be VOC-limited and ozone is predicted to be more sensitive to CCF emissions of terpenes.

Figure 10 For 18 July during (a) 2011 and (b) 2017 the predicted maximum increase in hourly ozone concentrations during daytime hours (06:00–18:00 LST) in blue and nighttime hours (18:00–06:00 LST) in black versus additional terpene emissions in Denver County. Also shown is the response in maximum daily average 8 h ozone concentration (MDA8) in red.

4 Conclusions

This study provides the first VOC emission inventory to be compiled for the cannabis industry in Colorado, the first time such an analysis has been conducted anywhere in the USA. Given the current state of knowledge of emission rates and growing practices, there are considerable uncertainties in the basic parameters required to build such an inventory. Using realistic bounds on each parameter, we developed seven scenarios, which resulted in estimated emission rates that ranged over an order of magnitude. The highest emissions occur in Denver County, with rates ranging between 36 and 362 t yr⁻¹ for the different scenarios, from a total of 66–652 t yr⁻¹ across Colorado as a whole.

We included these additional terpene emissions in the Comprehensive Air Quality Model with extensions (CAMx), the model used by the state of Colorado for regulatory monitoring and projections. Taking the worst case (3_EC) and median scenario (5_DPW) we consider representative of the current uncertainty upper boundary and future industry expansion, we find that these projected increases in emissions lead to maximum increases in terpene concentrations of up to 5.0 ppb. The largest impacts were seen in locations with the highest terpene emissions coming from CCFs, i.e., in Denver County. We further found that these increases in terpene concentrations affected the local atmospheric chemistry and air quality with ground-level ozone concentrations increasing by as much as 0.34 ppb during the day and 0.67 ppb at night. In general, simulated nighttime increases were higher than those during the daytime, and we take the nighttime of 27–28 July as a case study to further investigate. By applying process analysis (PA), following the evolving plume of VOCs and ozone, we find that the initial reactions of the additional terpenes with OH, NO₃, and ozone result in increased formation of peroxyl radicals, which increases the NO-to-NO₂ conversion rate and also removes the NO_x to generate more NO_z product. This effectively reduces the loss of ozone by reaction with NO, increasing the total ozone concentration.

We acknowledge, however, the considerable uncertainties that surround our projections and call for the need for continued efforts to reduce these such that a more accurate assessment of the regional air quality implications of this industry can be made. Future studies that include ambient BVOC measurements are critical for comparisons with model predictions. Additionally, in the model chemical mechanism more accurate and mechanistic representation of terpene species that can reflect the current cannabis emission composition is needed. Currently, the model surrogate TERP, which represents all monoterpene species in the mechanisms, may not represent the precise rate constant for BVOC emissions from cannabis. Further data are needed to reduce uncertainties in emission inventory estimates, specifically those regarding CCF-specific information on plant counts and weight by cultivar and growth stage, coupled with information about the agronomic practices of cannabis cultivation in CCFs. Additional measurements of emission capacities of different cannabis strains at different growth stages are also needed. Further, the emission inventory version is for the year 2011; it may not be suitable to estimate the ozone impacts by the CCF industry.

We chose to focus on ozone since Denver is a moderate non-attainment area with an ozone State Implementation Plan (SIP) (ENVIRON and Alpine Geophysics, LLC, 2017a, b; Colorado, 2018) in accordance with the EPA regulations. But assessments of the impact of these additional terpene emissions on particulate matter (PM_2.5) are warranted given the high secondary organic aerosol (SOA) yields of terpenes from 0.3 to 0.8 (Iinuma et al., 2009; Lee et al., 2006; Fry et al., 2014; Slade et al., 2017). It should also be borne in mind that investigations of indoor air quality are needed given the findings of Martyny et al. (2013) and Southwellb et al. (2017) that indoor terpene concentrations reached 50–100 ppb in growth rooms and 30–1600 ppb in flowering rooms, likely initiating intense photochemistry under the powerful growing lamps in use in CCFs.