Transport phenomena governing nicotine emissions from electronic cigarettes: Model formulation and experimental investigation

ABSTRACT

Electronic cigarettes (ECIGs) electrically heat and aerosolize a liquid-containing propylene glycol (PG), vegetable glycerin (VG), flavorants, water, and nicotine. ECIG effects and proposed methods to regulate them are controversial. One regulatory focal point involves nicotine emissions. We describe a mathematical model that predicts ECIG nicotine emissions. The model computes the vaporization rate of individual species by numerically solving the unsteady species and energy conservation equations. To validate model predictions, yields of nicotine, total particulate matter, PG, and VG were measured while manipulating puff topography, electrical power, and liquid composition across 100 conditions. Nicotine flux, the rate at which nicotine is emitted per unit time, was the primary outcome. Across conditions, the measured and computed nicotine flux were highly correlated (r = 0.85, p < .0001). As predicted, device power, nicotine concentration, PG/VG ratio, and puff duration influenced nicotine flux (p < .05), while water content and puff velocity did not. Additional empirical investigation revealed that PG/VG liquids act as ideal solutions, that liquid vaporization accounts for more than 95% of ECIG aerosol mass emissions, and that as device power increases the aerosol composition shifts towards the less volatile components of the parent liquid. To the extent that ECIG regulations focus on nicotine emissions, mathematical models like this one can be used to predict ECIG nicotine emissions and to test the effects of proposed regulation of factors that influence nicotine flux.

Copyright © 2017 American Association for Aerosol Research

EDITOR:

• Warren Finlay

1. Introduction

Electronic cigarettes (ECIGs) are a class of products that use an electrically powered heating element to aerosolize a liquid that usually contains propylene glycol (PG), vegetable glycerin (VG), flavorants, water, and the psychomotor stimulant drug nicotine (Hajek et al. 2014). ECIGs are controversial for a variety of reasons, including uncertainty regarding their ability to substitute completely for tobacco cigarettes, the health effects associated with long-term use, and their use by otherwise nicotine-naïve adolescents (Breland et al., 2016). Another controversial topic is the extent to which these products should be regulated to limit harm to individual and public health (Hajek et al. 2014; Pisinger 2014; McNeill et al. 2015; Tobacco Advisory Group 2016). Because nicotine is a drug that supports dependence (i.e., difficulty controlling drug use after repeated administration), understanding ECIG nicotine emissions becomes a key regulatory issue.

We have previously proposed “nicotine flux” (Shihadeh and Eissenberg 2015), the rate at which nicotine is emitted from an ECIG per unit time (e.g., mg/s), as a measure of potential regulatory value. Mapping the relationships between nicotine flux and ECIG design features and operating conditions can reveal various regulatory levers to control the potential ranges of nicotine intake by users (Shihadeh and Eissenberg 2015). Such a map could be constructed empirically by studying numerous product-use condition combinations. However, an engineering modeling approach can elucidate more efficiently key product parameters and provide an understanding of the underlying interdependent physical processes that govern ECIG nicotine flux. Models that allow prediction of nicotine flux may present an opportunity to constrain emissions by regulating the factors that control them, including device design features and liquid composition (Talih et al. 2015). Without such models, regulation may focus inappropriately on one factor, rather than the simultaneous operation of all factors (European Commission 2014). Here we develop a generic model that can be used to predict nicotine emissions based on a few ECIG design and operating conditions. Such a model could be used as a tool to predict emissions of existing and future products, and to analyze effects of regulations proposed for controlling ECIG nicotine emissions. Before describing the model, we review the basic phenomena involved in the ECIG aerosol production process.

1.1. Principles of ECIG operation

While ECIGs vary widely in design, the essential characteristic that defines this product category is the use of an electrical heating coil to vaporize a nicotine-containing liquid (Figure 1).

Figure 1. Schematic of an ECIG “tank” or “clearomizer” unit (battery not shown).

When the heater coil is activated during a puff, the liquid surrounding the coil heats up and vaporizes, and the vapors (and thermal energy) are convected away from the vicinity of the coil by the air drawn through the device. As the hot vapors come into contact with the cool air drawn into the device, vapors re-condense to form an aerosol mist that visually resembles tobacco smoke (Figure 2). Typically through the action of a wick, fresh liquid is transported to the coil to replace the liquid lost during vaporization.

Figure 2. Schematic of the ECIG aerosol generation process.

