Using energy balance to determine pore-scale wettability

a r t i c l e i n f o


Introduction
The wettability of a solid surface controls multiphase flow in porous media which has a wide range of applications including oil recovery, carbon capture and storage, and water flow in the https://doi.org/10.1016/j.jcis.2020.03.074 0021-9797/Ó 2020 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). unsaturated zones [1][2][3][4]. Based on the Young equation, the wettability is defined using contact angle: the angle that the two-fluid interface forms at the solid surface.
Contact angle is traditionally measured from the shape of a droplet placed on a solid surface, which is typically performed using a flat, smooth solid surface composed of a pure mineral [5]. However, the contact angle inside a porous material can differ locally depending on the roughness, mineral composition and degree of wettability alteration of the surface [6][7][8]1]. This wettability alteration is of particular importance for a crude oil-brine-rock system since the rock surface can be altered from its original water-wet to more oil-wet states through the contact with surface active components in the crude oil, resulting in a mixed-wet state where contact angle varies spatially [7,9].
High-resolution three-dimensional X-ray micro-computed tomography (micro-CT) has made it possible to directly observe multiphase flow in porous media. AlRatrout et al. [10] developed an automated method to geometrically measure contact angle from the configuration of the surface and fluid-fluid interface obtained from micro-CT images. The method provides typically a million measurement points along the three-phase contact lines from a billion-voxel image and it has been successfully applied to reservoir rocks to characterize different wettability states [11].
Although the method gives useful information on the spatial distribution of contact angle, there are some limitations. Since it uses a snapshot of the images during displacement, the measured angle may represent a hinging value on pores where displacement has not occurred. Indeed, Akai et al. [12] showed that to reproduce the experimentally observed fluid configurations [11], a higher contact angle than that geometrically measured had to be assigned in direct numerical simulations on the same samples. The other concern is that the accuracy of the measurement is highly dependent on the quality of the images because the measurement is performed at the three-phase contact line where uncertainty in image segmentation can be high.
Blunt et al. [13] have recently proposed a method to measure a thermodynamically consistent contact angle-thermodynamic angle-from micro-CT images based on the energy balance for displacement. The method determines the thermodynamic contact angle from two image datasets in capillary equilibrium. It provides a contact angle at the onset of displacement. Hence, this is the angle to use in the modeling of multiphase flow, specifically, in Mayer-Stowe-Princen calculations of threshold capillary pressure [14][15][16][17]. The method has demonstrated that the thermodynamic contact angle provides a more sensitive characterization of wettability of two experiments performed on sandstones with different wettability states [18,19], compared to the geometrical measurement of contact angle [10].
However, there are two remaining concerns with the use of the thermodynamic contact angle. First, in the original derivation of the method, viscous dissipation was ignored, which means that all the work done was assumed to be converted to surface energy. Second, the method only provides single representative angle for an entire rock sample and it is not clear if the method can give spatial information of wettability states, which is of great importance for the accurate characterization of mixed-wet media [20][21][22]12,13].
The validity of the assumption to ignore viscous dissipation depends on the mode of displacement. For drainage processes, previous experimental studies [23][24][25] have reported that up to 84% of the applied work was dissipated: in this case, the assumption is clearly invalid. For water injection processes, although Morrow [23] reported that the dissipated energy was less than 10% based on imbibition experiments on a bead pack, it is not clear if this high energy efficiency can also be obtained for natural rocks with smaller porosity and pore size. Furthermore, for water injection in a mixed-wet state, there has been no study on viscous dissipation and its impact on energy balance.
In this work we study these two concerns with the thermodynamic method proposed by Blunt et al. [13], using two-phase direct numerical simulation on pore-scale images of complex 3D porous media. We show the validity of the method considering viscous dissipation for water injection in water-wet and mixed-wet porous media which is a common situation during oil recovery from petroleum reservoirs. Furthermore, we demonstrate the use of the method on a pore-by-pore basis to obtain the spatial distribution of wettability.
First, the original method presented in Blunt et al. [13] is extended to account for viscous dissipation in Section 2.1. In Section 3.2, the amount of viscous dissipation is measured through the energy balance of displacement processes with a known prescribed contact angle used in the simulations. Then, the amount of viscous dissipation obtained is used to understand its significance on the calculation of thermodynamic contact angle in mixed-wet states in Section 3.3. Furthermore, in Section 3.4, we demonstrate the applicability of the method to inform the spatial distribution of wettability states by applying the method on a pore-by-pore basis, using simulation results performed for mixed-wet systems where the input contact angle differs poreby-pore.

