Calculation of Joule – Thomson and isentropic expansion coef ﬁ cients for two-phase mixtures

. Joule – Thomson (JT) and isentropic expansion coef ﬁ cients describe the temperature change induced by a pressure variation under isenthalpic and isentropic conditions, respectively. They are commonly used to model a variety of processes in which either ﬂ uid compression or expansion is involved. While a lot of work has been devoted to inferring the JT coef ﬁ cient from an equation of state when the ﬂ uid is a single phase, little attention has been paid to multiphase ﬂ uids, where phase equilibrium has to be taken into account; previous work has only addressed the construction on the JT inversion curve. In the present paper, we describe and implement an approach to calculate these two coef ﬁ cients for multi-component ﬂ uid systems, including when they form two different phases, liquid, and vapor, in thermodynamic equilibrium. The only ingredients are an equation of state and expressions for the ideal part of the speci ﬁ c heats of the ﬂ uid components. We make use of cubic equations of state, but any thermodynamic model can be used in the proposed framework. Calculations conducted with typical geo ﬂ uids, some of them containing CO 2 , show that these coef ﬁ cients are discontinuous at phase boundaries (where enthalpy and entropy variations exhibit angular points), as expected with any thermodynamic quantity built from ﬁ rst-order derivatives of state functions, and cannot be simply inferred from the coef ﬁ cients of the liquid and vapor phases.


Introduction
Depending on whether the enthalpy H or entropy S of a given fluid system is kept constant, the variation of temperature induced by a pressure variation is described by, respectively, the Joule-Thomson (JT) coefficient and the coefficient of isentropic expansion These coefficients are the slopes of the isenthalps (H = constant) and isentrops (S = constant) in the p-T (pressuretemperature) plane. A positive (resp. negative) value means the fluid experiences cooling (resp. heating) when submitted to depletion. As a rule, the coefficient of isentropic expansion is positive, except in rare instances, e.g., liquid water at temperatures near 0°C or liquid He 3 near 0 K [1], whereas JT coefficients commonly exhibit positive (usually for gases) or negative (for liquids and supercritical fluids) values, depending on the p and T conditions. The Joule-Thomson Inversion Curve (JTIC) in the p-T plane separates the cooling and the heating regions; it connects points where the JT coefficients are equal to zero (or where they change sign when a discontinuity arises). This curve delineates the T and p conditions where cooling occurs upon depletion (l JT > 0) from conditions where heating occurs (l JT < 0). These coefficients are key for describing a variety of processes in which a variation of fluid pressure, whether a compression or an expansion, generates a change in temperature, which is either exploited or avoided, depending on the application. Refrigeration and cryocooling processes exploit the Joule-Thomson effect, by passing a gas through a valve or a throttling device [2]. The efficiency of compressors relies on the isentropic expansion coefficient l S of the refrigerant used, which should indeed not be too high [3]. Hydrocarbon mixtures and other fluids such as CO 2 rapidly flowing to or from an underground porous reservoir, and passing through the wellbore and pipelines and the numerous chokes and valves, experience rapid pressure variations that generate temperature changes. These changes may be so large that a second fluid phase, or even a solid phase (e.g., hydrates [4] or paraffin [5,6]) appears, leading to flow impairment [7]. Temperature well testing is another area where the knowledge of the JT coefficients helps interpret the data in terms of the composition of the produced fluids-gas, oil, or water [8]. Finally, it is worth mentioning here the recent surge in interest for the JT coefficients of CO 2 -rich fluids in the context of CO 2 permanent storage in underground reservoirs. On the one hand, these coefficients play a role if CO 2 is injected in a depleted reservoir, because Joule-Thomson cooling may cause ice or hydrate formation and therefore severe injectivity problems [9]. On the other hand, the same effect provides a safety factor in case of CO 2 leakage from the reservoir, provided the temperature is low enoughbelow 10°C, the maximum temperature of CO 2 hydrate stabilityas is commonly encountered in offshore conditions [8,10].
These coefficients are difficult to measure experimentally, and most often are estimated from an Equation of State (EoS) or from molecular simulations [11]. Kortekaas et al. [12] calculated the JT coefficients of North Sea gas condensates in the single-phase region using conventional cubic EoS. These calculations were however limited to conditions where the fluid forms a single-phase in the p-T domain investigated.
To date and the best of our knowledge, JT and isentropic expansion coefficients have been calculated for single-phase fluids only, except in the case of one-component fluids where simple closed-form expressions exist [12]. In the case of multi-component fluids, the calculations have been limited to the determination of the JTIC in the p-T plane, see e.g. Nichita and Leibovici [13] and Refs. [5,6].
The purpose of this paper is to describe and implement on a few multi-component fluid examples simple calculation schemes based on an EoS for estimating the JT and isentropic expansion coefficients of fluids, both in the singleand two-phase domains. The multi-component fluids chosen are representative of geofluidsnamely, hydrocarbon-rich fluids of underground reservoirs, some of them containing CO 2 . One important expected feature is the discontinuity of these coefficients at the crossing of twophase boundaries, which stems from the fact that these coefficients are related to first-order derivatives of state functions. This feature has been analyzed for one-component fluids by Sychev [14]. In this work, discontinuities for two-phase, multicomponent systems are calculated.
The paper is structured as follows. The next section is a brief reminder of how these two coefficients are related to each other and to other thermodynamic parameters, which are readily calculated from an EoS. The calculation procedure is outlined for a fluid system containing two phases in thermodynamic equilibrium. Examples are given for a variety of oils and gas condensates, i.e., hydrocarbon mixtures exhibiting bubble or dew points, as well as a more complex fluid system. Conclusions are presented at the end of the paper. Three appendices give the expressions of required thermodynamic functions calculated with an EoS, ideal gas heat capacities, and mixture compositions and component properties used in test examples.