Key to the aerosol production process is the energy balance, or the disposition of the electrical energy input to the coil. While all the electrical energy initially is converted into thermal energy within the heater coil filament, this thermal energy in turn (1) heats liquid surrounding the coil, (2) converts liquid into vapor, (3) heats the air passing over the coil, and (4) flows by conduction along the wick and the electrical leads into the ECIG body. In general, the greater the energy dissipated from the coil to the ECIG body, the less will be available to produce vapor. Another important aspect of this process is its transient nature. When the heater coil is first energized, its temperature rises until the various thermal energy dissipation routes mentioned above collectively attain a rate that balances the rate at which electrical energy is input to the coil. The time taken to attain this steady state is a measure of the “thermal inertia” of the coil and its immediate surroundings; all else being the same, a more massive coil assembly will take longer to attain its steady state temperature.

ECIG liquid thermodynamics also play an important role in constraining the vaporization process. Typically, ECIG liquids are composed of a PG/VG mixture with varying amounts of nicotine, water, flavorants, and other constituents (Cheah et al. 2014; Pellegrino et al. 2012; Uryupin et al. 2013). If the temperature in the vicinity of the coil is below the bubble nucleation temperature (i.e., the “boiling point”) of the ECIG liquid, then the vaporization rate will be governed by evaporative convection. The evaporative convection rate is a function of the vapor pressures, mass diffusivities, and mole fractions of the evaporating species, the mass transfer coefficient (a function of the flow rate, geometry, and chemical species), and the area of the surface where evaporation occurs. If the coil reaches the boiling temperature of the solution, the vaporization rate will be governed strictly by the rate at which heat is delivered to the liquid, irrespective of flow conditions. Whether the system attains boiling or not depends on the electrical power input, the puff duration, and the thermal inertia of the coil/wick assembly.

Due to the differing volatilities of the evaporating species, the liquid composition in the heated portion of the wick will tend to become enriched in the less volatile species during a puff, causing the boiling temperature to creep upward in time towards the boiling temperature of the least volatile species present in the liquid. Thus during an individual puff, the vaporization process may switch from an evaporation-limited to boiling-limited regime, and the aerosol composition will change from being enriched in the most volatile species to being enriched in the least volatile species, relative to the parent liquid.

In addition to aerosol generation by vaporization, some liquid may be aerosolized by the formation and bursting of vapor bubbles near the ECIG heating coil surface, as suggested by the crackling sound that is often audible during ECIG operation, and the presence of nonvolatile nicotine salts in ECIG aerosols (El-Hellani et al. 2015).

Thus, predicting aerosol composition, including nicotine, requires treating overlapping heat and mass transfer phenomena, as well as mixture thermodynamics. While in principle these phenomena can be represented in fine detail for each individual ECIG device, the goal of the present work was to develop a model general enough to capture the major features of the relevant processes while minimizing the need for highly detailed product description (e.g., as would be required for computational fluid dynamics simulations).

1.2. Model formulation

In this work, a single control volume (CV) model is formulated to compute the vaporization rate of individual species from the coil-wick assembly. Vapors are emitted by evaporation or boiling from the surface of the CV, and fresh liquid enters through the wick boundaries. The material, thermal, electrical, and enthalpy fluxes through the CV are computed in discrete time steps to determine its instantaneous temperature and composition, which are assumed uniform everywhere. The temperature outside the CV is assumed to be that of the surroundings, , and the liquid composition outside the CV is assumed to be constant and equal to that of the parent liquid.

During puffing, fresh liquid is assumed to enter the cylinder at a rate equal to the instantaneous liquid vaporization rate,, such that the total CV mass remains constant. Fresh liquid is assumed to mix instantaneously with the contents of the CV, and the composition of the liquid in the CV is computed by species conservation. Between puffs, it is assumed that species diffusion is rapid, and that by the start of the next puff the liquid composition in the CV matches that of the parent liquid in the reservoir.

The principle of energy conservation is employed to compute the instantaneous temperature of the CV, with the thermal energy flows accounted for in the model shown in Figure 3. They include the conversion of input electrical energy into thermal energy, the heat loss to the air flowing over the cylinder, the heating and vaporization of liquid flowing through the CV, heat dissipated by conduction to the surroundings, and the heating up of the CV itself. The instantaneous temperature of the cylinder, T (K), is found by numerical integration of the unsteady energy equation:[1]