Thermodynamic contact angle
Energy balance has been used to find the thermodynamic contact angle, h t , by Blunt et al. [13]. In their work, the change in surface energy during a displacement was matched by the work done upon injection of fluids, but the viscous dissipation was ignored. It was assumed that all the work done by the injected fluids was converted to surface energy. In this work, we extend the original equation to account for the impact of viscous dissipation.
Following a thermodynamics approach by Morrow [23], we consider energy balance of two capillary equilibrium states for an isothermal irreversible process in an arbitrary control volume [23,26]: where DW ext is the applied external work to the control volume during displacement between the two states, DE surf is the change in surface energy, E surf , and DE vis is the energy dissipated by viscous forces during displacement between the two states. Since DE vis scales as the square of the flow rate [1], this term can be negligible under the assumption of macroscopically slow displacement. However, even in such a case, locally rapid pore-filling events called Haines jumps [27] can yield a non-negligible positive value of DE vis [25,26], which is lost as heat [1]. Since Eq. (1) is valid for the summation of a series of capillary equilibrium states, and the dissipated energy, DE vis , is always positive, we obtain: The surface energy of the control volume is given by: where A is an area, while the subscripts 1, 2 and s refer to phase 1, phase 2 and the solid, respectively. These subscripts will be used to refer to water, oil, and solid, respectively in Section 3. r is the interfacial tension between these two phases while r 1s and r 2s represent the interfacial tensions between the solid surface, and phases 1 and 2 respectively. The total solid area A s ¼ A 1s þ A 2s is fixed. Hence, Eq.
(3) becomes: The Young equation can be written as [1]: where De is the change in surface energy per unit area and h t is the thermodynamic contact angle. This can be derived from an energy balance at fixed volume, to find the most favorable h t with a constant volume of fluid residing on a solid surface, while the easiest, and standard, way to derive this is from a force balance at the three-phase contact line [23,1]. Inserting Eq. (5) into Eq. (4), the surface energy can be written as: Finally, the change in the surface energy, E surf , is given by: where DA is the change in the surface area. The external work, DW ext , is obtained based on pressurevolume work as: where P and DV are the fluid pressure and the change in the volume, P c is the capillary pressure defined by P c ¼ P 2 À P 1 ; j is the total curvature, / is the porosity, V is the volume of the control volume, and DS 1 is the change in the saturation of phase 1. Here we have used DV 1 ¼ ÀDV 2 ¼ /VDS 1 and the Young-Laplace equation: where r 1 and r 2 are the principal radii of curvature of the interface. We represent the viscous dissipation energy as: where R is the ratio of the viscous dissipation energy to the surface energy. Here, DE vis is the total dissipated energy from all the displacement events between two consecutive capillary equilibrium states. Inserting Eqs. (7), (9), and (11) into Eq. (1), the energy balance is described as: Rearranging Eq. (12), the thermodynamic contact angle, h t , is obtained from: The derivation is strictly correct for infinitesimal changes in equilibrium fluid configurations. When R ¼ 0, this equation reduces to the expression provided in [13]. In experiments, h t and R in Eq. (13) are undermined parameters, while the others, the saturation, the interfacial area, and the interfacial curvature, can be measured from fluid configurations obtained from micro-CT images during displacement. To find h t from experimentally obtained images, we need to assume a value of R, normally that R ¼ 0. This assumption has to be validated. One of the key objectives of this work is that we validate this assumption using numerical simulations with different wettability states where all the terms in Eq. (13) are independently determined.
In simulations, we geometrically impose an input contact angle, h i , which is the angle between the normal vectors of pore wall surface and the fluid interface during displacement -a receding angle for a receding interface and an advancing angle for an advancing interface. When we use a uniform input contact angle, h i , we can determine R as we know h t ¼ h i . Specifically, R is determined from the cross-plot between P DW ext and P DE surf obtained from simulation results, using Eqs. (7) and (8). Assuming R is constant for a series of capillary equilibrium states, from Eqs. (2) and (11), we obtain: Although R is not necessarily constant during displacement, we use this relationship to find an average viscous dissipation ratio, R, from the slope of the cross-plot. This will be discussed in Section 3.2 based on our simulation results.

