Impact of NAPL architecture on interphase mass transfer: A pore network study

Interphase mass transfer in porous media is commonly modeled using Sherwood number expressions that are developed in terms of fluid and porous medium properties averaged over some representative elementary volume (REV). In this work the influence of sub-grid scale properties on interphase mass transfer was investigated using a two-dimensional pore network model. The focus was on assessing the impact of (i) NAPL saturation, (ii) interfacial area (iii) NAPL spatial distribution at the pore scale, (iv) grain size heterogeneity, (v) REV or domain size and (vi) pore scale heterogeneity of the porous media on interphase mass transfer. Variability of both the mass transfer coefficient that explicitly accounts for the interfacial area and the mass transfer coefficient that lumps the interfacial area was examined. It was shown that pore scale NAPL distribution and its orientation relative to the flow direction have significant impact on flow bypassing and the interphase mass transfer coefficient. This results in a complex non-linear relationship between interfacial area and the REV-based interphase mass transfer rate. Hence, explicitly accounting for the interfacial area does not eliminate the uncertainty of the mass transfer coefficient. It was also shown that, even for explicitly defined flow patterns, changing the domain size over which the mass transfer process is defined influences the extent of NAPL bypassing and dilution and, consequently, the interphase mass transfer. It was also demonstrated that the spatial variability of pore scale parameters such as pore throat diameters may result in different rates of interphase mass transfer even for the same pore size distribution index.     

Quantification of Advective Transport Phenomena to Better Understand Dispersion in the Field

Observation of dispersion in field situations has left three issues that may be better understood by applying advective transport phenomena. (1) In some experiments, the longitudinal dispersivity becomes constant with increasing pathlength and in other cases it remains growing. (2) Dispersivities reported from multiple comprehensive observations at a single site differ at similar pathlength in some cases more than a factor two. (3) The observed difference between the plume fronts and plume tails is not represented in the reported parameters. The analytic equations for advective transport phenomena at macroscale of De Lange (2020) describe the thickness of the affected flow-tube and the spread of the plume front and tail. The scale factor defines the size of the averaging domain and so of the initial phase. The new macroscale correlation coefficient relates the growth of the longitudinal dispersivity beyond the initial phase to the aquifer heterogeneity. Using stochastic parameters for the aquifer heterogeneity, the parameters are quantified at 14 field experiments in the United States, Canada and Europe enabling the comparison of calculated and reported final dispersivities. Using the quantified parameters, 146 reported and calculated dispersivities along the traveled paths show a good match. A dispersivity derived from the local plume growth may differ a factor of two from the aquifer-representative value. The growths of plume fronts and tails between two plume stages are assessed in 14 cases and compared to calculated values. Distinctive parameters for the plume front and tail support better understanding of field situations. A user-ready spreadsheet is provided.     

A practical model for mobile, residual, and entrapped NAPL in water-wet porous media

Flow of nonvolatile nonaqueous phase liquid (NAPL) and aqueous phases that account for mobile, entrapped, and residual NAPL in variably saturated water-wet porous media is modeled and compared against results from detailed laboratory experiments. Residual saturation formation in the vadose zone is a process that is often ignored in multifluid flow simulators, which might cause an overestimation of the volume of NAPL that reaches the ground water. Mobile NAPL is defined as being continuous in the pore space and flows under a pressure gradient or gravitational body force. Entrapped NAPL is defined as being occluded by the aqueous phase, occurring as immobile ganglia surrounded by aqueous phase in the pore space and formed when NAPL is replaced by the aqueous phase. Residual NAPL is defined as immobile, nonwater entrapped NAPL that does not drain from the pore spaces and is conceptualized as being either continuous or discontinuous. Free NAPL comprises mobile and residual NAPL. The numerical model is formulated on mass conservation equations for oil and water, transported via NAPL and aqueous phases through variably saturated porous media. To account for phase transitions, a primary variable switching scheme is implemented for the oil-mass conservation equation over three phase conditions: (1) aqueous or aqueous-gas with dissolved oil, (2) aqueous or aqueous-gas with entrapped NAPL, and (3) aqueous or aqueous gas with free NAPL. Two laboratory-scale column experiments are modeled to verify the numerical model. Comparisons between the numerical simulations and experiments demonstrate the necessity to include the residual NAPL formation process in multifluid flow simulators.     