Figure 3. One-zone model (dashed line). Heat and mass fluxes across the control volume are computed in a time-resolved manner, as is its mean temperature and composition.

where is the electrical power input to the device (W); is the rate of heat transfer by convection to the air flowing over the cylinder (W), where h is the heat transfer coefficient of a cylinder in cross flow (W m⁻² K⁻¹), A is the surface area of the cylinder (m²), and is the ambient air temperature; is the rate of heat transfer by conduction through the wick and electrical leads, where K is the lumped thermal resistance to the surroundings (W K⁻¹); is the latent heat associated with change of phase from liquid to vapor, where is the liquid vaporization rate (kg s⁻¹), and the latent heat of vaporization, , is taken as the mass-weighted average of the latent heats of the liquid constituents (J kg⁻¹); is the energy expended heating the liquid entering the CV, where is the specific heat of the liquid (J kg⁻¹ K⁻¹); C is the effective heat capacity of the liquid-wick system given by (kJ K⁻¹), where V is the volume (m³), and ρ is the density (kg m⁻³), f is an empirical constant used to account for the non-uniform heating of the wick and liquid in the vicinity of the heating coil, and is derived by minimizing error between predictions and observations (see the online supplementary information [SI]) for methodology and sensitivity analyses).

The boiling temperature of the liquid, (K), is evaluated based on the instantaneous composition ignoring trace constituents, assuming ideal solution behavior. See the SI for empirical validation of the ideal solution assumption for PG-VG mixtures.

When the CV temperature is below boiling, , the vaporization rate of species i, is computed assuming evaporative partitioning governed by convective mass transfer:[2] where is the partial pressure of the evaporating species (Pa), is the mole fraction of i in the liquid, and is the vapor pressure of the evaporating species computed at . M is the molar mass of the evaporating species (kg mol⁻¹), R_u is the universal gas constant (J mol⁻¹ K⁻¹), h_m is the mass transfer coefficient (m² s⁻¹) and is calculated according the Churchill and Bernstein correlation (Churchill and Bernstein 1977) for cross flow over a cylinder. The total vaporization rate, , is then calculated as .

If the cylinder temperature reaches the computed , the system is switched to a boiling regime in which the temperature of the cylinder is forced to , and the liquid vaporization rate is then computed by balancing the resulting energy flows using Equation (1[1] ):[3]

In either vaporization or boiling regimes, the nicotine mass flux is computed as the product of the nicotine mass fraction in the liquid and the total liquid vaporization rate. The mass fractions of solvent species (PG, VG, water) in the vapors are computed on the basis of the mole fraction and vapor pressure of each species, in accordance with Raoult's Law. The instantaneous composition of the liquid remaining in the cylinder is computed using a species mass balance given by:[4] where and are the mass and mass fraction of species i in the cylinder, is the mass of liquid in the cylinder (kg), and is the mass fraction of species i in the parent liquid. It should be noted that the fluxes computed using Equation (4[4] ) represent the potential emissions from the device. In reality, some portion of the vaporized material will re-condense on the inner surfaces of the ECIG device, depending on the flow geometry and operating conditions. Thus for clarity, computed fluxes are referred to below as “potential emissions” or “potential yield.” “Yield” refers to the integrated flux (i.e., the total mass) emitted over a given number of puffs.

In summary, the model provides a time-resolved description of the vaporization processes during ECIG operation, accounting for the changing temperature and composition of the liquid in the heated zone resulting from the intermittent powering of the heating coil and the airflow drawn over it. The model utilizes one empirical constant, f, to account for the non-uniform heating of the wick. f is optimized to provide the best correlation between computed potential and measured nicotine yields over a range of liquid compositions, powers, and puff parameters, as described in the SI.

A simplified steady-state model of nicotine flux can be derived by assuming that the cylinder instantaneously attains the liquid boiling temperature when the coil is energized, that the vapor composition is equal to the parent liquid composition and is constant in time, and that heat losses from the cylinder to the surroundings can be neglected. Under these conditions, the vaporization rate during a puff may be computed directly as , where .

2. Materials and methods