The Joule-Thomson and isentropic expansion coefficients
The Joule-Thomson and isentropic expansion coefficients are related to the partial derivatives with respect to pressure p and temperature T of, respectively, the enthalpy H and entropy S of the fluid system by and H and S also depend on the number of moles n i of the various constituents i ¼ 1; . . . ; nc, which are kept constant and for simplicity omitted in the above expressions. In the case of a multicomponent fluid system, with molar content n ¼ ðn 1 ; n 2 ; . . . ; n nc Þ, the enthalpy and entropy are split into an ideal gas term, H ig and S ig , and a departure term directly related to the EoS, H dep and S dep and The above expressions (Eqs. (3)-(6) hold when substituting H and S and their ideal and departure parts with their molar counterparts h = H/n, s = S/n, h id ¼ H id =n, etc., where n ¼ n 1 þ n 2 þ Á Á Á þ n nc is the total number of moles. For cubic EoS, simple closed-form expressions exist for the departure parts h dep ¼ H dep =n and s dep ¼ S dep =n as a function of compositions or mole fractions, z i = n i /n. These expressions are given in Appendix A, for a general form of two-parameter cubic EoS, containing the Soave [15] and Peng-Robinson [16,17] EoS. Note that any thermodynamic model can be used; only the expressions of enthalpy and entropy departures are required. On the other hand, the ideal parts h id ¼ H id =n and s id ¼ S id =n have partial derivatives with respect to temperature expressed as and @s ig @T where C pi T ð Þ is the ideal gas isobaric heat capacity per mole of component i (see Appendix B for the expressions used in this work). The ideal gas enthalpy and entropy do not depend on pressure.
When two distinct phases (liquid and vapor, hereafter identified with subscripts L and V) are present, the enthalpy and entropy of the two-phase fluid system to consider in equations (3)-(6) are the sum of the enthalpies and entropies of the liquid and vapor phases, i.e., H t = H V + H L and S t = S V + S L or, expressed on a molar basis: and where h is the molar fraction of the vapor phase, and x and y are the vectors of liquid and vapor phase compositions, respectively (subscript t stands for total enthalpy and entropy). We assume in this work that thermodynamic equilibrium is ensured, which is legitimate when there are strong enough interdispersion and small domains of the two phases: these conditions hold for instance for a fluid in a porous medium, at least under incipient phase separation, i.e., near bubble or dew point conditions [1]. The quantities h, x, and y are then readily obtained from an EoS by a flash calculation for any given T and p. The partial derivatives of interest are obtained by numerical derivation of the enthalpy and entropy (whose expressions are given below for the ideal parts and in Appendix A for the departure part) along constant-pressure or constanttemperature paths around T and p.
The isobaric ideal gas heat capacity of one mole of the two-phase fluid is The ideal gas part of the molar enthalpy is and the ideal part of entropy is where H ig 0 and S ig 0 are the ideal gas enthalpy and entropy at the reference conditions, respectively. These ideal parts depend only on temperature (not on pressure), hence the partial derivatives of the total enthalpy H t with respect to pressure and temperature (see Eq. (3)) can be expressed as follows: and @H t @T Similarly, the partial derivatives of the total entropy S t with respect to pressure and temperature (see Eq. (4)) are: and @S t @T Note that the partial derivatives with respect to pressure and temperature in equations (14)- (17) are at constant n (mixture composition). In this work, these derivatives are calculated numerically. An analytical treatment would require a framework similar to that presented in Ref. [13], relating the partial derivatives at constant n to those at constant n k (phase compositions) obtained from the EoS, by differentiation of equations (A10) and (A11). Another route to calculate the JT and isentropic expansion coefficients makes use of the differential expression of the enthalpy H and the Maxwell relation leading to @H @p and therefore, inserting equations (20) into (3) and equations (19) into (4), where V is the volume and a is the isobaric thermal expansivity and the isobaric heat capacity. The two coefficients are therefore related as follows: Equations (21) and (22) are equivalent to equations (3) and (4), and offer an alternative route for estimating these coefficients. These equations hold at two-phase pressure and temperature conditions if the extensive thermodynamic functions involved are the sum of the contributions of the two equilibrium phases. When the fluid system contains two phases (liquid and vapor), the volume is the total volume with Z L and Z V the compressibility factors of the liquid and vapor phases (obtained from an EoS, see Appendix A, and a flash calculation), respectively. The two-phase expansivity and isobaric heat capacity are obtained by carrying out flash calculations at various temperatures near the temperature of interest along a constant pressure path. The single-phase JTIC is calculated analytically, as detailed in Refs. [19,20]. In the two-phase region, the location of the JTIC is calculated by tracking sign changes of the numerators (the denominators are always positive) in equations (3) or (21) on a given isotherm or isobar. Partial derivatives are calculated numerically by finite differences. It is incorrect to calculate the JTIC at feed composition in the two-phase domain.

Results and discussion
The single-phase and two-phase Joule-Thomson and isentropic expansion coefficients of several geofluids reported in the literature are calculated here with the Peng-Robinson EoS [16,17] using the approaches described above. Some of the fluid examples and the pressure and temperature conditions chosen here are close to those considered in Ref. [13] where the focus was on the JTIC. The JT coefficients are therefore expected to be small, whether the fluid considered is an oil or a gas condensate at reservoir temperature (a few K/bar), much smaller than the typical values of gases at low pressure. All input data, that is, feed compositions z i , component critical properties (temperature and pressure) acentric factors x i , molecular weights M Wi ; as well as Binary Interaction Parameters (BIPs) k ij are taken from the literature. Riazi [21] and Kesler and Lee correlations [22] are used for the ideal gas isobaric heat capacity, for the light and heavy components, respectively. In the two-phase domain, the partial derivatives of enthalpy and entropy with respect to pressure and temperature are calculated numerically using a second-order central finite difference method (the pressure and temperature derivation steps used were Dp = 0.0001 bar and DT = 0.00001 K, respectively; the sensitivity on the results to perturbations used in the numerical differentiation was carefully checked).
The phase envelope and the JTIC are drawn in Figure 1. The JTIC exhibits three distinct branches: a branch in the single-phase region and an "S-shaped" branch in the twophase region, which are connected by a segment of the bubble point curve. This allows up to three Joule-Thomson inversions on isotherms below the intersection of the twophase branch with the phase boundary. The minimum inversion temperature located in the two-phase region is around 185 K and the maximum inversion temperature located in the single-phase region is about 1625 K. The variation of the Joule-Thomson coefficient in the P-T plane for temperatures from 225 K to 675 K and for pressures between 100 and 500 bar is depicted in Figure 2 (an arrow indicates the location of the phase boundary). Figure 3 plots the reduced departure enthalpy H dep /RT (note that the ideal part of the enthalpy does not depend on pressure) of this reservoir oil at four isotherms: T = 400 K, T = 540 K, T = 615 K, and T = 650 K, in the pressure interval from 100 to 400 bar. The Bubble Point (BP) pressures calculated with the PR EoS [16] for these isotherms are 305.4 bar, 275.65 bar, 240.28 bar, and 220.23 bar, respectively. The figure shows that the departure enthalpy exhibits angular points at the bubble point, i.e. a discontinuity of its derivative with respect to temperature. The single-and two-phase Joule-Thomson coefficients, calculated on the above isotherms in the pressure interval from 100 to 400 bar are plotted in Figure 4. The figure shows that the Joule-Thomson coefficient is negative at 400 K, it changes sign once at 540 and 650 K and three times at 615 K. The results agree with those from Ref. [12]. Except for the first isotherm where it varies in a monotonous manner, the Joule-Thomson coefficient exhibits discontinuities (corresponding to the angular points of the enthalpy in Fig. 3).   The isentropic expansion coefficient of this reservoir fluid is represented in Figure 5 vs. pressure on four isotherms: T = 400 K, T = 540 K, T = 615 K, and T = 650 K, in the pressure interval from 100 to 400 bar. As expected, the isentropic expansion coefficient is positive at all four isotherms, therefore it presents no inversion. In a similar manner to the Joule-Thomson coefficient, the isentropic expansion coefficient varies in a monotonous manner, except at phase boundaries where it exhibits discontinuities (this corresponds to the angular points of the entropy variation with pressure).
The phase envelope and the JTIC are depicted in Figure 6. This mixture exhibits an open-shape phase envelope with no critical point. A dotted line indicated the extrapolation of the single-phase JTIC in the two-phase domain; it has no physical significance since the stable state of the mixture is a two-phase one and it is incorrect to calculate the JT inversion this way (as done for instance in Ref. [13]). There are three distinct branches of the JTIC: a single-phase branch, a portion of the dew point curve, and a two-phase branch. In a small temperature interval, there are three inversion pressures on a given isotherm. The variation of the Joule-Thomson coefficient in the P-T plane for a temperature range from 390 K to 470 K and a pressure range from 350 to 500 bar is drawn in Figure 7 (the location of the dew point curve is indicated by an arrow).
The reduced enthalpy departures of this mixture along three isotherms (T = 400 K, T = 435 K, and T = 470 K) are plotted in Figure 8 against pressure in the interval 400-600 bar. The dew-point pressures calculated with the PR EoS for the three isotherms are 545 bar, 495 bar, and 439 bar, respectively. The departure enthalpies exhibit angular points at dew-point pressures, corresponding to a discontinuity of their derivatives with respect to pressure. The results match those reported in Ref. [13].
The Joule-Thomson coefficients are plotted against pressure in Figure 9 for the above three isotherms in the pressure interval from 400 to 600 bar, with discontinuities at the phase boundary. At T = 430 K, Figure 9b shows three inversion pressures on this isotherm. The isentropic expansion coefficient is drawn vs. pressure for the same three isotherms in Figure 10, exhibiting discontinuities at dew points; it is always positive with no inversion.

Bakken fluid
The third example, referred to as the Bakken fluid, is described by 8 components (four low-molecular-weight alkanes, two intermediate, and two heavy pseudo-components). The composition, component properties (given in Tab. C3) are taken from Nojabaei et al. [24]. The non-zero binary interaction parameters in the EoS are: k 1À2 ¼ 0:005, The phase envelope and the JTIC calculated with the PR EoS are drawn in Figure 11   The Author(s): Science and Technology for Energy Transition 78, 20 (2023) temperatures. This behavior is a reminiscence from pure fluids, for which the Joule-Thomson inversion takes place on the vapor pressure curve at low temperatures, as described earlier by Nichita and Leibovici [13] (see also Sychev [14]). Figure 12 shows the variation of the Joule-Thomson coefficient in the P-T plane for temperatures from 200 K to 500 K and for pressures between 40 and 250 bar.
The calculations are carried out along the isotherm T = 389.3 K in the pressure interval from 100 to 400 bar, where the bubble point pressure is 197.69 bar. The results for the JT and isentropic expansion coefficients are depicted in Figures 13 and 14, respectively. When the pressure decreases along this isotherm and bubble point conditions are crossed, both coefficients experience a discontinuous increase: in the case of the JT coefficient, the jump is from a negative to a positive value.

Conclusions
We have defined and implemented two different but equivalent approaches to extend the calculation of both Joule-Thomson and isentropic expansion coefficients of    The Author(s): Science and Technology for Energy Transition 78, 20 (2023) fluid mixtures from the single-phase to the two-phase region. If one of these coefficients is known, the other can easily be obtained (by adding or subtracting 1/C p , where C p is the total isobaric heat capacity at two-phase conditions). For both coefficients, a discontinuity occurs at liquid/vapor phase boundaries, which correspond to angular points of the enthalpy and entropy variations, respectively.
The JTIC may exhibit several branches: a branch in the single-phase region (down to the intersection with the phase boundary), a branch in the two-phase region, and a portion of the phase boundary itself. Thus, several inversion temperatures/pressures may occur at given pressure/temperature conditions. A rule seems to emerge from JTIC calculations for the test mixtures in this work and confirmed for many other mixtures (not reported here), that is, a two-phase branch exists for open-shaped (or S-shaped) phase envelopes). For closed-phase envelopes, the JT inversion takes place at the phase boundary at low temperatures, as a reminiscence from the behavior of pure fluids, where the JT inversion takes place on the vapor pressure curve at low enough temperatures.
A cubic equation of state is used here, but any thermodynamic model can be used; beyond a two-phase flash calculation routine, only the expressions of the ideal part and of the residual (or departure) part (which is specific to a given EoS) of the enthalpy and entropy are required. In the two-phase region, the required partial derivatives can be calculated numerically, by a finite-difference scheme with the inputs obtained from flash calculations along a constant temperature or constant pressure path.
This work is a first step of a larger project, aimed to study the JT effect and other second-order derivative properties for CO 2 injection in depleted oil and gas reservoirs, as well as the difference between using conventional (at constant pressure and temperature) and isobar-isenthalpic phase equilibrium calculations in compositional reservoir simulation. From practical reasons, in hydrocarbons-brine-CO 2 systems, it is sufficient to use a robust multiphase flash routine and an accurate evaluation of enthalpies in the  The values of X a and X b and the expression of the function m(x i ) are specific to a given EoS. For the PR EoS X a = 0.45724, X a = 0.0778 and m(x i ) is [15,16] m x ð Þ ¼ 0:37464 þ 1:54226x À 0:26992x 2 ; x < 0:49; ðA6Þ m x ð Þ ¼ 0:379642 þ 1:48503x À 0:164423x 2 þ 0:016667x 3 ; x ! 0:49: For the SRK EoS, X a = 0.42748, X b = 0.08664 [14] and m x ð Þ ¼ 0:48508 þ 1:55171x À 0: The implicit form (in compressibility factor Z) of the EoS is obtained by substituting A = ap/R 2 T 2 , B = bp/RT, and Z = pv/RT into equation (A1) The fugacity coefficients are, for i = 1, nc and k = V, L From the general form of the cubic equations of state, h dep k is expressed as follows: where p 0 is the reference pressure (taken equal to 1 bar).