Unified Theory of Ultimate Hydrocarbon Recovery for Primary and Cyclic Injection Processes in Ultratight Reservoirs

This paper presents a simple method to estimate ultimate recovery factors (URF) of ultratight reservoirs based on equilibrium by diffusion in which URF is only a function of changes in hydrocarbon density between initial and final states. URF is defined at infinite time and therefore does not depend on the transient behavior. Although URF may not be achievable during the life-cycle of the field development and production, it provides valuable insights on the role of phase behavior. We show that equilibrium phase behavior defines the absolute upper-bound for URF during primary production and explains the poor recovery from shale oil reservoirs compared to the high recovery factor in shale gas reservoirs in a unifying way. Further, we quantify how injected solvent compositions (CH4, CO2, N2, and C2H6) during huff’n’puff enhanced oil recovery (EOR) improve recovery based on density reduction and compositional dilution, and show that the largest percentage increase in recovery occurs for heavier oils. Our calculations provide a practical means to define the URF from primary production as a function of reservoir fluid composition, temperature, and pressure drawdown. In addition, our calculations articulate incremental URF (IURF) of solvent huff‘n’puff based on net solvent transfer into ultratight rock, which is a key design consideration. The results illustrate that solvent transfer dilutes the hydrocarbons in place, thus maximizing long-term hydrocarbon recovery. Net mass transfer can be improved by enhancing the diffusion of solvent into the matrix based on the huff‘n’puff design parameters including solvent composition, drawdown pressure, and the net amount of solvent injected based on optimal frequency and cycle duration.

huff 'n'puff (HnP) design parameters (injection pressure, solvent composition, solvent amount injected, cycle number and duration) for a given reservoir composition. Furthermore, knowing how to increase the ultimate recovery factors could yield economic increases in the transient recovery. There needs to be a clear framework to approach reservoir development as a unified process spanning primary production and solvent-based EOR/EGR for different fluid compositions.
Miscible gas injection based on the multi-contact-miscibility (MCM) process requires knowledge of which tie line controls miscibility (that is, when a tie line in the two-phase region obtains zero length first). Miscibility attained in this manner results from interaction of advective transport (flow governed by relative permeability) and phase behavior such that the forward contact (oil tie line), backward contact (gas tie line), or a cross-over tie line controls miscibility 17,18 . In an ultratight reservoir matrix, however, miscibility becomes a function of solubility because advection-dominated flow conditions required for the MCM process are no longer relevant 6 . Instead, transport within a tight matrix is anticipated to be a diffusion-dominated process 6,19 . Concentration gradients within the matrix near the matrix-fracture interface control mass transport in these reservoirs for the primary production and HnP process due to the significant contact area created during hydraulic fracturing between the matrix and fracture network.
A unified theory applicable to primary recovery for any fluid type and for any cyclic injection sequence is needed to quantify practical expectations for the maximum ultimate recovery factor achievable. In this paper, we first develop such a unified theory based on diffusion principles, where equilibrium is ultimately achieved between an initial and final state. We evaluate ultimate mass recovery factors for the n-alkane series as a compositional proxy for a range in unconventional reservoir compositions. Next, we compare the potential use of four different solvents (CH 4 , CO 2 , C 2 H 6 , and N 2 ) to enhance ultimate recovery factor for various reservoir compositions, reservoir conditions, and design parameters (amount of solvent injected and producing pressure). We then evaluate recovery (primary and HnP) using several Eagle Ford compositions. Finally, we discuss the theoretical and practical implications of these results for a multiple-cycle HnP process.

Theory and Methods
In this section, we first develop the underlying theory of ultimate recovery for the primary and HnP processes. We then present a sensitivity analysis using a compositional proxy that captures a range of reservoir fluid types from liquid-like to vapor-like.
Theory. Diffusion should equilibrate overall compositions if given sufficient time. Thus, ultimate hydrocarbon mass recovery is strictly a function of the equilibrium hydrocarbon pressure-volume-temperature (PVT) behavior (phase densities, compositions, and saturation) observed at initial (t = 0) and final (t → ∞) reservoir conditions. In a reservoir system with n c components and n p fluid phases, the ultimate mass recovery factor for each component, . Equation (1) allows for any number of phases to be present in the reservoir at initial or final conditions, but assumes no-flux outer boundaries, constant pore volume, spatially uniform initial pressure and temperature distributions, and no gravity effects. (1) is a general equation for URF that applies to primary and/or HnP processes. Equation (1) shows that ultimate recovery for hydrocarbon components will be improved by any process capable of reducing C i,F below the initial reservoir concentration 19 .