To validate model predictions and to explore ECIG emissions over a wide range of operating conditions, yields of nicotine, total particulate matter (TPM), PG, and VG were measured while manipulating puff duration, flow rate, electrical power, and liquid composition (PG/VG ratio, water content, nicotine concentration) using analytical grade liquids and the VaporFi Platinum™ (VF) tank system. Measurements were compared to model predictions for the same variables, under the same conditions. In order to examine the applicability of the model to other ECIG technologies, two additional data sources were used to compare model predictions to measurements. The first source was previously un-reported data from measurements made using a re-fillable unbranded dual-coil atomizer with commercial ECIG liquids. The second source was data reported previously (Talih et al. 2015), using single-coil prefilled cartridges (Vapor4Life CoolCart®, 18 and 36 mg/ml nominal nicotine concentrations). Thus a total of 100 test conditions were examined for the above mentioned yields, as outlined in Table 1.

Table 1. List of conditions tested for nicotine yield (a more detailed list is available in the online supplementary information).

ECIG type/brand	Puff duration (s)	Flow rate (L/min)	Power (W)	PG/VG (vol/vol)	Nicotine concentration (mg/mL)	Water Content (vol%)	Liquid type
Cartomizer/V4L	2, 4, 8	1, 2	3, 8	84/16	9, 16	8	Pre-filled
Dual coil carto/NB	2, 4	1, 2	7, 18	100/0 to 17/83^b	9, 18, 36	0	Commercial
Tank/VF Platinum	1, 2, 4, 8	1, 2	1 to 12^a	100/0 to 0/100^d	8, 18, 36	0 to 20^c	Analytical

^a1, 2, 3, 4, 8, 12; ^b100/0, 74/26, 58/42, 17/83; ^c0, 0.5, 1, 5, 20; ^d100/0, 90/10, 80/20, 70/30, 60/40, 30/70, 20/80, 10/90, 0/100.

In addition to these test conditions, an exploratory investigation was conducted to examine whether bubble fragmentation could be a significant source of aerosol emissions. To investigate this possibility, NaCl was dissolved in ECIG liquid that was then aerosolized. Because NaCl is non-volatile, its appearance in the aerosol phase serves as a marker of aerosol generation through a mechanism other than evaporation/boiling.

Finally, the ideality of the PG/VG system was examined by measuring the boiling points of varying PG/VG ratio solutions and comparing to boiling points computed assuming ideal solution behavior (Section S2 in the SI).

2.1. Materials

VF tanks and refillable unbranded (NB) dual-coil cartomizers were procured from internet vendors in the United States during the period 11/2013 to 01/2015. The electrical resistances of the atomizers were measured before and after each use, using a standard laboratory ohmmeter at 22°C and found to be 2.4 ± 0.12 Ω for the VF and 1.6 ± 0.02 Ω for the NB.

Analytical grade PG (≥99.5%, CAS 57-55-6), VG (99.0–101.0%, CAS 56-81-5), and nicotine (≥99%, CAS number 54-11-5) were procured from Sigma-Aldrich Corporation and used to prepare liquids of varying PG/VG ratio, nicotine concentration, and water concentration for the VF tank measurements. De-ionized water (18.2 MΩ.cm) was obtained from an in-house laboratory water supply system.

For the refillable dual-coil ECIGs, PG, VG and nicotine solutions were procured from an internet vendor, and mixed to produce liquids of varying composition. Concentrations of nicotine (102.02 ± 3.16 mg/mL), and water (PG base: 0.27 ± 0.00 wt%, VG base: 0.10 ± 0.02 wt%) for these solutions were verified using gas chromatography (nicotine) and Karl Fisher technique (water), and were equal to labeled values.

Because the composition of the liquids found in the V4L cartridges used in Talih et al. (2015) was not characterized previously, randomly selected V4L cartridges (N = 3 for each labeled nicotine concentration) that were procured during that study were analyzed. Nicotine concentrations were found to be 8.53 ± 0.71 for the 8 mg/mL labeled cartridges and 15.73 ± 1.21 mg/mL for the 18 mg/mL labeled cartridges. Water content was 7.68 ± 2.56 wt% and PG/VG ratio was 84/16 vol/vol.

Internal standards for nicotine (Hexadecane, ≥99%, CAS 544-76-3) and PG/VG (β-Citronellol, 95%, CAS 106-22-9) quantification were procured from Sigma-Aldrich.

2.2. Aerosol generation and sampling