A novel experimental procedure to investigate the biodegradation of NAPL under unsaturated conditions

Soils need to be thoroughly investigated regarding their potential for the natural attenuation of non-aqueous phase liquids (NAPL). Laboratory investigations truly representative of degradation processes in field conditions are difficult to implement for porous media partially saturated with water, NAPL and air. We propose an innovative protocol to investigate degradation processes under steady-state vadose zone conditions. Experiments are carried out in glass columns filled with a sand and, as bacteria source, a soil from a diesel-fuel-polluted site. Water and NAPL (n-hexadecane diluted in heptamethylnonane (HMN)) are added to the porous medium in a two-step procedure using ceramic membranes placed at the bottom of the column. This procedure results, for appropriate experimental conditions, in a uniform distribution of the two fluids (water and NAPL) throughout the column. In a biodegradation experiment non-biodegradable HMN is used to provide NAPL mass, while keeping biodegradable n-hexadecane small enough to monitor its rapid degradation. Biodegradation is followed as a function of time by measuring oxygen consumption, using a respirometer. Degradative activity is controlled by diffusive transfers in the porous network, of oxygen from the gas phase to the water phase and of n-hexadecane from the NAPL phase to the water phase.     

Estimation of water saturation dependence of dispersion in unsaturated porous media: experiments and modelling analysis

In the dispersion theory, a linear relationship has been verified between the coefficient of hydrodynamic dispersion and water velocity, both in saturated and in unsaturated porous media. But for unsaturated soils the variability of flow directions and microscopic velocities can be larger than in saturated soils because of the lower degree of water saturation. This leads to an increased dispersion. Therefore, relationships between water content and relative water velocity fluctuations and water content together with the coefficient of dispersivity in unsaturated porous media respectively have been investigated systematically by displacement experiments in glass beads and coarse-textured sandy soil columns. The breakthrough curves (BTCs) of chloride showed that an increase of solute mixing with a decrease of water content was caused by an increase of flow velocity fluctuations for different pathways. In order to explain the observed tailing effect in unsaturated flow, two mathematical models were used to fit theoretically derived nonlinear functions of water content dependent dispersivities for both porous media. The close agreement between the observed and computed results suggests that the theoretical model of hydrodynamic dispersion can be extended to transport in unsaturated porous media, providing that BTCs of the effluent water are used to estimate representative dispersivity parameters of soils.     

Influence of mass transfer resistance on detection of nonaqueous phase liquids with partitioning tracer tests

Partitioning interwell tracer tests (PITTs) are a relatively new technique for measuring the amount of nonaqueous phase liquid (NAPL) within saturated porous media. In this work we examined the influence of mass transfer limitations on the accuracy of measured NAPL from PITTs. Two mathematical models were used along with laboratory column experiments to explore the influence of tracer partition coefficient, tracer detection limit, and injected tracer mass on NAPL measurements. When dimensionless mass transfer coefficients were small, NAPL measurement errors decreased with decreasing tracer partition coefficient, decreasing tracer detection limit, and increasing injected tracer mass. Extrapolating breakthrough curves exponentially reduced but did not eliminate systematic errors in NAPL measurement. Although transport in a single stream tube was used in the mathematical models and laboratory experiments, the results from this simplified domain were supported by data taken from a three-dimensional computational experiment, where the NAPL resided as large pool. Based on these results, we suggest guidelines for interpreting tracer breakthrough data to ascertain the importance of mass transfer limitations on NAPL measurements.     