Recovery processes. Equation
Primary production. Primary recovery results mainly from the density reduction (fluid expansion) that occurs in the reservoir as pressure declines. Thus, increasing the amount of pressure drawdown (p I − p wf ) during primary production will increase the ultimate recovery factor (attained after an infinite period of production). Equation (1) is used to estimate the ultimate mass recovery [fraction of original component i mass] following primary production, , , where subscript P replaces F to denote "primary". C i,P [lbm/ft 3 ] is the mass concentration of component i at equilibrium conditions following an infinitely-long period of primary production (T  p  z t  , , , is the overall equilibrium composition in the reservoir following primary production. → z P will be constant through the primary recovery process as long as all components can move coherently through the nm-sized pores. However, filtration 20 , size exclusion 21 , and adsorption 21 effects may selectively hinder component transport within small pores or those lined with organic matter. Chromatographic separation 22 could also be important in tight/shale gas reservoirs. Equation (1) can be simplified when the overall composition of the fluid is constant during the primary recovery process and the pressure does not go below the bubble point. Confinement tends to suppress bubble point pressure within nanopores 23,24 . In this case, the primary ultimate recovery factor, URF P , is controlled by the initial (ρ I ) and final (ρ F ) fluid densities within the reservoir, where URF P = 1 − ρ I /ρ F . www.nature.com/scientificreports www.nature.com/scientificreports/ Huff 'n'puff (cyclic solvent) process. A second way to increase URF is the use of a Huff 'n'Puff (HnP) process to reduce C i by decreasing ρ j terms (fluid expansion) and decreasing ω ij (selective dilution of hydrocarbon components) in the reservoir. Solvent-hydrocarbon mixing during the huff and soak causes a local change in fluid composition in the fracture network and adjacent matrix blocks contacting the fracture. We consider two conceptual limits for solvent-hydrocarbon mixing. The lower limit, which we refer to as the early-time limit (ETL), considers perfect solvent mixing throughout the wellbore and fractures by a combination of advection (Darcy's law) and diffusion. Mixing in the early-time volume is not instantaneous, but the ETL assumes no mass-transfer between ultratight matrix and the early-time volume owing to the large conductivity contrast. The upper limit, which we refer to as the late-time limit (LTL), considers perfect solvent mixing throughout the entire reservoir. After a long period of time (t → ∞), diffusion will bring the ultratight matrix into equilibrium with the early-time volume. Figure 1 schematically illustrates the ETL and LTL tank models.
Equation (1) is used to determine the URF's for the HnP process. The fluid composition used in Eq. (1) is assumed to be the equilibrium composition obtained from material balance at the end of the huff and soak period, which is determined by the mass of solvent injected during the huff and the total mass of hydrocarbon remaining in the reservoir after primary. The final densities and compositions are then used in Eq. (1)  = − , where subscript H replaces F to denote "huff ".
In the case of multiple HnP cycles, one must account for the mass of fluid and components (a mixture of hydrocarbon and solvent) produced back during each puff cycle to determine the net amount of solvent (n s [lbmol]) and remaining hydrocarbon (n hc [lbmol]) in the reservoir, We assume equilibrium conditions and compute recovery using compositions from a material balance as a function of the net mole fraction of solvent in the reservoir, z s , given as z s = n s /(n s + n hc ). Values for z s are physically restricted between 0 and 1 [mole fraction]. These limiting values span the continuum between zero solvent injection (z s = 0, primary only) and perfect hydrocarbon dilution/replacement with solvent (z s = 1) in the reservoir. The practical upper-bound for z s from HnP is much less than unity because operating constraints restrict how much solvent can be injected during a given huff cycle. Thus, multiple cycles are typically needed to increase z s , although incremental increases in z s from each additional cycle likely become less. The importance of z s to recovery argues also to increase the injection pressure to its maximum field value.
The values used in Eq. (2) come from either production data or analytical/numerical simulation. Therefore, the actual masses of components produced and injected during each HnP cycle at a well(s) could also be used to directly calculate the expected URF after all HnP cycles are complete. In this paper, we first consider only one HnP cycle to illustrate the fundamentals underlying the recovery process. In the discussion, we extend observations and findings to multiple cycle HnP.
By design, solvent IURF is maximized in the LTL and minimized in the ETL. The performance of a real HnP cycle is expected to fall in between these two limits. We quantify HnP performance based on incremental ultimate recovery factor IURF i , [fraction of original component i mass] because EOR and/or EGR operations are ) is injected into V wf and subsequently produced back in the huff with no mass-transfer between V wf and the ultratight matrix. In the LTL, perfect mixing is achieved (z sf = z sm = z s ) by diffusion during the soaking period (as t → ∞) before starting the huff. n sf , n sm , and n s are the amounts of solvent [lbmol] present in the V wf , ultratight matrix, and the entire reservoir, respectively, and z sf , z sm , and z s are the corresponding solvent mole fractions.
www.nature.com/scientificreports www.nature.com/scientificreports/ judged according to the economic value of incremental hydrocarbon production compared with the reference case (recovery without EOR/EGR). Accordingly, the solvent IURF corresponding to the ETL (IURF i,ETL ) is calculated by,  3 ] is the pore volume of the entire reservoir.
Lumping the wellbore and fractures together into a cup-mixed tank with pore volume V wf [ft 3 ] is reasonable for ultratight reservoirs that act on completely different time scales. Highly permeable layers (e.g. carrier beds) connected to the fractures could, perhaps, also be lumped into V wf . One can then calculate the minimum injected solvent (Δn s , [lbmol]) required to pressurize V wf to a target bottom-hole injection pressure (p H [psia]) using a compositional volume balance method 25,26 , is the cup-mixed overall composition in the fractures, V fluid [ft 3 ] is the water, oil, and gas phase volumes at the specified pressure, temperature, and composition.
Counter-diffusion occurs during the huff and soak because solvent diffuses from fractures into matrix while hydrocarbon simultaneously diffuses into the fractures 6,19,26 . However, in this study, we assume a constant pore volume and ignore mass transfer between V wf and ultratight matrix when using Eq. (4). Furthermore, VBE was implemented when calculating the amount of injected solvent required to pressurize the wellbore-fracture volume to a target pore pressure in the absence of diffusive mass transfer between matrix and V wf . The wellbore-fracture pore volume is fixed when pressure is specified. However, the volume of the hydrocarbon-solvent mixture depends on the mass of injected solvent. The amount of injected solvent is iteratively improved via Newton-Raphson with respect to VBE until volume balance is achieved, i.e., | In other words, the amount injected was iteratively improved using pore volume minus fluid volume as the objective function to minimize in a Newton-Raphson update based on n n VBE n VBE n where VBE′ is the derivative of the VBE with respect to the moles of injected solvent (Δn s ). We conducted a sensitivity study and found that relative error asymptotically stabilizes for any convergence criterion less than 10 −4 . We chose 10 −6 to be conservative and align with typical reservoir simulation convergence criteria.
The solvent IURF corresponding to the LTL (IURF i ) is calculated by, are the only terms that can be modified in field applications because C i,I is fixed. As discussed earlier, C i,P is decreased using pressure-drawdown. C i H LTL , is decreased further during HnP by changing the composition of fluids in the reservoir (compositional dilution). Thus, solvent composition, amount of solvent, and pressure drawdown (initial reservoir pressure minus wellbore pressure) are the three adjustable factors explored in this paper. The difference between Eqs (3) and (5) is the pore volume fraction (V wf /PV Total ) multiplier and the composition vector used to compute mass concentrations in the huff. We assume V wf /PV Total = 0.01 [−] based on flowback analysis 27 .
Compositional proxy system. To demonstrate how the proposed unified theory can explain ultimate hydrocarbon recovery, we use a normal alkane (n-alkane) series from C 1 to n-C 25 as a compositional proxy to continuously represent a range of real reservoir compositions from natural gas (light alkanes) to oil (heavy alkanes). This is because real oils contain a complex mixture of paraffins, naphthenes, and aromatics making it difficult to draw any conclusions without simplifying assumptions. We chose n-C 25 as the upper end of the series because production becomes minimal for long hydrocarbons through the small pore-throats in shale reservoirs. This conclusion is supported by comparing the composition of wellstream samples versus those obtained via solvent extraction of crushed core materials 28 .
Experimental property values do not exist for pure n-paraffins heavier than approximately n-C 25 29 . We supplemented experimental data 30 using the Magoulas-Tassios 31 critical property correlations (T c , p c , ω), which are suitably accurate 32 for the n-alkanes we investigate. We used Whitson's best-fit 33 to Yarborough's correlation 34 to obtain standard condition specific gravity (γ [−]) values for liquid n-alkanes and compute Péneloux volume shift 35 parameters for the Peng-Robinson 36 equation of state (PR EoS). Binary interaction parameters were from Whitson 29 for CO 2 -and N 2 -hydrocarbon pairs and the Chueh-Prausnitz equation 37 for C 1 -C 7 + pairs. All other hydrocarbon-hydrocarbon pairs were zero. Figure 2 shows these input parameters as a function of carbon number (CN). This approach to the input parameters gives consistent and continuous values. Experiments 30 show that critical pressures monotonically decrease for n-alkanes as CN increases (Fig. 2). Ethane is the exception because a lack of compliance in the single carbon-carbon bond prohibits ethane from packing as efficiently as the heavier alkanes containing multiple carbon-carbon bonds 38