The color-gradient lattice Boltzmann method
To obtain fluid configurations during two-phase displacement in complex 3D porous media, we used the color-gradient LB method [28] constructed for the D3Q19 lattice velocity model with the Multiple-Relaxation-Time (MRT) collision operator [29,30]. The distribution of the phases was tracked using a color function q N defined by: where q r ðx; tÞ and q b ðx; tÞ are the density of oil and water at position x and time t, respectively. q N takes the value of À1 6 q N 6 1, while the exact location of the oil/water interface is given by The unit normal vector of the oil/water interface, n, is given by: The wettability of the solid surface was modeled using the wetting boundary condition presented in [31] which imposes the correct direction of the phase normal vector, n, at the solid/fluid interface lattice nodes based on the wall normal vector, n s , and input contact angle, h i . This wetting boundary condition can accurately model wettability in complex geometries, and has been validated against analytical solutions [31,32] and experimental results performed on mixed-wet carbonate rocks [33]. A detailed description of our LB model can be found in [31,34] In the following sections, the fluid configurations were obtained based on the value of q N at a lattice node, while the oil/water interface was extracted from the contour surface correspond to q N ¼ 0.

Results
In Section 3.1, the numerical simulations used in this work are described. In Section 3.2, we determine the viscous dissipation ratio, R, for uniform input contact angle cases from the mismatch between the applied external work and surface energy using Eq. (14) with h t ¼ h i . In Section 3.3, the thermodynamic contact angle, h t , for mixed-wet cases where both h t and R are unknown is discussed. R obtained from the uniform h i cases is used to investigate the significance of R on the estimation of h t in mixed-wet media. Lastly, in Section 3.4, we demonstrate the applicability of the thermodynamic contact angle to obtain the spatial distribution of wettability in mixed-wet cases using Eq. (13) on a pore-by-pore basis.

Numerical simulations
Simulations were performed on synthetic images of a bead pack and micro-CT images of a Bentheimer sandstone. These images were composed of a porous domain of 288 Â 288 Â 288 cubic voxels of size 3.58 lm (a size of $1 mm 3 ). The pore structures are shown in Fig. 1. The porosity of these structures were 36% for the bead pack, and 21% for Bentheimer sandstone. The pore spaces of these structures were divided into pore regions, resulting in 630 and 414 pore regions for the bead pack and Bentheimer sandstone, respectively. A pore is defined as a wide region in the pore space; its center is the center of the largest sphere that can fit in the pore space. A pore region contains all the void voxels associated with a pore [1,35]. The mean pore diameter (defined as the diameter of the largest sphere) which accounted for 50% of the pore volume was 77 lm and 63 lm for the bead pack and Bentheimer sandstone, respectively. The distribution of the pore diameters of these structures is shown in Fig. S1 in the Supplementary Material.
Similar to the simulations presented in [36], drainage and water injection were performed with water-and oil-wet porous plates attached to the Àx and þx faces of the porous domain, respectively. These porous plates consisted of a mesh of squares 5 voxels in width and 20 voxels in length with contact angles of 30 and 150 , respectively. These porous plate domains were followed by 10 slices of a complete void space as a buffer domain. This simulation condition mimics a laboratory porous plate capillary pressure measurement. The details of the simulation model are shown in Fig. S2 in the Supplementary Material.
The grid size and time step used in the simulations were Dx ¼ 3:58 lm and Dt ¼ 4:27 Â 10 À7 sec, respectively. Identical densities and viscosities of water and oil were used: The interfacial tension between oil and water was set to r ¼ 25 mN/m. Drainage (DR) was performed with a uniform input contact angle of h i ¼ 45 . Oil was introduced into the porous domain initially filled with water by gradually increasing the capillary pressure. This was performed by increasing the pressure of the oil buffer domain, while maintaining the pressure of the water buffer domain constant, applying constant pressure boundary conditions [37] at x ¼ 0 and x ¼ 348, respectively.
Water injection was performed for three wettability states: uniformly water-wet (WW) cases with h i ¼ 45 , uniformly oil-wet (OW) cases with h i ¼ 135 , and non-uniform mixed-wet (MW) cases. The oil-wet cases will be used to evaluate the significance of viscous dissipation in the mixed-wet cases in Section 3.3, because these oil-wet cases should give the maximum likely viscous dissipation. For the mixed-wet cases, different contact angle values ranging from 45°to 165°were assigned for each pore region based on the oil saturation in a pore region after the drainage simulations -a higher contact angle was assigned for a pore with higher oil saturation (fraction of the void voxels filled with oil in each pore region) after drainage. The input contact angle of the mixed-wet cases had a pore volume weighted average of 90°with a standard deviation of 30°. The contact angles of the mixed-wet cases are shown in Fig. 2. For the water-wet cases, water injection was started immediately after the end of drainage. For the oil-wet and mixed-wet cases, an aging period was simulated before water injection. This aging was performed by continuing the simulations for 20,000 time steps while maintaining the pressure boundary condition at the end of the drainage after altering the input contact angle. This allowed a new position of capillary equilibrium to be attained before water injection.
The capillary number was defined by Ca ¼ lq t =r, where q t is the total flow rate of oil and water. In both drainage and water injection, at each capillary pressure, simulations were performed until Ca < 1 Â 10 À6 to ensure capillary equilibrium. The increment of capillary pressure was chosen to maintain slow displacement, resulting in an average capillary number during displacement smaller than 10 À5 for all the simulations. The capillary number and imposed capillary pressure as a boundary condition are shown in Figs. S3 and S4 in the Supplementary Material. The simulated capillary pressures are shown in Fig. 3, while the simulated capillary pressure for drainage and water injection on the water-wet Bentheimer sandstone is compared with the experimentally measured capillary pressure on the same rock by Raeesi et al. [38] in Fig. S5 in the Supplementary Material. The simulation output was obtained with 14 saturation points for drainage and 21 saturation points for water injection.