A multiphase flow and multispecies reactive transport model for DNAPL-involved Compound Specific Isotope Analysis

This study presents a multiphase flow and multispecies reactive transport model for the simultaneous simulation of NAPL and groundwater flow, dissolution, and reactive transport with isotope fractionation, which can be used for better interpretation of NAPL-involved Compound Specific Isotope Analysis in 3D heterogeneous hydrogeologic systems. The model was verified for NAPL-aqueous phase equilibrium partitioning, aqueous phase multi-chain and multi-component reactive transport, and aqueous phase multi-component transport with isotope fractionation. Several illustrative examples are presented to investigate the effect of DNAPL spill rates, degradation rate constants, and enrichment factors on the temporal and spatial distribution of the isotope signatures of chlorinated aliphatic hydrocarbon groundwater plumes. The results clearly indicate that isotope signatures can be significantly different when considering multiphase flow within the source zone. A series of simulations indicate that degradation and isotope enrichment compete with dissolution to determine the isotope signatures in the source zone: isotope ratios remain the same as those of the source if dissolution dominates the reaction, while heavy isotopes are enriched in reactants along groundwater plume flow paths when degradation becomes dominant. It is also shown that NAPL composition can change from that of the injected source due to the partitioning of components between the aqueous and NAPL phases even when degradation is not allowed in NAPL phase. The three-dimensional simulation is presented to mechanistically illustrate the complexities in determining and interpreting the isotopic signatures with evolving DNAPL source architecture.     

Effect of NAPL Source Morphology on Mass Transfer in the Vadose Zone

The generation of vapor‐phase contaminant plumes within the vadose zone is of interest for contaminated site management. Therefore, it is important to understand vapor sources such as non‐aqueous‐phase liquids (NAPLs) and processes that govern their volatilization. The distribution of NAPL, gas, and water phases within a source zone is expected to influence the rate of volatilization. However, the effect of this distribution morphology on volatilization has not been thoroughly quantified. Because field quantification of NAPL volatilization is often infeasible, a controlled laboratory experiment was conducted in a two‐dimensional tank (28 cm × 15.5 cm × 2.5 cm) with water‐wet sandy media and an emplaced trichloroethylene (TCE) source. The source was emplaced in two configurations to represent morphologies encountered in field settings: (1) NAPL pools directly exposed to the air phase and (2) NAPLs trapped in water‐saturated zones that were occluded from the air phase. Airflow was passed through the tank and effluent concentrations of TCE were quantified. Models were used to analyze results, which indicated that mass transfer from directly exposed NAPL was fast and controlled by advective‐dispersive‐diffusive transport in the gas phase. However, sources occluded by pore water showed strong rate limitations and slower effective mass transfer. This difference is explained by diffusional resistance within the aqueous phase. Results demonstrate that vapor generation rates from a NAPL source will be influenced by the soil water content distribution within the source. The implications of the NAPL morphology on volatilization in the context of a dynamic water table or climate are discussed.     

The influence of wettability on NAPL dissolution fingering

Laboratory experiments and numerical simulations in homogeneous porous media were used to investigate the influence of porous medium wettability on the formation and growth of preferential dissolution pathways, dissolution fingers, during nonaqueous phase liquid (NAPL) dissolution. As the porous medium became increasingly NAPL-wet, dissolution fingers grew wider and slower. This result was observed in physical experiments with 0% and 100% NAPL-wet conditions and confirmed with numerical simulations at these and intermediate wettabilities. A previously derived expression for an upscaled mass transfer rate coefficient that accounts for the growth of dissolution fingers was used to quantify the effect of fingering on overall NAPL removal rates. For the test cases evaluated, NAPL dissolution fingering controlled the overall rate of NAPL dissolution after the dissolution front moved 4 cm in 0% NAPL-wet conditions and 18 cm in 100% NAPL-wet conditions. Thus, even in completely NAPL-wet media dissolution fingering may control the overall rate of NAPL dissolution after relatively short travel distances. The importance of NAPL dissolution fingering in heterogeneous systems with spatially varying NAPL saturations, though, remains an important question for future work.     