Results
In this section, we evaluate changes in primary ultimate recovery factor (URF P ) as functions of reservoir composition (n-alkane), reservoir pressure, drawdown pressure, and reservoir temperature for the compositional proxy system as well as Eagle Ford fluid. Next, we repeat the analysis for HnP where we include the effects of solvent composition and the amount of injected solvent.
Compositional proxy system. Figure 3 shows the effect of temperature on primary ultimate recovery factor (URF P ) versus n-alkane CN as p wf is lowered from p wf = p i − 500 = 7500 to p wf = 500 psia at T res = 150, 200,  www.nature.com/scientificreports www.nature.com/scientificreports/ 250, and 300°F. Contours were generated by successively solving Eq. (1) at every integer (pure n-alkane) CN and 5 psi increment for p wf .
Ultimate recovery dramatically improves whenever an initially liquid-like (L) composition becomes more vapor-like (V) or enters the L-V region (vapor pressure curve if pure-component) at final conditions (p wf , T res ). Closely-spaced near vertical URF P contours denote the transition between vapor-to liquid-like compositions as n-alkane CN increases at fixed p wf . The progressive flattening of URF P contours in Fig. 3 confirms the transition from V to L with increasing CN. Closely-spaced horizontal URF P contours correspond to crossing a pure-alkane vapor pressure curve at fixed CN. Due to higher compressibility in vapor-like phases versus liquid-like phases, the calculated URF P values (which are based on density reduction) are in agreement with our qualitative expectations from compressibility trends.
URF P values systematically improve for all n-alkanes (C 1 through n-C 25 ) with increasing pressure drawdown (p I − p wf ) and monotonically decrease as n-alkanes become heavier (larger CN) at any T res value (Fig. 3) within the pore pressure range of interest. Figure 4 illustrates the variation in URF P with respect to T res for a few representative compositions (C 1 for gas, C 2 , C 3 , and n-C 4 for gas-condensate, n-C 10 for volatile oil, and n-C 20 for black oil) at constant p I = 8000 and p wf = 1000 psia.
The results in Fig. 4 show that URF P systematically improves with increasing temperature and decreasing CN (recall MW ≈ 14CN + 2). In addition, URF P is more sensitive for C 2 , C 3 , and n-C 4 with respect to temperature than the other compositions, a result shown as well in Fig. 3.
Next, we compare the incremental ultimate recovery factor (IURF) using C 1 , C 2 , CO 2 , and N 2 as the injected solvent at different T res (150 and 250°F) and p wf (2000 and 1000 psia) values while keeping the initial reservoir pressure constant at 8000 psia (Fig. 5). IURF is calculated using the late-time limit (LTL) for solvent mixing and the assumption that primary and puff periods have the same p wf . The axes correspond to the amount of solvent retained in the reservoir with perfect mixing (z s ) versus pure n-alkane CN. For practical considerations, we restricted z s to a maximum of 20 mole % since perfect hydrocarbon dilution/replacement (z s = 100 mole %) will never be reached.
Based on our previous discussion of recovery theory, the solvent IURF analysis in Fig. 5 captures the recovery behavior as the solvent mole fraction approaches its conceptual limits. Solvent IURF is zero for z s = 0 mole % (primary only) and approaches large finite values as z s tends to 100 mole %. Solvent IURF must be zero when solvent and n-alkane compositions are identical, which is the case for the C 2 (asymptote of 0.0% for CN = 2) and C 1 (bounding asymptote of 0.0% for CN = 1) cases. Closely-spaced horizontal IURF contours in Fig. 5b,d reveal the pronounced vaporization effect of nitrogen solvent on n-alkanes that would ordinarily be liquid-like at the given conditions (p wf , T res ) but are instead vapor-like (or a two-phase mixture) after the introduction of solvent. Propane (CN = 3) broadly experiences the greatest absolute improvement to recovery factor at any given z s and solvent composition based on finding