Aerosols were machine-generated using a custom-designed digital puff production machine at the American University of Beirut (Talih et al. 2015). For each experiment, the mouth end of the ECIG cartridge was connected to a polycarbonate filter holder that contains a Gelman type A/E 47 mm glass fiber filter. Total particulate matter (TPM) was determined by weighing the filter pad and holder before and after each sampling session. Filters were also analyzed for nicotine, propylene glycol, and glycerin content. Because of the large particle mass concentrations of these analytes relative to their saturation vapor concentrations, less than 1% of the mass partitions to the vapor phase and no attempt was made to collect and analyze vapors. That the vapor phase mass was negligible was verified in preliminary measurements by a nicotine mass balance; it was found that nicotine collected on the filters equaled that lost by the ECIG liquid to within experimental error.

2.3. Chemical analysis

Nicotine in the aerosol and liquid was determined by GC-MS analysis of samples extracted in an ethyl-acetate solvent as described in (Talih et al. 2015). An extracted calibration curve with concentrations ranging from 1 to 20 ppm and spiked with the internal standard hexadecane was used to interpret the resulting chromatograms. Spiked filter assays of nicotine in PG, PG/VG (50/50 and 70/30) and VG solutions showed recoveries of >90%.

PG and VG concentrations in the aerosol were measured by dissolving the filter in 4 mL of ethyl acetate and sonicating the solution for 30 min. A 10-fold diluted solution was spiked with the internal standard β-Citronellol and run on GC-FID. The quantification was done against a calibration curve of standard 20–315 ppm of PG or VG in ethyl acetate solutions spiked with the same internal standard β-Citronellol (50 ppm). Recoveries of standard solutions prepared in the lab and spiked on filters were >90%.

2.4. Salt tracer investigation

A stock aqueous solution of NaCl was prepared (320.2 mg/mL) and 150 μL of this solution was diluted to 10 mL using a 70/30 and a 0/100 (vol/vol) PG/VG solution. Aerosols were generated from the prepared liquids under the following conditions: 4W and 11W power, each with 4 and 8 sec puff duration (N = 3 per condition). NaCl analysis in the liquid and aerosol was performed using a Metrohm 850 Ion Chromatogram equipped with A Supp 7–250/4.0 IC column specific for anion separation. Quantification was done using a calibration curve in the range 0.1–20 μg/mL. The percentage of mechanically generated aerosol for a given condition was computed as the ratio of the measured NaCl concentration in the aerosol to that in the liquid.

2.5. Statistical analysis

Average nicotine flux for each condition was determined by dividing the nicotine yield by the total puff duration (number of puffs × puff duration). Linear regression of predicted and measured variables was performed using SPSS version 23.0 (IBM, Armonk, NY, USA). Statistical significance was p < 0.05.

2.6. Numerical method

Equations (1)–(4) were solved simultaneously using a forward Euler method coded in the Matlab computing environment. Time steps of 0.01 s were used; independence of time step size was assured by comparing simulation results derived using a time step of 0.001 s. The initial temperature condition was set to ambient (27°C).

3. Results

3.1. Nicotine yield, flux, and aerosol concentration

Across conditions, nicotine flux ranged from 0.09 to 103.03 μg/s, corresponding to a yield range of 0.01 to 6.65 mg in 15 puffs. Nicotine concentration in the aerosol ranged from 3.84 to 34.29 μg Nic/mg TPM, and was nearly equal to the concentration of the parent liquid across conditions; nicotine yield could be reliably predicted from measurement of the TPM and knowledge of the nicotine concentration of the liquid, as shown in Figure 4.

Figure 4. Nicotine yield versus TPM multiplied by the nicotine mass fraction. r is the correlation coefficient, N is the sample size.

3.2. Effect of power, puff duration, and liquid composition on nicotine flux: Measurements and model optimization

Model predictions of the effects of liquid composition, power, and puff duration on potential nicotine flux are shown in Figure 5, along with the measured values. Model predictions replicated all the trends found with the measurements: increases in flux with power, puff duration, and nicotine concentration, no sensitivity of nicotine flux to air flow rate or water content, and decreases in nicotine flux with increasing VG content.

Figure 5. Effects of power, puff duration, and liquid composition (nicotine concentration, PG/VG ratio, and water content) on normalized nicotine flux. Unless otherwise shown, experimental conditions were: 4 Watts, 4 sec puff duration, 1 L/min flow rate, 10 s interpuff interval, zero water content, and 8 mg/mL nicotine concentration. Lines represent model predictions, while points represent measurements. Error bars: 95% CI. r is the correlation between the measured and predicted nicotine fluxes, N is the sample size.