Matrix diffusion-derived plume attenuation in fractured bedrock

Matrix diffusion can attenuate the rate of plume migration in fractured bedrock relative to the rate of ground water flow for both conservative and nonconservative solutes of interest. In a system of parallel, equally spaced constant aperture fractures subject to steady-state ground water flow and an infinite source width, the degree of plume attenuation increases with time and travel distance, eventually reaching an asymptotic level. The asymptotic degree of plume attenuation in the absence of degradation can be predicted by a plume attenuation factor, beta, which is readily estimated as R' (phi(m)/phi(f)), where R' is the retardation factor in the matrix, phi(m) is the matrix porosity, and phi(f) is the fracture porosity. This dual-porosity relationship can also be thought of as the ratio of primary to secondary porosity. Beta represents the rate of ground water flow in fractures relative to the rate of plume advance. For the conditions examined in this study, beta increases with greater matrix porosity, greater matrix fraction organic carbon, larger fracture spacing, and smaller fracture aperture. These concepts are illustrated using a case study where dense nonaqueous phase liquid in fractured sandstone produced a dissolved-phase trichloroethylene (TCE) plume approximately 300 m in length. Transport parameters such as matrix porosity, fracture porosity, hydraulic gradient, and the matrix retardation factor were characterized at the site through field investigations. In the fractured sandstone bedrock examined in this study, the asymptotic plume attenuation factors (beta values) for conservative and nonconservative solutes (i.e., chloride and TCE) were predicted to be approximately 800 and 12,210, respectively. Quantitative analyses demonstrate that a porous media (single-porosity) solute transport model is not appropriate for simulating contaminant transport in fractured sandstone where matrix diffusion occurs. Rather, simulations need to be conducted with either a discrete fracture model that explicitly incorporates matrix diffusion, or a dual-continuum model that accounts for mass transfer between mobile and immobile zones. Simulations also demonstrate that back diffusion from the matrix to fractures will likely be the time-limiting factor in reaching ground water cleanup goals in some fractured bedrock environments.     

A field study of the coupled effects of aquifer stratification,fluid density,and groundwater fluctuations on dispersivity assessments

A number of experimental studies have tackled the issue of solute transport parameter assessments either in the laboratory or in the field. But yet, the behavior of a plume in the field under density driven forces, is not well known due to possible development of instabilities. Some field tracer tests on the fate of plumes denser than native groundwater such as those encountered under waste disposal facilities, have pointed out the processes of sinking and splitting at the early stage of migration. The process of dispersion was widely investigated, but the range of dispersivity values obtained from either experimental tests, or numerical and theoretical calculations is still very large, even for the same type of aquifers. These discrepancies were considered to be essentially caused by soil heterogeneities and scale effects. In the meantime, studies on the influence of sinking and fingering have remained more scarce. The objective of the work is to analyze how transport parameters such as dispersivities can be affected by unstable conditions, which lead to plume sinking and fingering. A series of tracer tests were carried out to study under natural conditions, the transport of a dense chloride solution injected in a shallow two-layered aquifer. Two types of experiments were performed: in the first type, source injection was such that the plume could travel downward from one layer to the other of higher pore velocity, and in the second one, the migration took place only in the faster layer. The results suggest some new insights in the processes occurring at the early stages of a dense plume migration moving in a stratified aquifer under groundwater fluctuations, which can be summarized through the following points: (i) Above a stability criterion threshold, a fingering process and a multi modal plume transport take place, but local dispersivities can be cautiously derived, using breakthrough curves matching. (ii) When water table is subject to some cycling or rising, the plume can be significantly distorted in the transverse direction, leading to unusual values of the ratio between longitudinal and transverse dispersivities. (iii) Under stable conditions, for example in the case of straightforward injection in the faster aquifer layer, longitudinal dispersivity is greater than the transverse component as usually encountered, and the obtained transport parameters are closed to macro dispersivity values, which reach their asymptotic limit at very short distances. (iv) The classical scale effect about the varying dispersivity at short distances could be a process mainly due to the distance required for a plume stabilization.     