s I wf res
∂ ∂ = from graphical analysis (Fig. 5). However, this does not mean that reservoirs with propane-like compositions are "optimal" candidates for solvent HnP.
A practical question is what solvent will be the "best" to inject into a reservoir and how much solvent should be injected. Solvent IURF performance is a function of pressure, temperature, original fluid composition, and the amount and type of solvent injected.
In the context of EOR/EGR, a design is considered "optimal" when it maximizes project economics 17 . On an equal energy content basis, liquid hydrocarbons are more commercially valuable than gaseous hydrocarbons, something which is referred to as the "liquid premium". Primary recovery factors are very large in the first place for light n-alkanes (Figs 3 and 4), which makes HNP less commercially advantageous. In contrast, even a very modest solvent IURF (~1-2% OOIP Mass) would constitute an enormous relative improvement over primary alone (URF P ~3% OOIP Mass) in a heavy oil (Figs 3 and 5). Therefore, the (IURF/URF P ) ratio may provide a quick screening tool for a more formal economic-based assessment.
In a single-phase mixture, the ρω i product yields component mass concentration. The introduction of solvent initiates two competing processes; compositional dilution (improves solvent IURF by decreasing ω i ) and mixing-induced density increases (reduces IURF by increasing ρ). The compositional dilution effect occurs whenever the solvent MW is larger than the hydrocarbon MW, which causes ω hc to decrease faster than z s increases as www.nature.com/scientificreports www.nature.com/scientificreports/ solvent is added, which is readily apparent for C 1 hydrocarbon (Fig. 6a). However, this effect is limited to CN < 3 because the MW of n-alkanes (14CN + 2) quickly exceeds the solvents (Fig. 6b).
If there were no volume change on mixing, compositional dilution would entirely explain differences in IURF among solvents. However, real fluids experience some volume change on mixing due to changes in the density and number of phases. Thus, solvent IURF performance for any hydrocarbon composition is not strictly ordered based on the solvent MW at a given solvent mole fraction and pressure/temperature. Instead, the controlling factor for IURF is the equilibrium phase density (and the number of phases) at the puff conditions. Using C 1 hydrocarbon as an example, Fig. 6b shows IURF with C 2 > CO 2 > N 2 instead of CO 2 > C 2 > N 2 as predicted by Fig. 6a. C 2 and CO 2 solvent behave similarly across the n-alkane range. However, N 2 becomes dominant at lower pressures for heavy n-alkanes due to strong vaporizing effects (closely spaced horizontal solvent IURF contours in Fig. 5).  Mass fraction of (a) C 1 and (b) C 2 in a hydrocarbon-solvent mixture versus solvent mole fraction for different solvents, C 1 , N 2 , C 2 , and CO 2 . Compositional dilution (gray-filled region) occurs rapidly for C 1 with all solvents because the mass fraction of C 1 changes faster than solvent mole fraction for dilute mixtures. In contrast, rapid compositional dilution only occurs for dilute C 2 -solvent mixtures using CO 2 due to MW effects. (2019) 9:10706 | https://doi.org/10.1038/s41598-019-47099-3 www.nature.com/scientificreports www.nature.com/scientificreports/ Eagle Ford case study. We selected the Eagle Ford (South Texas, USA) for two reasons. First, the Eagle Ford spans the entire PVT window from dry-gas, wet-gas, retrograde condensate, volatile oil, to black oil with CH 4 mole fractions ranging from 72-38% 39 . Second, the Eagle Ford is one of the most prolific/actively developed ultra-tight oil reservoirs in North America with the extensive literature.
In this section, we first analyze primary ultimate recovery factor (URF P ) values for Eagle Ford compositions ranging from dry gas, lean condensate, rich condensates, to oils reported by Orangi et al. 40  Results for pure n-alkanes (Figs 3 and 4) demonstrate that URF P is sensitive to reservoir temperature and initial pressure. Therefore, we standardized reservoir conditions (p I = 8000 psia, T res = 285°F) to reflect the oil and condensate averages and eliminate all factors unrelated to fluid composition and producing pressure when comparing URF P values (Fig. 7). We use equivalent hydrocarbon molecular weight, = ∑ ∑ MW zMW z / eq i i i [lbm/ lbmol], to express fluids as single points in hydrocarbon composition space (the summation excludes non-hydrocarbon components). Figure 7 shows that URF P values are inversely related to MW eq . Recovery improves for all compositions with increasing drawdown (p − p wf ), but heavier compositions are more sensitive due to lower fluid compressibilities. Practical (time-dependent) recovery factor (RF) values will be lower due to preferential production of lighter components, limited production timeframes, and un-accessed reservoir volume.
Next, we evaluate and compare the HnP performance of four different solvent compositions; CH 4 , C 2 H 6 , CO 2 , and N 2 , to improve the URF beyond primary values. We express improvement based on incremental ultimate recovery factor, IURF [OOIP mass %], summed over all components ( = ∑ IURF IURF i n i c ). In Fig. 8, we compare primary solvent IURF [OOIP Mass %] for a single cycle of solvent injection in the early-time limit (ETL) defined by Eq. (3) and the late-time limit (LTL) defined by Eq. (5). In practice, the field-scale HnP process is designed based on injecting a fixed volume of solvent or injecting solvent at a fixed maximum bottom hole pressure 26 . The results show that IURF is proportional to the total amount of solvent injected. Positive IURF values indicate an overall reduction in density (successful compositional dilution) while negative IURF values indicate that density is higher in the puff than it was during primary. Incremental recovery may be positive for individual components (IURF i > 0) even if the overall IURF is negative. C 1 was beneficial to recovery for all Eagle Ford compositions and the heavier MW solvents (CO 2 and C 2 ) were generally detrimental for the condensate compositions. N 2 benefited almost all compositions. However, N 2 was slightly inferior to C 1 for all compositions except the lowest volatility (gas oil ratio (GOR) = 500 SCF/STB). In general, the oil compositions benefited from using any solvent. Larger slugs are preferred whenever additional solvent tends to improve recovery.
In Fig. 9, we compare solvent IURF performance in the early-time and late-time limit based on the same pressure reached at the end of the huff, but assuming different levels of primary reservoir pressure depletion before starting injection (p I = 8000 psia, p I − p wf = 3000, 6000, 7000, and 7500 psi). The total amount of solvent injected progressively increases with drawdown because there is more compressive storage within the early-time volume. In addition, the introduction of solvent tends to increase the fluid compressibility, which allows further amounts of solvent to be injected. Total injected solvent is highest for the gas-like compositions and lowest for the heavy oil compositions. Figure 10 shows the relative RF improvement (IURF/URF P − 1) [%] for the cases in Fig. 9. According to Fig. 10, the ratio of IURF (late-time limit)/IURF (early-time limit) changed by a factor of 2 at p I − p wf = 6000 psi, to a factor of 4 at p I − p wf = 7000 psi, to a factor of 10 at p I − p wf = 7500 psi for oils. The small drawdown scenario (p I − p wf = 3000 psi) showed negligible differences between ETL and LTL IURF. This is because the changes in density and pore volume in Eqs (3) and (5) roughly cancel each other out due to the small amount of solvent injected. In contrast, dry gas and condensate were sensitive to solvent composition, with only C 1 and N 2 solvents providing positive IURF values over the range in p wf and soak duration limits. CO 2 and C 2 solvents performed poorly for condensates, except perhaps very rich condensates. This makes sense because these solvents are similar www.nature.com/scientificreports www.nature.com/scientificreports/ to condensates. The large relative improvement in recovery using N 2 seen in Fig. 10 for the heavy oil (MW eq = 110 for p I − p wf = 6000 psi) and volatile oil (MW eq = 60 for p I − p wf = 3000 psi) indicate that the solvent-hydrocarbon mixture crossed a phase envelope boundary (or changed from one side of a critical point to another).
In summary, solvent IURF is beneficial whenever hydrocarbon density is lower during the puff than primary production. The magnitude of this benefit depends on the initial fluid composition, solvent type, amount injected, pressure, temperature. The ETL and LTL equilibrium fluid compositions have very different solvent strengths (mole fraction of solvent), and thus, very different amounts of density change. Recovery in the LTL tends to be substantially larger than the ETL. The reason is that a small density change in the entire reservoir in the LTL outperforms a very large density change in a small fraction of reservoir volume in the ETL. However, Fig. 9 suggests that HnP is viable (IURF > 0) in the ETL for oils with any solvent (miscible or immiscible), provided that the fractures contained sufficient pore volume to make the effort worthwhile. Perfect mixing minimizes solvent recycling while V wf -only mixing maximizes it. Solvent recycling is a key concern for gas injection projects in fractured reservoirs with low permeability matrix for both technical and economic reasons 17,41 . Thus, soak duration would largely be driven by time-value of production considerations and the costs associated with obtaining solvent, recovering solvent from the produced hydrocarbon, and (re)-injection.