Energy efficiency and viscous dissipation
We evaluated the amount of viscous dissipation for the entire porous domain during drainage and water injection using the uniform input contact angle cases (DR, WW, and OW cases) in which the change in surface energy can be estimated based on Eq. (7) with h t ¼ h i -for the mixed-wet cases a representative input contact angle for the system was not known a priori as the input contact angle varied pore-by-pore.
The change in the surface energy, DE surf , was obtained from the simulation results as shown in Fig. 4. The fluid configurations were used to obtain DV w in Eq. (7), while the oil/water and water/solid interfaces were used to obtain DA ow and DA ws in Eq. (7). The oil/water interface was obtained by extracting a contour surface corresponding to q N ¼ 0 from the simulated color function, while the water/solid interface was obtained by applying 600 iterations  of Laplacian smoothing to the boundary surface of water and solid nodes which originally had a staircase shape. This can provide sufficient smoothness while maintaining the original shape of a surface [32]. Using Eq. (7) with h t ¼ h i for sequential simulation outputs, the summation of the change in the surface energy, RDE surf , was obtained.
To obtain the applied external work, DW ext , the simulated fluid pressure and fluid configurations were used. The simulated capillary pressure, P sim c , was obtained with P sim  Fig. 3, because P sim c includes disconnected phases, while the imposed capillary pressure is only responsible for connected phases. The summation of the externally applied work, RDW ext , was obtained using Eq. (8). When micro-CT images during flooding experiments are used, in which the fluid pressure is typically not available, this capillary pressure needs to be estimated from the measured curvature of the oil/water interface using Eq. (10). A description of the measurement of capillary pressure from the interface curvature and an assessment of its accuracy can be found in [32,36].  5 shows the cross-plot between RDW ext and RDE surf during drainage and water injection. We assumed that any mismatch between RDW ext and RDE surf was consumed as viscous dissipation.
These computed values fell in the region expressed by Eq. (2) as shown in transparent blue in the figure. Furthermore, a linear trend between RDW ext and RDE surf suggests that although each local displacement event might result in different amounts of dissipation depending on the size and speed of the event, the total dissipated energy between two capillary equilibrium states was approximately proportional to the change in surface energy. From the figure, we can see that drainage is a process where positive applied work is converted to stored surface energy in the system, while water injection in the water-wet cases is a process where stored surface energy after the drainage is converted to work. Finally, water injection in the oil-wet cases is the same process as drainage in a water-wet state.
Based on Eq. (14), the ratio of viscous dissipation to surface energy, R, was determined from the slope of the cross-plot between RDW ext and RDE surf (Fig. 5) as shown in Table 1. When the slope approaches unity, R approaches 0, meaning a high energy efficiency. In both porous media, for drainage and water injection in the oil-wet state, the amount of viscous dissipation was more than 30% of the surface energy, indicating a relatively low energy efficiency. Bentheimer sandstone with smaller porosity and a smaller pore size showed more viscous dissipation in these cases. For water injection in the water-wet cases, we obtained negative values of R, because the sign of the change in surface energy, DE surf , is negative in these cases, while the sign of the energy dissipated, DE vis , is always positive. The viscous dissipation for the bead pack and Bentheimer was 10% and 3% of the surface energy, respectively, indicating a high energy efficiency. The larger viscous dissipation in the bead pack resulted from faster displacement; the average capillary number during water injection in the water-wet state was Ca ¼ 6 Â 10 À6 for the bead pack, whereas Ca ¼ 2 Â 10 À6 for Bentheimer sandstone.
To compare our simulations with the experimental results obtained with random packings of equal spheres in [23], we also measured the energy efficiency of the processes based on the same definition as used in [23]. Our simulation results were consistent with the experiments of Morrow [23] -the experimentally  measured energy efficiencies of the bead pack in a water-wet state were 79% for drainage and 92.5% for water injection [23], while those obtained from our simulations of drainage and water injection in the water-wet bead pack were 75% and 90%, respectively, which is presented Table S1 in the Supplementary Material.