Lessons Learned from 25 Years of Research at the MADE Site

Field studies at well‐instrumented research sites have provided extensive data sets and important insights essential for development and testing of transport theories and mathematical models. This paper provides an overview of over 25 years of research and lessons learned at one of such field research sites on the Columbus Air Force Base in Mississippi, commonly known as the Macrodispersion Experiment (MADE) site. Since the mid‐1980s, field data from the MADE site have been used extensively by researchers around the world to explore complex contaminant transport phenomena in highly heterogeneous porous media. Results from field investigations and modeling analyses suggested that connected networks of small‐scale preferential flow paths and relative flow barriers exert dominant control on solute transport processes. The classical advection‐dispersion model was shown to inadequately represent plume‐scale transport, while the dual‐domain mass transfer model was found to reproduce the primary observed plume characteristics. The MADE site has served as a valuable natural observatory for contaminant transport studies where new observations have led to better understanding and improved models have sprung out analysis of new data.     

Fractal generation of surface area of porous media

Many natural porous geological rock formations, as well as engineered porous structures, have fractal properties, i.e., they are self-similar over several length scales. While there have been many experimental and theoretical studies on how to quantify a fractal porous medium and on how to determine its fractal dimension, the numerical generation of a fractal pore structure with predefined statistical and scaling properties is somewhat scarcer. In the present paper a new numerical method for generating a three-dimensional porous medium with any desired probability density function (PDF) and autocorrelation function (ACF) is presented. The well-known Turning Bands Method (TBM) is modified to generate three-dimensional synthetic isotropic and anisotropic porous media with a Gaussian PDF and exponential-decay ACF. Porous media with other PDF's and ACF's are constructed with a nonlinear, iterative PDF and ACF transformation, whereby the arbitrary PDF is converted to an equivalent Gaussian PDF which is then simulated with the classical TBM. Employing a new method for the estimation of the surface area for a given porosity, the fractal dimensions of the surface area of the synthetic porous media generated in this way are then measured by classical fractal perimeter/area relationships. Different 3D porous media are simulated by varying the porosity and the correlation structure of the random field. The performance of the simulations is evaluated by checking the ensemble statistics, the mean, variance and ACF of the simulated random field. For a porous medium with Gaussian PDF, an average fractal dimension of approximately 2.76 is obtained which is in the range of values of actually measured fractal dimensions of molecular surfaces. For a porous medium with a non-Gaussian quadratic PDF the calculated fractal dimension appears to be consistently higher and averages 2.82. The results also show that the fractal dimension is neither strongly dependent of the porosity nor of the degree of anisotropy assumed.     

Fractal generation of surface area of porous media

Many natural porous geological rock formations, as well as engineered porous structures, have fractal properties, i.e., they are self-similar over several length scales. While there have been many experimental and theoretical studies on how to quantify a fractal porous medium and on how to determine its fractal dimension, the numerical generation of a fractal pore structure with predefined statistical and scaling properties is somewhat scarcer. In the present paper a new numerical method for generating a three-dimensional porous medium with any desired probability density function (PDF) and autocorrelation function (ACF) is presented. The well-known Turning Bands Method (TBM) is modified to generate three-dimensional synthetic isotropic and anisotropic porous media with a Gaussian PDF and exponential-decay ACF. Porous media with other PDF's and ACF's are constructed with a nonlinear, iterative PDF and ACF transformation, whereby the arbitrary PDF is converted to an equivalent Gaussian PDF which is then simulated with the classical TBM. Employing a new method for the estimation of the surface area for a given porosity, the fractal dimensions of the surface area of the synthetic porous media generated in this way are then measured by classical fractal perimeter/area relationships. Different 3D porous media are simulated by varying the porosity and the correlation structure of the random field. The performance of the simulations is evaluated by checking the ensemble statistics, the mean, variance and ACF of the simulated random field. For a porous medium with Gaussian PDF, an average fractal dimension of approximately 2.76 is obtained which is in the range of values of actually measured fractal dimensions of molecular surfaces. For a porous medium with a non-Gaussian quadratic PDF the calculated fractal dimension appears to be consistently higher and averages 2.82. The results also show that the fractal dimension is neither strongly dependent of the porosity nor of the degree of anisotropy assumed.     