3.3. Effect of power on aerosol composition for binary liquid solutions: Measurements and predictions

For a given liquid composition and puff duration, at low power (e.g., below 4 Watts) the measured aerosol composition was found to be depleted of the lower volatility species relative to the liquid. For the low power conditions, the lower volatility species (VG) were found in lower concentrations in the aerosol than in the parent liquid. As power was increased, the composition of the aerosol approached that of the parent liquid. These trends were reproduced by the model (r = 0.75, p < .001), as shown in Figure 6.

Figure 6. Effect of power on aerosol composition for binary liquid solutions. PG fraction in the aerosol (wt%) for a 50/50 PG/VG liquid. (Error Bars: 95% CI.) 4 sec puff duration, 10 s interpuff interval and 1 L/min flow rate.

3.4. Measured versus computed nicotine flux across all conditions

Measured nicotine flux was found to be highly correlated to computed potential flux, across all conditions tested (Figure 7, r = 0.85, p < .0001). Computed potential flux was a factor of 2 times the measured yields (see Figure A4 for Bland-Altman plot), indicating that approximately half the material vaporized in the heating zone does not exit the mouthpiece.

Figure 7. Measured versus predicted nicotine flux for three ECIG devices with varying liquid composition, power, and puff duration. (1:2 parity line shown for reference. Error bars: SD.)

3.5. Simplified steady-state model

Measured and computed nicotine fluxes using the simplified steady state model were highly correlated across the range of tested conditions (r = 0.84, p < .0001). Relative to the complete model, the steady-state model predicted equally accurately the effects of electrical power and nicotine concentration on flux. The simplified model erroneously predicted a null effect of puff duration on nicotine flux, and erroneously predicted a significant decrease in nicotine flux with increasing water content. The simplified model also under-predicted by a factor of 3 the sensitivity of nicotine flux to PG/VG ratio at the 4W condition.

3.6. Salt tracer measurements

Salt concentration in the aerosol was not significantly different across the tested conditions (p = 0.48). Across conditions, the concentration of NaCl in the aerosol was equal to 3.98 ± 2.18% of that of the parent liquid. Particle generation by means other than vaporization thus accounted for less than 5% of the aerosol mass emitted from the ECIGs under any condition tested.

4. Discussion

The primary aims of this study were to elucidate over a wide range of plausible use conditions the phenomena that govern ECIG nicotine emissions, and to develop a physics-based mathematical model that can be used to predict those emissions given device design, liquid composition, and operating conditions. Consistent with previous reports (Kosmider et al. 2013; Talih et al. 2015), we found that nicotine flux increases with power output, puff duration, and liquid PG and nicotine concentration. Also consistent with previous reports (Talih et al. 2015), we found that flow rate—and therefore puff volume for a fixed puff duration—does not affect nicotine emissions. We also learned that water content does not impact nicotine flux, that the PG/VG system can be modeled as an ideal solution, and that nicotine yield can be estimated reliably from knowledge of the nicotine concentration in the liquid and the change in mass of the ECIG. Finally, we found that mechanical generation of aerosol particles via bubble burst could account for no more than 5% of the emitted aerosol mass for any condition studied, indicating that a vaporization-centered model is appropriate for predicting ECIG emissions. Limitations of this study include the fact that model predictions were compared to measurements for a small number of ECIG devices, and that all tests were conducted using a steady periodic puffing regimen using rectangular puff velocity waveforms. An additional limitation of the model is the assumption that liquid transport through the wick is always sufficiently rapid as to provide a consistently wet heated zone, i.e., vaporization is never rate limited by the wicking velocity.

The model was able to predict the individual effects of varied liquid composition, electrical power input, and puff topography variables on nicotine flux. Furthermore, over a wide range of use conditions in which these factors varied simultaneously, model predictions of nicotine flux were well correlated to measured nicotine, regardless of ECIG technology (tank, single-coil cartomizer, dual-coil cartomizer) and regardless of whether the liquids were constituted in the lab from analytical standards or were purchased as ready-to-use ECIG liquids. Taken together, these findings indicate that the formulated model reasonably represents the salient physical phenomena that underlie ECIG aerosol production.

4.1. Aerosol composition