A representative thermodynamic contact angle
We investigated a representative thermodynamic contact angle, h t , for the mixed-wet cases where R in Eq. (13) could not be directly determined from the simulations. This was conducted through the comparison between the angles obtained without considering the viscous dissipation (R ¼ 0) and those obtained with R for the uniformly oil-wet cases shown in Table 1 which should give a maximum bound for the mixed-wet cases. Fig. 6 shows the computed h t as a function of normalized water saturation, S Ã w , which is defined by:   where S w is the water saturation in the porous domain, S wir is the irreducible water saturation at the beginning of water injection, and S or is the residual oil saturation at the end of water injection. In the figure, the computed angles for 10% 6 S Ã w 6 90% are shown, because the values at the beginning and the end of water injection could not be determined accurately due to the influence from the water-wet and oil-wet porous plates. The computed h t for the uniformly water-wet and oil-wet case are also shown in the figure. As expected, the computed h t showed a good agreement with the input contact angles as indicated by the dashed lines in the figure, because we determined the dissipation ratio, R, so as to satisfy the energy balance with the input contact angle, h i , in Section 3.2.
For these cases, we can make two observations. (1) For water injection in water-wet states which is represented by spontaneous imbibition with an increase of water saturation with positive capillary pressure, h t can be determined even ignoring viscous dissipation owing to the high energy efficiency of the displacementassuming R ¼ 0 resulted in a difference of less than 5°compared to values obtained accounting for viscous dissipation. (2) For water injection in oil-wet states which is represented by forced displacement with an increase of water saturation with negative capillary pressure, viscous dissipation needs to be taken into account to determine h t -ignoring this with R ¼ 0 resulted in undetermined h t with j cos h t j > 1.
For the mixed-wet cases where each pore region had different input contact angles ranging from 45°to 165°, we observed only a small impact of the assumed value of R on the resultant h t with less than 10°difference between the angles obtained with R ¼ 0 and these obtained with R of the oil-wet case. This can be understood from Eq. (13). In the mixed-wet cases, the displacement occurred at P c $ 0, hence j $ 0 (see Fig. 3). In this case, the contribution of the viscous dissipation ratio, R, becomes less significant compared to the term of DA ow in Eq. (13).
In our simulations of the mixed-wet cases, displacement occurred at P c $ À1000 Pa which corresponds to the total curvature j ¼ À40 mm À1 with our input interfacial tension of r ¼ 25 mN/m (see Eq. (10)). This means that when micro-CT images during displacement in mixed-wet porous media are used to find the thermodynamic contact angle, if the capillary pressure normalized by the interfacial tension obtained from the measurement of the total curvature is greater than À40 mm À1 (a negative value whose absolute value is smaller than 40 mm À1 ), Eq. (13) with R ¼ 0 estimates a representative h t within 10°. Furthermore, although the measurement of the interfacial curvature in experiments gives additional uncertainty in Eq. (13), with an image resolution D ¼ 3:5 lm/voxel, this magnitude of curvature (jD ¼ À0:14 voxel À1 ) can be measured accurately-the measurement of the interfacial curvature can be performed within AE11% error for a dimensionless total curvature jjDj < 0:2 voxel À1 [36,12].
For both mixed-wet cases, h t increased with water saturation. This trend has been also observed experimentally where h t was calculated based on X-ray images of water flooding on mixedwet rocks [39,40]. This trend can be explained from the sequence of displacement. For mixed-wet media, water first fills water-wet pores, followed by more oil-wet elements with a decrease of capillary pressure. Since the thermodynamic contact angle is obtained from the difference of two capillary equilibrium states, the angle obtained captures the displacement between these states. Fig. 7 shows the sequence of oil displacement for the mixed-wet cases during water injection. Here, all the pore regions were classified into five groups based on their input contact angle, which ranged from 45°to 165°. For each group, the number of oil filled voxels were counted and normalized by its value at the beginning of water injection. As expected, pores with a lower contact angle filled first at a low water saturation, followed by invasion of more oil-wet pores at a high water saturation. Hence, we first observed a smaller contact angle corresponding to filling in water-wet pores, then a higher contact angle for invasion of more oil-wet pores at a high water saturation. Furthermore, there was less displacement even by the end of water flooding for the more oil-wet pores where more negative entry capillary pressure is required, since displacement in these regions is less favored.