A review of visualization techniques of biocolloid transport processes at the pore scale under saturated and unsaturated conditions

Field and column studies of biocolloid transport in porous media have yielded a large body of information, used to design treatment systems, protect water supplies and assess the risk of pathogen contamination. However, the inherent “black-box” approach of these larger scales has resulted in generalizations that sometimes prove inaccurate. Over the past 10–15 years, pore scale visualization techniques have improved substantially, allowing the study of biocolloid transport in saturated and unsaturated porous media at a level that provides a very clear understanding of the processes that govern biocolloid movement. For example, it is now understood that the reduction in pathways for biocolloids as a function of their size leads to earlier breakthrough. Interception of biocolloids by the porous media used to be considered independent of fluid flow velocity, but recent work indicates that there is a relationship between them. The existence of almost stagnant pore water regions within a porous medium can lead to storage of biocolloids, but this process is strongly colloid-size dependent, since larger biocolloids are focused along the central streamlines in the flowing fluid. Interfaces, such as the air–water interface, the soil–water interface and the soil–water–air interface, play a major role in attachment and detachment, with significant implications for risk assessment and system design. Important research questions related to the pore-scale factors that control attachment and detachment are key to furthering our understanding of the transport of biocolloids in porous media.     

Dependence of lumped mass transfer coefficient on scale and reactions kinetics for biologically enhanced NAPL dissolution

Residual dense non-aqueous liquids (NAPLs) in aquifers constitute a great challenge for groundwater cleanup. Active engineered treatment of regions that contain residual NAPLs is often required to shorten the long-term impact of NAPLs on groundwater quality. Enhanced residual NAPL cleanup can be achieved by promoting biodegradation of NAPL components in the aqueous phase, thereby increasing contaminant fluxes from the NAPL phase. Reaction-enhanced NAPL dissolution is often mathematically simulated under the assumption that lumped mass transfer coefficients, used to describe the dissolution behavior of the NAPL phase, are independent of the reactions. However, this assumption is not warranted because reactions occurring near the water–NAPL interface can reduce characteristic mass transfer lengths, which tend to enhance mass transfer over the no-reaction case.In this study, we mathematically investigated the connections between lumped mass transfer coefficients and reaction kinetics over an idealized residual NAPL domain. Since mass transfer is frequently a scale-dependent process, we also examined the influence of system extent on mass transfer coefficients. For our idealized domain with an assumed first-order decay reaction, the results show that lumped mass transfer coefficients depend on reaction kinetics and system scale. The mass transfer coefficient derived from the non-reactive case cannot properly represent the mass transfer process under the reactive conditions. When the advection time scale is long in comparison to the transverse dispersion time scale in the system, a fast reaction can increase significantly the lumped mass transfer coefficient. The mass transfer coefficient used for simulation was also found to be affected by the nature of the numerical scheme used.     

Fluid dispersion effects on density-driven thermohaline flow and transport in porous media