Discussion
In the case of primary production, Cronin et al. 6 illustrated why the ultimate recovery factor from primary production (URF P ) is independent of hydraulic fracture spacing, matrix-fracture contact area, and matrix transport coefficients. These parameters control recovery rates, but they do not control the equilibrium hydrocarbon density that would be reached after a long time (as t → ∞). Vapor-like hydrocarbon fluid compressibilities are approximately 100 times greater than liquid-like hydrocarbon compositions. As a consequence, pressure depletion has a much smaller effect on reducing hydrocarbon density for single-phase oils versus gases (Figs 3, 8). While it may be poor conventional reservoir management, URF improves with the development of a gas-like second phase (L → L + V) due to depletion or compositional changes following solvent injection. For hydrocarbons close to critical points (p c , T c ), small temperature changes control whether the composition behaves more liquid-like versus vapor-like since T c is proportional to MW (Fig. 2). The progressive shift in the URF P contours corresponds to a transition from liquid-like to vapor-like, which happens abruptly in composition (Fig. 3) and temperature space (Fig. 4). www.nature.com/scientificreports www.nature.com/scientificreports/ The compositional proxy (n-alkane series) qualitatively predicted the observed trends for primary production in the Eagle Ford with respect to increasing drawdown, increasing temperature, and increasing equivalent MW. The compositional proxy also revealed the solvent IURF potential for heavier oils. Specifically, the compositional proxy showed solvent assisted IURF was large (Figs 5 and 6) compared to primary recovery (Fig. 3), which means www.nature.com/scientificreports www.nature.com/scientificreports/ that the relative recovery improvement over primary recovery alone (IURF/URF P ) would be substantial for oils, which is precisely the behavior observed in Fig. 10 for oil compositions.
Maximizing long-term hydrocarbon recovery requires net solvent transfer deeper into matrix to compositionally dilute the hydrocarbons in place and assist with volume expansion. However, it has been clearly demonstrated at the field-scale that relatively short injection-soak-puff cycles have improved incremental oil recovery in the South Texas Eagle Ford 42,43 . Single-cycle IURF results from the early-time limit (Figs 8-10) provide evidence that these short-term improvements could be explained based solely on the rapid pressurization and production of fluid in the fractures during multiple cycles of solvent injection and puff. This is conceptually equivalent to www.nature.com/scientificreports www.nature.com/scientificreports/ operating the wellbore-fractures like a set of bellows that are quickly inflated and deflated during multiple cycles of solvent injection and puff. Solvent IURF calculations assumed the most pessimistic assumptions possible; uniform depletion prior to injection and no diffusive transfer between matrix and fractures for the early-time limit (ETL) scenario. In practice, diffusive influx of hydrocarbons would replenish the fractures due to the favorable concentration gradients created due to fluid expansion and compositional dilution caused by the introduction of solvent. Recall that any primary production history would establish density gradients facilitating diffusive recharge of the fractures. Thus, field conditions would be more favorable than the lower limits our calculations reflect.
Optimal cycle frequency, soak duration, and solvent injection pressure is an actively studied problem 1,7,8,11,26,42,[44][45][46][47] . However, we argue that the concept of "optimal" merely reflects the designer's criteria for maximum tolerable solvent recycling and time-value of production, neither of which changes the underlying fact that long term recovery from matrix (where hydrocarbon primarily resides) requires net solvent transfer deeper into matrix. In this way, our state-function (time-independent) solvent IURF methods provide clear guidance on problems that have previously only been addressed from time-dependent perspectives.
Perfect mixing would be impossible if geologic barriers/dislocations exist within a reservoir to prevent transport (even after an infinite time) within a particular gross rock volume. Therefore, inaccessible pore volume should be excluded from the URF calculation. We chose to assume a non-deforming reservoir with a well-connected constant pore volume when presenting the unified theory for URF (primary and solvent) to emphasize the importance of equilibrium phase behavior and simplify the presentation and scope of the study. In addition, we omitted explicit reference to the effects of compaction/consolidation (changes in porosity and/or the bulk reservoir volume), adsorption, gravity, and spatially variable reservoir properties, pressure, and temperatures. Trivial modifications to Equation (1) expand the generality of the URF procedure, where ρ s [lbm/ft 3 ] is the solid (grains) density, ω is [-] is the mass fraction for component i adsorbed in the solid phase s, V b is bulk volume, with a volume integral over the gross reservoir volume honoring the spatial dependencies of the bracketed term at initial I(x, y, z) and final F(x, y, z) conditions with respect to composition, pressure, temperature, confinement effects on phase equilibria and properties, gravity effects, reservoir deformation, and reservoir heterogeneity. For completeness, multiple solid phase constituents (e.g. inorganic versus organic matrix domains 48 ) would be incorporated by modifying the adsorption term is the volume fraction, ρ k [lbm/ft 3 ] is the density, and ω ik [−] is the mass fraction for component i adsorbed within the k th solid-phase constituent and n m is the total number of solid-phase constituents.