Spatial distribution of contact angle
Since Eq. (13) is valid for an arbitrary control volume, we applied the equation to each pore region by taking each pore region as a control volume, to compute h t on a pore-by-pore basis. We used two sequential simulation results to compute h t for pores where displacement occurred between the two states. This was performed from the first to the last simulation output to cover as many pore regions as possible. Then, h t for each pore region was determined by an average value of the computed angles from different sets of two sequential results. The viscous dissipation was not considered in this calculation-we used R ¼ 0 in Eq. (13). In cases where there were rapid pore filling events in a pore between two consecutive capillary equilibrium states, h t resulted in an undetermined value of j cos h t j > 1, because the dissipated energy by these events was large compared to the surface energy in the pore -these values were discarded in the analysis. This effectively Fig. 7. The sequence of oil displacement for the mixed-wet cases computed for (a) the bead pack and (b) Bentheimer sandstone. All the pore regions were classified into five groups based on their input contact angle, which ranged from 45°to 165°. The number of oil filled voxels on the y-axis was normalized by its value at the beginning of water injection. eliminated rapid pore filling events, while h t for these pores were determined from the displacement at other times. Fig. 8 shows the computed h t for each pore region. Here, similar to the treatment in Section 3.3, we discarded the pores which were in direct contact with the porous plates. The input contact angle, h i , and AE15 difference from h i are shown by the solid and dashed lines, respectively. A quantitative analysis on the difference between h t and h i is shown in Table 2. In both the bead pack and Benthemier sandstone, h t was determined for pores accounting for more than 68% of the pore volume, of which more than 66% of the values were within AE15 of h i . Among the six water injection cases, the mixed-wet state of the bead pack showed the lowest accuracy. The main source of the error came from the overestimate of h t for water-wet pore regions where h i ¼ 45 (see Fig. 8b). Fig. 9 shows an example of pore regions where a significant discrepancy between h t and h i was observed. In this pore region, the fluid meniscus spanned multiple neighboring pores whose h i was greater than 45°. As a result, the threshold capillary pressure for filling was affected by the higher contact angles in these neighboring pore regions during a cooperative pore filling event [41]. Therefore, the thermodynamic angle of this pore was over-estimated. In contrast, for the Bentheimer sandstone which had a lower porosity than that of the bead pack, the oil/water interface existed mostly within single pore regions. Hence, this case resulted in a better accuracy with 83% of the estimated values of h t laying within 15°o f the input values. An example of pore regions of Bentheimer sandstone where the thermodynamic contact angle was accurately determined is shown in Fig. S6 in the Supplementary Material.
For water injection in water-wet and mixed-wet states, which are our main focus in this study, the thermodynamic approach provided a good estimate of contact angle on a pore-by-pore basis, since there was less viscous dissipation in the water-wet cases, while viscous dissipation has only a minor impact on the energy balance for mixed-wet media as discussed in Section 3.3.
Mascini et al. [42] have recently proposed a force based approach to obtain contact angle. In their approach, the contact angle for drainage was determined for each displacement event based on the threshold capillary pressure obtained from images, assuming a cylindrical pore with a circular cross-section. They showed that the contact angles determined with their force based approach gave a narrower distribution than those determined from the thermodynamic approach discussed in this work. In fact, the local displacement events they chose for the analysis are not suitable for the calculation of the thermodynamic contact angle, because these events cause significant viscous dissipation. However, we can see a reasonable agreement between h i and h t in water injection in an oil-wet state where displacement occurs in a drainage mode, Fig. 8c and f, when we discard the undetermined values due to significant viscous dissipation. Hence, the thermodynamic approach can give a reasonable estimate for contact angle on a pore-by-pore basis even for water injection in oil-wet media.
The thermodynamic and force based approaches can be viewed as complementary. The former can determine contact angle regardless of displacement mechanism, but cannot be used for rapid pore filling events as seen in drainage. On the other hand, the latter can determine contact angle when a displacement mechanism can be described by a simple force balance, as in drainage.
Overall, we have shown that the application of Eq. (13) with R ¼ 0 to find the thermodynamic contact angle can be used to determine the wettability on a pore-by-pore basis. It provides an accurate estimate of the spatial distribution of the contact angle when the size of each oil/water meniscus involved in a displacement was confined to single pore regions, while the estimated angle was biased by higher contact angles of neighboring pores when the meniscus spanned multiple pore regions. Choosing a control volume for the thermodynamic calculation that accommodates the meniscus involved in each displacement event might improve the determination of the spatial distribution of the wettability.