A key observation in this regard was the fact that for a given liquid, electrical power affected not only the amount of aerosol emitted, but also its composition. In particular, as shown in Figure 6, at low power the aerosol produced by a 50/50 PG/VG liquid solution is composed almost entirely of PG, the more volatile component of the solution. This result is consistent with evaporative partitioning theory for ideal solutions (i.e., Raoult's Law), which states that the evaporation rate of each species is proportional to the product of its intrinsic vapor pressure and its mole fraction in the liquid phase. Because the vapor pressure of PG is much greater than that of VG, vapors evaporating from a 50/50 PG/VG solution must be composed mainly of PG. However, as power was increased, Figure 6 shows that the composition shifted towards that of the parent liquid, i.e., towards a 50/50 PG/VG composition.

Application of Raoult's law shows that a liquid composed of 7/93 PG/VG is required to produce an aerosol of 50/50 PG/VG composition. In other words, at the higher power conditions, the local liquid composition in vicinity of the coil must differ greatly from the 50/50 PG/VG parent liquid found in the tank reservoir. These observations point to the conclusion that resistance to species transport by diffusion within the wick gives rise to a concentration gradient between the tank and heated zone. As a result, accurate prediction of nicotine flux over a range of ECIG operating conditions requires a model that allows the composition of the liquid near the heating coil to vary in time, even if the parent liquid composition in the tank remains unchanged.

The above dynamics are observable in the model simulation of a single puff at low and high power, shown in Figure 8. For both conditions, at the start of the puff the mass fractions of liquid PG and VG (Figure 8b) correspond to those of a 50/50 PG/VG liquid (by volume). Once the puff starts, the temperature of the heated zone begins to rise (Figure 8a), and evaporation commences. During the early phase of the puff the more volatile PG is driven off the liquid at a far higher rate than is the less volatile VG (Figure 8c), and as a result the vapors are almost exclusively composed of PG (Figure 8b). As the puff proceeds, the temperature continues to rise, causing the vaporization rate to accelerate (Figure 8c), and the liquid composition in the wick to shift towards VG. Eventually, this shift is sufficient to increase the VG presence in the vapors as well. This shift occurs far more rapidly in the high power condition, resulting in a vapor composition that more rapidly approaches that of the parent liquid.

Figure 8. Model predictions of temperature, liquid composition, and vaporization rate for a 50/50 PG/VG liquid at two ECIG power levels (2W, left; 11 W right). 4 s puff duration, 10 s interpuff interval and 1 L/min flow rate.

4.2. Transient and steady state operation

The temperature profile of the high power condition shown in Figure 8a also illustrates key phases in the evolution of a puff: a transient phase during which temperature, local liquid composition, and boiling temperature evolve, and a steady-state phase during which the composition and temperature of the system stabilize. During the transient phase, the vaporization regime may switch from evaporation to boiling if the boiling temperature of the liquid is attained (Figure 8a2). Because the liquid rapidly depletes of PG during boiling, the boiling point, and therefore liquid temperature, rises considerably, and a high liquid vaporization rate is attained during this phase (Figure 8c2), despite the system not yet having attained its maximum temperature. As the PG in the liquid depletes, the system approaches steady state and the vaporization rate decreases from the peak value attained during the transient phase. The combination of the greater liquid temperature and greater VG mole fraction causes the VG to begin evaporating at a rate similar to that of the PG (Figure 8c2). At this point an equilibrium is attained in which the mass flux of the evaporating species is exactly balanced by their replenishment from incoming fresh liquid, and the composition of the liquid in the heated zone—and therefore the boiling point—stabilize, such that a steady-state temperature is also attained.

4.3. Effect of puff duration and thermal inertia on nicotine flux

Because total vaporization rate increases until the system attains steady state, longer puff durations translate to a greater proportion of the time that the system operates in a high vaporization rate mode. That is, one 4-second puff will produce more nicotine emissions than two 2-second puffs. The greater the vaporization rate, in turn, the greater will be the rate at which fresh liquid—and with it, additional nicotine—is wicked into the heated zone. These factors explain the observed trend in nicotine flux versus puff duration (Figure 5b); longer puffs exhibit greater nicotine fluxes, but with diminishing effect due to the fact that a maximum is attained once the system reaches steady state. In the limit, as puff duration approaches infinity, the average flux approaches that of the steady state value. Importantly, the time to reach steady state is a function of the thermal inertia of the system; all else being equal, a more massive coil/wick assembly, for example, will take more time to reach steady state, and therefore for fixed puff duration, will result in a lower average nicotine flux.

4.4. Synthesis: An overarching principle governing nicotine emissions

Because TPM yield is to negligible error equal to the mass of liquid vaporized, Figure 4 can be interpreted as stating that nicotine yield is ultimately determined by how much liquid of a given concentration has been vaporized during a puffing session. Thus, any factor that influences ECIG nicotine emissions must do so by influencing the rate at which fresh liquid, and therefore nicotine, is fed to the heated zone. This simple overarching principle can be invoked to explain any of the observations made in the current study.

Take, for example, the observation that reducing the PG/VG ratio considerably reduced the nicotine flux (Figure 5d). Model simulations showed that with the 100/0 condition, nicotine, whose vapor pressure falls between that of PG and VG, is initially slow to evaporate compared to PG. As a result nicotine initially builds up in the liquid, and eventually begins to evaporate rapidly only once its concentration has sufficiently increased. After some time a steady condition is attained in which the nicotine concentration in the liquid is such that the liquid inflow and vaporization rates balance. Since as a trace constituent nicotine does not alter the gross material transport rate in the wick, the liquid inflow rate is governed strictly by the rate at which PG vaporizes. In other words, the availability of nicotine for vaporization from the heated zone is controlled by the rate at which PG vaporizes.

At the other extreme, a 0/100 PG/VG ratio, the higher temperature attained with VG (due to its higher boiling point relative to PG) causes nicotine to rapidly vaporize and deplete from the heated zone. More nicotine becomes available for vaporization only as additional VG arrives there to make up for what has vaporized. Thus, regardless whether the nicotine is more or less volatile than the solvent, its vaporization rate—i.e., the nicotine flux—is governed by the rate at which the solvent vaporizes. Because of PG's greater volatility and similar enthalpy of vaporization and specific heat to VG, PG-rich liquids vaporize more rapidly than VG-rich liquids do, resulting in a greater nicotine flux for the former, as demonstrated in Figure 5d.

Conversely, we can also understand the null effect of water on nicotine emissions (Figure 5e) as one in which water content must not affect the gross liquid vaporization rate. Although water is more volatile than PG, suggesting a more rapid vaporization rate at a given temperature, it also has a larger specific heat and heat of vaporization than PG, both of which tend to reduce the temperature of the system. However, the effect of greater volatility is offset by the effect of greater heat of vaporization of water. This results in an approximately constant liquid vaporization rate, regardless of the water content of the liquid. As a result, there is also a constant flow of nicotine into the wick across the various water concentration conditions, and therefore constant nicotine emissions.

Similarly, puff velocity did not affect nicotine flux because, as reported in (Talih et al. 2015), while higher air flow rate increases the mass transfer coefficient, it also reduces the temperature of the liquid, resulting in a lower vapor concentration in the boundary layer when the system is in the evaporation operating regime. As a result, the net effect of changing air flow on liquid vaporization rate, and therefore nicotine flux, is null. In the boiling regime, the added convection heat transfer at higher air flow rates represents a negligibly small additional heat loss compared to the latent heat flux, and the vaporization rate—and therefore the nicotine flux—is not significantly affected.

In summary, the physical phenomena that govern bulk liquid vaporization rate from the wick, and therefore the transport of fresh liquid to the heated zone, ultimately determine the availability of nicotine for aerosolization.

5. Conclusions

Nicotine flux can vary widely for a given ECIG device, depending on liquid composition, electrical power, and puff topography. We found that device power, nicotine concentration, PG/VG ratio, and puff duration influenced nicotine flux, while water content and puff velocity did not. The physics-based model described in this study captures the key physical phenomena that govern nicotine flux, providing accurate predictions of the overlapping effects of these factors. It also provides a window into the physics and thermodynamics of ECIG operation.

ECIGs are proliferating and evolving at a greater pace than the knowledge needed to develop policies that address individual and public health concerns. Performance standards are commonly written in consumer product regulations in the United States and elsewhere, and with the appropriate scientific basis, can be developed for ECIGs. To the extent that proposed performance standards will involve nicotine emissions, models like the one described can be used to predict nicotine emissions of any existing or future ECIG product, without recourse to costly testing in an analytical laboratory. These models can also be used to test the effects of proposed regulation on factors that influence nicotine flux.

Funding

This work was supported by the National Institute on Drug Abuse of the National Institutes of Health (grant number P50DA036105) and the Center for Tobacco Products of the U.S. Food and Drug Administration. The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health or the Food and Drug Administration. The authors have no conflict of interest to declare.