This study introduces the dispersive fluid flux of total fluid mass to the density-driven flow equation to improve thermohaline modeling of salt and heat transports in porous media. The dispersive fluid flux in the flow equation is derived to account for an additional fluid flux driven by the density gradient and mechanical dispersion. The coupled flow, salt transport and heat transport governing equations are numerically solved by a fully implicit finite difference method to investigate solution changes due to the dispersive fluid flux. The numerical solutions are verified by the Henry problem and the thermal Elder problem under a moderate density effect and by the brine Elder problem under a strong density effect. It is found that increment of the maximum ratio of the dispersive fluid flux to the advective fluid flux results in increasing dispersivity for the Henry problem and the brine Elder problem. The effects of the dispersive fluid flux on salt and heat transports under high density differences and high dispersivities are more noticeable than under low density differences and low dispersivities. Values of quantitative indicators such as the Nusselt number, mass flux, salt mass stored and maximum penetration depth in the brine Elder problem show noticeable changes by the dispersive fluid flux. In the thermohaline Elder problem, the dispersive fluid flux shows a considerable effect on the shape and the number of developed fingers and makes either an upwelling or a downwelling flow in the center of the domain. In conclusion, for the general case that involves strong density-driven flow and transport modeling in porous media, the dispersive fluid flux should be considered in the flow equation.     

Magnetic Resonance Imaging of Nonaqueous Phase Liquid Distribution in Unconsolidated Soils Using Discriminatory Freezing

Organic contaminants present as nonaqueous phase liquids (NAPLs) in the subsurface often pose a long-term risk to human health and the environment. Investigating the distribution of NAPLs in porous media remains a major challenge in risk assessment and management of contaminated sites. Conventional soil coring and monitoring wells have been widely used over past decades as the primary means of subsurface investigation to determine NAPL extent. Known limitations of conventional approaches have led us to explore an alternative or a complementary technique to provide high-quality information of NAPL source zone architecture. This work advances an imaging tool for a variety of organic NAPL contaminants in unconsolidated soils through magnetic resonance imaging (MRI) of frozen cores. Using trichloroethylene (TCE) and o-xylene as model species, we illustrate that discriminatory freezing of water, while keeping the NAPL in a liquid state, enables high-resolution qualitative delineation of NAPL distribution within porous media. This novel approach may help improve site conceptual models and consequentially lead to highly tailored, more efficient remedial measures.     

Unifying NAPL Drawdown and Transmissivity Testing in Unconfined,Confined, Perched,and Fractured Settings using the Z-Factor and MH Principles

Light nonaqueous phase liquid (LNAPL) flow in in fractured rock is governed by the same physics as porous media, but LNAPL discharge to a well from fractured rock is subject to the unique geometry of the fractures within the rock and the degree of interconnectivity between the factures. Previous conceptualization and definition of drawdown of nonaqueous phase liquids (NAPL) has employed a single drawdown value to represent the entire vertical interval of mobile NAPL. Application of the single drawdown model may result in erroneous calculation of NAPL transmissivity in fractured rock settings. This work illustrates how drawdown in multiphase systems can be variable over the vertical interval of mobile NAPL. In settings with discrete fracture networks, it is clear that consistently applying a single drawdown value will not accurately represent the pressure gradients. This work presents the multiphase head (MH) model, which is proposed as a comprehensive methodology for evaluating NAPL drawdown in fractured rock, and unconsolidated porous media. The MH model utilizes fluid statics and physical principles to accurately represent pressure differences in the formation and convert those into NAPL drawdown for discrete elevations. This first principles approach to describing how drawdown varies with NAPL-production zone elevations and fluid levels, resulting in a more accurate representation of discharge vs. fluid elevation behavior. Application of the MH model to various scenarios has identified that dissimilar scenarios can represent similar behavior during recovery from a NAPL removal event or baildown test. The resulting understanding improves the selection of representative portions of baildown test data to use in NAPL transmissivity analysis. Proper conceptualization of drawdown in bedrock identifies an alternate analysis method, the Z-factor, to estimate NAPL transmissivity. The resulting drawdown calculations and transmissivity analysis method result in a comprehensive approach to calculating NAPL transmissivity in both bedrock and unconsolidated porous media.