Conclusions
We proposed an approach to unify our understanding of the ultimate recovery from primary production and solvent-based processes across the full range of hydrocarbon compositions. The unifying theory is based on equilibrium density changes between the initial and final state and practically adjustable factors such as drawdown and solvent composition/amount injected. We have drawn the following conclusions: 1. Density reduction (fluid compressibility) explains why shale oil reservoirs have extremely low recovery factor (<10%) while gas reservoirs may have recovery factors >60%, with condensates in between. 2. Ultimate hydrocarbon recovery is improved by any process that reduces mass concentration via fluid expansion or compositional dilution. 3. Equilibrium phase behavior provides a transport independent method to calculate the absolute theoretical upper-bound value for mass recovery as t → ∞. 4. Maximizing long-term hydrocarbon recovery requires net solvent transfer deeper into the matrix to dilute the hydrocarbons in place and assist with volume expansion (which eventually requires long soaks). Multiple huff 'n'puff cycles at field pressure will allow for greater net solvent transfer into the matrix. 5. Short term incremental recovery benefits can be attained with minimal soaking, particularly in oil compositions with enough fluid in the early-time volume (wellbore and fractures). 6. The late-time limit (LTL) and early-time limit (ETL) define the theoretical bounds for IURF as a function of how much solvent-hydrocarbon mixing occurs in the reservoir. The IURF performance of a real HnP cycle will reside between these two limits. 7. C 1 and N 2 were the most beneficial solvents over the entire range of fluid compositions (gas, condensates, and oils) based on equal injection pressure and equal injected volume comparisons. CO 2 and C 2 solvent were sensitive to producing pressure and tended to benefit only oil compositions or rich condensates at large drawdown values.
We emphasize that our path independent URF calculations do not replace traditional reservoir simulation. Instead, our calculations provide a practical means to define the maximum recovery attainable from primary production or cyclic solvent processes based on reservoir fluids and the key adjustable factors (drawdown and solvent composition/amount). Indeed, URF calculations explain the astounding recovery contrast from shale gas versus shale oil reservoirs. Although URF may not be achievable during the life-cycle of the field development and production, it provides valuable insights on the role of phase behavior on hydrocarbon recovery from ultratight reservoirs.