Conclusions
The original equation to calculate the thermodynamic contact angle, h t , from energy balance proposed by Blunt et al. [13] has been extended to account for the impact of viscous dissipation. Based on this, the amount of viscous dissipation has been quantified using two-phase direct numerical simulations performed on a bead pack and Bentheimer sandstone with prescribed uniform contact angles of 45°for a uniformly water-wet and 135°for a uniformly oil-wet case. Assuming a viscous dissipation between the extremes of these water-wet and oil-wet cases, the significance of viscous dissipation in mixed-wet media has also been studied.
For water injection in the water-wet media we studied, h t has been accurately determined even ignoring the viscous dissipation, because of the high energy efficiency of the process. On the other hand, for water injection in the oil-wet media, viscous dissipation had to be taken into account because of non-negligible viscous dissipation due to rapid pore-filling events such as Haines jumps.
For water injection in the mixed-wet media studied, we have shown that the viscous dissipation does not have the significant impact on the calculation of h t when the displacement occurs at P c $ 0. Assuming a viscous dissipation of the uniformly oil-wet case as the extreme, the values of h t obtained resulted in a change of less than 10°compared to the values obtained ignoring viscous dissipation.
An increasing trend in the computed h t during the displacement was observed in mixed-wet media, which is consistent with experimental observations [39,40]. This is because displacement first occurs in water-wet pores followed by invasion of more oil-wet pore regions. We also saw that there was less displacement overall in oil-wet pores, where the entry capillary pressure is more negative (less favored).
Furthermore, through the comparison between h t and h i for each pore region, we have shown that the method can be used on a pore-by-pore basis to characterize the spatial distribution of the wettability. The difference between these two angles becomes small when the meniscus involved in the displacement is confined to the control volume used for the calculation of h t .
The spatial distribution of wettability characterized by the thermodynamic contact angle on a pore-by-pore basis could provide key input data for predictive modeling of two-phase flow in mixed-wet media [20][21][22]12].
In our future work, we will test this pore-by-pore based thermodynamic contact angle measurement on the images of unsteady-state experiments imaged with time-resolved highresolution synchrotron X-ray (for instance [40]).

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.