Threshold capillary pressure of caprocks for CO2 storage: Numerical insight on the dynamic and residual method

.


Introduction
Every year around 40 billion tons of CO 2 are emitted into the atmosphere in the world, against a threshold of 18 billion tons that would allow to prevent further global warming (Pörtner et al., 2022).To tackle CO 2 emissions different strategies are proposed, such as capturing and storing carbon dioxide in geological formations (CCS).The technology of underground CO 2 storage is based on the same principles that nature has been employing to keep oil and gas underground for millions of years.Geological trapping of CO 2 is the primary mechanism that prevents CO 2 migration from the storage reservoir to the surface and relies on the existence of a low-permeability caprock formation acting as a seal.To achieve a higher storage volume, the carbon dioxide is compressed and injected in the reservoir in supercritical state scCO 2 (pressure, p ≥ 7.4 MPa and temperature, T ≥ 31.2 • C) and for this reason a minimum depth of 700 to 800 m is desired assuming hydrostatic conditions.
Geological storage should be able to store CO 2 for thousands of years and maintain 99 % of stored CO 2 for at least a thousand years (Bruant et al., 2002).A deep characterisation of the hydraulic and mechanical properties of the formation and potential leakage routes is then required.Possible risks associated with CO 2 leakage to the surface are: asphyxiation and suffocation (continued exposure to concentrations above 20 ÷ 30 % of CO 2 is associated with suffocation to humans and most air-breathing animals) and acidification, increase in calcium and hardness of aquifers (due to high concentrations of carbon dioxide, leading to a release of adsorbed contaminants, especially heavy metals as arsenic, lead and mercury) (Mortezaei et al., 2021).
The ability of caprock to prevent leakage must be evaluated before designing a geological storage operation in a depleted hydrocarbon reservoir, deep saline aquifer, or coal seam.In a CO 2 geological storage process, the CO 2 injected in the formation buoyantly rises until it reaches a low-permeability barrier acting as a seal.Leakage through the caprock can occur due to mechanical failure, capillary failure, or effects related to coupled hydro-chemo-mechanical processes.Mechanical and capillary failures are caused by the overpressure in the fluid saturating the host rock, while the hydro-chemo-mechanical processes are caused by the acidification of the environment by CO 2 that can trigger geochemical reactions, leading to mineralogical and structural alterations of the caprock.Capillary failure occurs when the capillary pressure p c at the interface between the caprock and the reservoir is greater than the threshold capillary pressure of the caprock p * c , i.e. p * c is the lowest capillary pressure at which a continuous advective flow of nonwetting fluid can take place through the overburden.Capillary failure is considered the most likely mode of caprock failure under moderate overpressure conditions, whereas mechanical failure tends to occur for high overpressures (Gluyas and Swarbrick, 2021).The latter may give rise to significant leakage rates in localised areas of the subsurface while the former leads to more pervasive and slower escapes.
Relevant research has been carried out to understand the mechanisms which govern the flow and the entrance (i.e.threshold capillary pressure) of scCO 2 in porous rocks, for instance concerning the role of fabric and scale heterogeneities (e.g.; Rabinovich andCheng, 2020, Moreno andRabinovich, 2021) and the dependency of the interfacial phenomena between the solid, wetting and non-wetting phases (Zhou et al., 2017).
The p * c has been given several names in literature, including entry pressure, threshold capillary pressure, breakthrough pressure, displacement pressure and critical pressure.The precise meaning of these terms varies according to the authors and the experiments carried out, and a good overview is available in Amann-Hildenbrand et al. (2002) and Wu et al. (2020).In this paper, the term threshold capillary pressure will be used.
Direct determination of the threshold capillary pressure in the laboratory requires a step-by-step increase of the non-wetting fluid pressure at the inlet of a water-saturated sample.As long as the difference between the pressure of the non-wetting fluid and the one of the water is smaller than the threshold capillary pressure, water drainage does not take place and no fluid flow is detected at the sample outlet.Conversely, when the difference between the non-wetting fluid pressure and the water one is higher than the threshold capillary pressure, the nonwetting fluid enters the sample and water drainage occurs (Thomas et al., 1968).However, for samples with permeability in the nanoDarcy range, the delay between the non-wetting fluid displacing the water at the inlet and water production at the outlet can be days or weeks and the flow rate can be so small to be hardly noticeable.Furthermore, the method requires a continuous outflow, which can take days to establish if the specific storage is high (Boulin et al., 2013).The procedure is thus very precise but lengthy, especially when small pressure increments are applied.The duration of the experiments depends on the hydraulic characteristics of the sample: for instance, the test on Opalinus Clay in Amann-Hildenbrand et al. (2015) (intrinsic permeability k = 10 − 21 m 2 ) lasted 150 days while it lasted 10 days in the experiment on a mudrock (k = 10 − 19 m 2 ) described in Ito et al. (2011).Testing times can be reduced by applying larger pressure increments, but this leads to overpredicting the threshold pressure (Boulin et al., 2013;Busch and Müller, 2011;Egermann et al., 2006;Li et al., 2005).
Different methods have been proposed for faster estimates of p * c in caprocks: the continuous injection method (Rudd and Pandey, 1973), the residual capillary pressure method (Amann-Hildenbrand et al., 2002), the dynamic threshold pressure method (Egermann et al., 2006) and the racking method (Meyn, 1999).The dynamic threshold capillary method (hereafter only dynamic method) and the residual capillary pressure method (hereafter only residual method) appear to be the most widely used (Kawaura et al., 2013;2014;Boulin et al., 2013;Egermann et al., 2006;Ito et al., 2011;Amann-Hildenbrand et al., 2013;2015;Rezaeyan et al., 2015;Kim and Santamarina, 2013, Minardi et al., 2021, Stavropoulou and Laloui, 2022, Gao et al., 2014, Wollenweber et al., 2010, Espinoza and Santamarina, 2017).In the dynamic method the sample is flushed firstly with water and then with CO 2 , and the threshold capillary pressure is estimated by accounting for the drop in the discharge occurring when the latter enters the sample.Testing times are significantly shorter than those of the step-by-step method, for instance   2013) was about 3 days against 48 days for the step-by-step method.In the residual method, two pressure vessels are connected to the ends of a water-saturated sample.The vessel at the inlet contains CO 2 which is allowed to flow into the sample, whereas the vessel at the outlet contains water and the hydraulic circuit downstream of the latter vessel is closed.The pressure difference between the two vessels evolves and then tends to a value considered a precautionary estimate of the threshold capillary pressure.The testing times are intermediate between those of the dynamic and the step-by-step method (e.g. 2 days for the mudrock in Ito et al. (2011) and 19 days for the claystone Boulin et al. (2013)).A preliminary and fast (< 1 day) evaluation is also possible with the tangent method (Schlömer andKrooss, 1997, Rashid et al., 2015) if the capillary pressure curve of the material for the fluids system of interest is determined through the results of mercury intrusion porosimetry tests.This method consists in plotting the tangent to the inflection point of the capillary pressure curve and an estimate of the threshold capillary pressure is identified as the value of such tangent line for wetting fluid saturation equal to one.As these tests are not designed to directly determine the minimum overpressure which allows the flow of the non-wetting fluid, their interpretation relies on assumptions, discussed in the following sections, which might have an impact on the experimental determinations.
According to experimental results from the literature, the threshold capillary pressure estimated with the dynamic method is relatively close to the one determined with the step-by-step method.The underestimation obtained with the dynamic method respect to the step-by-step one, expressed in terms of percentual error, is below 13 % (8 % on average) for samples of Tavel carbonate (Boulin et al., 2013), chalk, sandstone, and "caprock" (Egermann et al., 2006), mudrock (Kawaura et al., 2014).Overestimations were obtained for Ketzin mudstone (percentual error 25 %, Boulin et al., 2013) and carbonate (percentual error 3 %, Egermann et al., 2006) samples.The residual capillary pressure appears to underestimate the results of the step-by-step method, although to very different extents.According to literature results (Boulin et al., 2013;Egermann et al., 2006;Kawaura et al., 2014;Amann-Hildenbrand et al., 2015;Ito et al., 2011;Rezaeyan et al., 2015) the underestimation ranges from 8 % for the mudstone studied by Ito et al. (2011) to 78 % on the shale analysed by Rezaeyan et al. (2015).
Although specimens from the same sample were tested in all these studies, it shall be noticed that the discrepancies might depend as well on small sample heterogeneities and on the fact that in the step-by-step method the pressure increase is applied through finite increments, which leads naturally to (limited) overestimations of the actual threshold pressure.
Despite being evident from the literature that the actual value of the threshold capillary pressure is different from the ones estimated with the dynamic and the residual methods, this evidence seems to be built on the comparison of the experimental results rather than on a theoretical/ numerical analysis of the physical reasons for the difference observed.As only exception, numerical simulations of residual tests on Opalinus Clay were presented in Kivi et al. (2022), suggesting that for this material the residual pressure values depend on the sample length and vessel volume.As the boundary conditions applied to the sample are different, different two-phase flow processes occur in the two methods.Two-phase flow is mainly governed by the capillary pressure curves and the relative permeability of the material, which then also control the estimated value of the threshold capillary pressure.Unfortunately, the hydraulic properties of the samples are often missing in the experimental characterization provided alongside the results of the tests.Numerical sensitivity analyses might then aid in gathering an insight into the relevant aspects of the tests, considering also that a parametric analysis accounting for the hydraulic properties of the material and for the physical processes taking place during the tests has, to the best of our knowledge, not been published thus far.
This paper aims to provide insight into the physical processes taking place during tests performed according to the dynamic and the residual method.Numerical models were implemented through the Finite Element code Comsol Multiphysics® to support the interpretation of these tests, to explore the limitations and advantages of the two methodologies and to propose experimental measures in order to obtain the most accurate estimate of the threshold capillary pressure.The numerical sensitivity analyses carried out allow us to underline the effects of the hydraulic properties and the test conditions on the estimate.
In the following, the tangent, the dynamic and the residual methods are briefly recalled.The properties of different virtual (i.e.numerical) samples studied in this work are then introduced.Samples were chosen to be realistic and according to the tangent method they all have the same threshold capillary pressure.The predictions of the results that would be obtained with each method are shown together with the physical processes that occur during the test.A discussion section outlines the differences between the estimates of the threshold capillary pressure as would be determined through each method, accounting for the role of the material properties and testing set-up on the results, and how the best estimate could be obtained with both methods.

Tangent method
The tangent method (Schlömer and Krooss, 1997;Rashid et al., 2015) allows a preliminary evaluation of the threshold capillary pressure and it relies on the analysis of the primary drainage branch of the capillary pressure curve of the fluids system of interest.It consists in plotting the tangent to the inflection point of such a curve, where the capillary pressure p c is drawn in the logarithm scale.An estimate of the threshold capillary pressure, hereafter named p .. c , is identified as the value of such tangent line for a water saturation S w = 1.The capillary pressure curve can be obtained by the centrifuge or porous plate methods (Hassler and Brunner, 1945;Falode and Manuel, 2014) or, more simply and quickly, by mercury intrusion porosimetry (MIP) (Monicard, 1975).In the latter case, the capillary pressure curve obtained refers to the mercury-air system as well as the threshold capillary pressure p ..  c(Hg,air) which is obtained.To achieve the threshold capillary pressure concerning the fluid system of interest (e.g.scCO 2 -water, p .. c(scCO2,H2O) ), it will be necessary to convert p .. c(Hg,air) using the Washburn-Laplace equation: Hg,air) γ (scCO2,H2O) cosθ (scCO2,H2O) γ (Hg,air) cosθ (Hg,air) where γ is the interfacial tension between the wetting and the nonwetting fluid and θ is the contact angle between the wetting fluid and the solid.Some limitations concerning the representativeness of the MIP tests for field conditions make the estimated p .. c only roughly similar to p * c (Egermann et al., 2006;Minardi et al., 2021;Wu et al., 2020) since: it is not possible to impose the stress conditions of the site during the experiments, and only the entrance radius of pores is detected while their connectivity is ignored.Nevertheless, the MIP is a fast and simple technique that can provide a consistent order of magnitude of the p * c which can be used to design properly more accurate tests.In Fig. 1(a), the tangent method is applied to the capillary pressure curves of the virtual samples explored in this work, further discussed in Section 3. It can be noticed that, despite having different capillary pressure curves, the three samples were chosen so to have all the same value p .. c = 6.70 MPa.The adopted capillary pressure curves were chosen to be compatible with those of real caprocks, in particular the capillary pressure curve of sample B is very similar to the one of the Opalinus Clay Shale in Favero & Laloui (2018) (see also Fig. 1(b)), whose threshold capillary pressure was determined by Minardi et al. (2021).

Dynamic threshold pressure method
The dynamic method was firstly proposed by Egermann et al. (2006) and a scheme of the setup is reported in Fig. 2. The method consists of injecting the non-wetting fluid into the water-saturated sample at a constant pressure p inlet nw such that the pressure difference Δp t = p inlet nw − p outlet w imposed between the inlet and the outlet of the specimen is greater than the actual threshold capillary pressure p * c .Initially, the drainage pipes at the inlet of the sample are saturated with water and a single-phase flow of water is imposed.
The intrinsic permeability k is determined through the outgoing discharge Q w : where L is the specimen length, A is the specimen section area, μ w is the dynamic viscosity of water.
As the non-wetting fluid starts penetrating the specimen there is a significant decrease in the outgoing discharge (hereafter called effective discharge) Q eff w , due to the capillary pressure jump at the interface between water and scCO 2 .Egermann et al. (2006) implicitly assume a penetration front of the non-wetting fluid that divides the volume of the sample into a portion saturated with non-wetting fluid and one saturated with water.The entry of the non-wetting fluid into the pores involves a pressure jump at the front p ∘ c = p front nw − p front w , with p ∘ c assumed equal to the threshold capillary pressure of the porous medium p * c .The total pressure drop imposed across the specimen can then be decomposed as follows:  where Δp nw is the pressure drop along the sample volume saturated with the non-wetting fluid and Δp w is the pressure drop along the sample volume saturated with water (Fig. 3).
In determining p ∘ c , the Authors consider the initial instants in which the non-wetting fluid permeates the inlet of the sample (hence Δp nw ≅ 0), and they determine Δp w from the measurement of the effective discharge Q eff w : Combining the Eqs.( 2)-( 4) it follows: The method requires then only the value of Δp t and of the discharges Q w and Q eff w to be known (Fig. 4).It shall be noticed that the implicit assumption of Eq (3), i.e. of a 'piston like' displacement of water because of CO 2 penetration, implies that the sample is fully saturated by the CO 2 when the capillary pressure is larger than the threshold one.It also does not account for viscous fingering phenomena.Although desaturation is surely incremental and viscous fingering might take place as CO 2 enters the sample, the method should be applied at an instant when the penetration of CO 2 is negligible (Δp nw ≅ 0), when this condition has no practical relevance.

Residual capillary pressure method
The residual method was firstly proposed by Amann-Hildenbrand et al. (2002).It is based on the transient flow of the non-wetting fluid into the water-saturated sample.A scheme of the setup for the residual method is reported in Fig. 5. Two vessels, containing the non-wetting fluid and the wetting fluid respectively, are connected at the inlet and the outlet of the sample.The pressure in the two vessels is recorded during the entire experiment.At the beginning of the test, the pressure difference between the two vessels is set above the expected value of the threshold capillary pressure and flow of non-wetting fluid from the high-pressure vessel (at the inlet of the specimen) to the low-pressure one (at the outlet) takes place.As the non-wetting fluid enters the sample, a drainage process starts.The drainage line of the outlet vessel is closed.Such undrained condition causes an increase in the water pressure at the outlet, recorded during the test, which lately also induces an imbibition process that interrupts the flow of the non-wetting fluid through the sample.Two different procedures can be followed when running the test, in the following named procedure I and procedure II.According to procedure I, the non-wetting fluid pressure in the inlet vessel is kept constant throughout the test.According to procedure II, such pressure is kept constant for a short time, then the drainage line of the inlet vessel is closed and the inlet pressure decreases with time, as done in the Pulse Decay test (Brace et al., 1968).The pressure difference between the inlet and the outlet vessels in stationary conditions is called residual capillary pressure p ∧ c and it is used to estimate the threshold capillary pressure (Fig. 6).As both a drainage and an imbibition path occur during the test, Zweigel et al. (2005) suggested that the residual capillary pressure should be lower than the p * c , as found in Egermann et al. (2006), Boulin et al. (2013), Kawaura et al. (2014), Rezaeyan et al., 2015.Therefore it has been considered a prudent estimate of the threshold capillary pressure.Amann-Hildenbrand et al. ( 2012) claim that the residual capillary pressure is the capillary pressure for which, following the imbibition of the wetting fluid, the non-wetting fluid is no longer a continuous phase within the porous medium ('snap-off' pressure).

Numerical sensitivity analyses
The dependency of the results of the dynamic and the residual method on the capillary pressure curve and on the relative permeability of the materials was investigated through 1D numerical models implemented in the FEM code Comsol Multiphysics®.The role of some experimental aspects such as the scCO 2 pressure applied at the inlet of the specimen (both methods), the size of the vessels and the holding time of the scCO 2 pressure at the inlet before closing the hydraulic circuit (residual method) were investigated as well.
In the simulations, the mass balance equations for the wetting and the non-wetting fluid were solved at the macroscopic ('Darcy') scale.The two fluids were considered to be compressible and immiscible, whereas the sample volume and the pore network were assumed to remain constant during the tests, i.e. no volume strains and no opening of new microfractures were considered to take place.Although the intrinsic permeability of some clayey materials might depend on the saturating fluid (see e.g.Gupta et al., 1991;Behnsen and Faulkner, 2011;Romero et al., 2011) the reference intrinsic permeability at full saturation was assumed to be the same for both the wetting and the non-wetting fluid.This assumption was introduced as the dependency of the permeability on the pore fluid composition is highly related to the mineralogy, while this work focuses on more general aspects.
For the sake of clarity, the hydraulic coupled model, the modelled geometry, the mesh and the boundary conditions adopted in the simulations are provided in the Appendix, while the discussion in the present section is limited to the material and experimental conditions that were varied in the analyses.As a way to validate the model and then the results of the sensitivity analysis, the Appendix also provides the results of the simulations of a dynamic test in Boulin et al. (2013) and of the residual test in Minardi et al. (2021).These two tests were chosen as they are among the few in the literature for which information regarding the capillary pressure curves are also provided.2006) assumption (modified by Egermann et al., 2006).

Pore size densities and capillary pressure curves
Three synthetic materials having different pore size distributions and relative permeabilities were selected for the numerical analyses.The pore size density curves (PSDs) are presented in Fig. 7.They range from a rather heterogeneous distribution (material A) to a more uniform one (material C).
PSDs represent the fractional void ratio (ratio between the volume of voids and volume of solids, see e.g.Bear and Cheng, 2010) for a certain pore size and are used to visualize the relevance of the classes of pores.They are related to the relative frequency of pores having radius r, f(r), through the expression: where e = ϕ /(1-ϕ) is the total void ratio of the sample.The capillary pressure p c at which a non-wetting fluid can enter a cylindrical capillary tube of a given radius r depends on the interfacial tension γ between the wetting and the non-wetting fluid and on the contact angle θ between the wetting fluid and the solid (Washburn-Laplace equation): According to Millington and Quirk (1961) and Mualem (1976) the porous medium can be approximated with a bundle of capillary tubes having radii between R min (minimum radius) and R max (maximum radius).For a given capillary pressure the wetting fluid fills the tubes with radius r < R and the non-wetting fluid fills those with radius r ≥ R. The degree of saturation of the wetting fluid is then (Della Vecchia et al., 2015): The main drainage branches of the capillary pressure curves (S w , p c ) presented in Fig. 1(a) correspond to the PSDs of Fig. 7, and were obtained through Eqs. ( 7) and ( 8) .The main wetting curves were obtained by adopting the same PSDs and an advancing contact angle of the scCO 2 -H 2 O system equal to 50 • (Huang et al., 2019;Jafari and Jung, 2017).

Models of capillary pressure curve and relative permeability
The capillary pressure curves were fitted with the Brooks & Corey (1964) and the van Genuchten (1980) relationships.These relationships were chosen among others since they are widely used and because they are coupled with relative permeability formulations with very different trends.Brooks & Corey (1964) capillary pressure formulation reads: where λ is called the pore size distribution index, p b is the bubbling pressure, and Ŝw is the effective water saturation: where S wr is the irreducible water saturation, imposed equal to 0.03.van Genuchten (1980) formulation reads: where α is a parameter related to the entry capillary pressure of the nonwetting fluid and n and m are related to the pore size distribution of the material.
The values of the parameters of the capillary curves of the three materials with PSDs are provided in Table 1: A linear scanning equation was defined for the transition from the main wetting to the main drying curve, which is of interest for the simulation of the residual method: The unsaturated permeability to the wetting fluid k w and the nonwetting fluid k nw were defined by: where k is the intrinsic permeability, k r,w is the relative permeability to the wetting fluid, and k r,nw is the relative permeability to the nonwetting fluid.
The relationship between relative permeability and saturation is likely to be material-dependent.One of the purposes of the paper is to outline how the relative permeability influences the test results.To allow this, the relationships associated to the Brooks & Corey (1964) and the van Genuchten (1980) curves were adopted.Brooks & Corey expressions are as follows: For the wetting phase, van Genuchten (1980) suggested: whereas for the non-wetting fluid Parker et al. (1987) obtained the following expression by coupling Eq. ( 11) with the model in Mualem (1976): Fig. 8(a) shows the relationship between the relative permeability to the wetting fluid and the water saturation as predicted by Eqs. ( 15) and ( 17), while Fig. 8(b) shows the same relationship for the non-wetting fluid, as predicted by Eqs. ( 16) and (18).

Table 1
Parameters of capillary pressure curves (Brooks and Corey, 1966;van Genuchten, 1980).In the range of water saturation of interest for the performed numerical analyses (approximately S w > 0.80), the Brooks and Corey formulation always predicts a higher relative permeability to the wetting fluid and a lower relative permeability to the non-wetting one than the van Genuchten formulation.The difference between the two formulations increases as the heterogeneity in the pore size distribution increases, i.e. it is larger for material A in relation to material C. For both formulations, an increase in pore size heterogeneity leads to a decrease of the relative permeability to the wetting fluid and an increase of the one to the non-wetting fluid.
It shall also be observed that, according to Schowalter (1979), in homogeneous samples the threshold capillary pressure is linked to water saturation in the range 0.85 < S w < 0.95.For such values of S w , the relative permeability to the non-wetting fluid predicted by the van Genucthen relationship for material A is still quite high (0.08 < k vG r,nw < 0.25), suggesting then a relatively high flow of the non-wetting fluid even for such high water saturation.This case might thus be more representative of samples where a contribution to flow is also offered by a system of small fractures parallel to the direction of flow.

Testing conditions
Different pressures of the non-wetting fluid at the inlet of the specimen p inlet nw were considered to study its effect on the results: p inlet nw was set equal to a⋅p .. c + p 0 w and the analysed scenarios for both methods are reported in Table 2.In all analyses the initial pressure of wetting p 0 w and non-wetting pressure p 0 nw were imposed equal to 15 MPa and zero respectively.
For the residual method, the effect of the vessel volume V placed at the two ends of the specimen and of the holding time τ of the nonwetting fluid pressure at the inlet, before closing the hydraulic circuit, was studied considering the scenarios summarized in Tables 3 and 4. The vessel volumes were chosen to be consistent with the ones declared in Amann-Hildenbrand et al. (2002).All the analyses were set up to reproduce 21 days (504 h) of testing.

Dynamic threshold pressure method
To study the reliability of the dynamic method and the associated physical processes, 42 scenarios of experimental tests were simulated, which differ from each other in the properties of the caprock and the non-wetting pressure applied at the inlet (Table 2).The simulations reproduced the testing stages: water flow was imposed for 6 h, and after this time, the scCO 2 injection was simulated.For the sake of clarity, the hydraulic coupled model, the modelled geometry, the mesh and the boundary conditions adopted in the simulations are provided in the Appendix.
In the following, the physical processes taking place during the dynamic tests will be examined and discussed concerning the properties of material B with Brooks and Corey relative permeability.In particular, the influence of the pressure difference imposed across the sample will be examined, in terms of the temporal evolution of the effective discharge Q eff w (Fig. 9), of profiles of the wetting and non-wetting pressures (Fig. 10) and of scCO 2 saturation profiles (Fig. 11).Insights on the

Table 3
Reservoir size used in the sensitivity analyses (residual method).

Table 4
Holding time of the non-wetting fluid pressure at the inlet used in the sensitivity analyses (residual method).
Holding time of p inlet nw 0 6 role of the material properties on the results will be addressed in the Section 5.
The key parameter in the interpretation of the dynamic method (Eq. ( 5)) is the ratio between the effective discharge Q eff w and the initial discharge Q w (due to the single-phase flow of water).In all simulations, such a ratio was found not to be constant, as the discharge that occurs after that scCO 2 reaches the inlet evolves with time.Evaluating the time at which Q eff w Qw should be determined to allow a proper use of the method is then crucial.This aspect will be discussed in the following lines, considering that according to the assumptions of the method the scCO 2 penetration should be limited, so that the pressure drop of the nonwetting phase is Δp nw ≅ 0. Fig. 9 shows the time evolution of the ratio Q eff w Qw provided by the numerical analyses when different values of p inlet nw are imposed.It can be appreciated that the effective discharge decreases with time for all the pressures, and it tends to stabilize at different times depending on the applied pressure.Even when Δp t is very small (i.e.Δp t = 0.67 MPa, p inlet nw = 15.67 MPa), the ratio of the discharges goes to zero only after a transient phase which lasts about 6 h.This phenomenon can be explained by considering that when scCO 2 reaches the inlet water stops entering the sample.However, the water pressure inside the sample still exceeds the one applied at the outlet.Such overpressure dissipates with time and after some time the pressure is everywhere equal to the one at the outlet.Only at this time the discharge also drops to zero.Pore pressure dissipation requires the onset of a transient flow process which generally depends on the hydraulic and mechanical properties of the porous medium and of its constituents.When the capillary pressure at the inlet is well below the threshold pressure, as it is in this case, no scCO 2 enters the sample and under the assumption of no volume deformations the transient process is controlled only by the compressibility and the viscosity of water and by the permeability of the sample.The diffusion equation governs the phenomenon under the hypothesis of constant water compressibility and constant permeability (see e.g.Biot, 1941).
For values of p inlet nw smaller than 23.71 MPa the time evolution of Q eff w Qw is about the same and, in particular, the discharge at large times drops to zero.Small differences are only due to the effect of a small amount of scCO 2 which accumulates at the inlet without actually flowing through the sample (i.e. it only desaturates the sample surface).As confirmed by the numerical analysis, no appreciable penetration of scCO 2 occurs in these cases.
For higher values of p inlet nw the time evolution of Q eff w Qw follows the same trend of the previous cases only initially, whereas it departs from it at a certain time.The instant when the solutions diverge is a function of with respect to the previous cases and the persistence of an effective discharge at large times is due to the penetration of scCO 2 , which replaces a fraction of the water inside the pores.This analysis of the numerical data suggests that, to comply with the requirements of the dynamic method, Q eff w Qw should be evaluated at the time when it diverges from the value obtained applying a Δp t which does not exceeds the threshold capillary pressure, as this time indicates the onset of scCO 2 penetration.
The time at which the model assumptions hold can also be determined considering the time evolutions of the profiles of pressure and scCO 2 saturation within the sample.These are provided in Figs. 10 and 11 for the case p inlet nw = 23.71MPa (Δp t = 1.3 p .. c ).The water pressure drops linearly during the first stage of the test, until scCO 2 reaches the inlet (t inj = 0 min, where t inj is the time elapsed from the instant when scCO 2 reaches the inlet).At this time, the scCO 2 pressure is everywhere equal to zero.During the first 1.83 h from injection, (t inj = 1 min ÷ 110 min) the water pressure decreases all along the sample while a scCO 2 overpressure accumulates only in the inlet zone (see also scCO 2 saturation profiles in Fig. 11).As long as no CO 2 is flowing, the drop in water pressure is governed by the diffusion equation as commented above.The steep slope in the scCO 2 pressure and saturation profiles starts progressing into the sample between t inj = 110 min and t inj =120 min.The capillary pressure of the inlet at t inj =110 min, evaluated as the difference between the scCO 2 and water pressure, is therefore equal to the actual threshold pressure of the sample, which is then p * c = 6.91 MPa.After this time, scCO 2 further penetrates, and 24 h after injection the CO 2 stretches along the first 5 mm of the sample.
By looking at Fig. 9, it can also be observed that the time t inj = 110 min is very close to the one at which the Q eff w Qw ratio of the case at hand diverges from the ones obtained for smaller values of p inlet nw .Introducing the effective discharge of this time into Eq.( 5), an estimate p ∘ c = 6.80 MPa is obtained, very close to p * c = 6.91 MPa determined by the pressure profiles.It can also be observed that p * c is larger than 6.70 MPa, a Δp t for which according to the interpretation provided in Fig. 9 no flow of CO 2 occurs.
The determination of the time at which the scCO 2 begins flowing inside the sample could be based then either on profiles of pressures and saturation, or through the plots of Fig. 9. Experimentally, one might think to repeat the test imposing two different values of Δp t , one well below the threshold capillary pressure only to determine the evolution Qw which occurs without scCO 2 flow, and another above the threshold to determine the same ratio at the time of divergence.Although this determination is quite easy in the numerical analyses it is not in usual laboratory tests, where an accurate measurement of very small discharges is problematic, especially considering that at the time of concern those are far from being constant.
Because of these limitations, the conventional assumption made in the literature is to consider the about stable Q eff w , which as shown in Fig. 9 is obtained at later times (Boulin et al., 2013;Kawaura et al., 2014;Egermann et al., 2006).As such value of Q eff w is smaller than the one for which the scCO 2 has not yet penetrated inside the specimen and it is only at the inlet, the p ∘ c overestimates the actual threshold pressure of the caprock.Overestimation is also related to the scCO 2 saturation profiles at the considered times (Fig. 11), as they are different from the 'piston like' assumption of the method (Fig. 3).The results of the numerical analyses suggest that the overestimation of the threshold pressure, referred by Busch & Müller (2011) and Zhang & Wang (2022), might be due to the use of an effective discharge determined at a time when the method assumptions are not satisfied, rather than to an intrinsic limitation of the method.
The drop in Q eff w occurring between the time at which the method assumptions are satisfied and the time of flow stabilization increases with the inlet pressure.Therefore, also the p ∘ c overestimation increases with the inlet pressure.Estimated values of p ∘ c using Eq. ( 5) with Q eff w Qw at 24 h are presented in Table 5, together with the error of each estimate with respect to the actual value.
While the estimates of the threshold capillary pressure provided in Table 5 depend on the applied Δp t , the analysis of these results still allows stating that p * c is between 6.70 MPa (pressure difference Δp t for which the effective discharge at 24 h is zero) and 7.84 MPa (p ∘ c obtained for the smallest Δp t of Table 5 which ensures that the effective discharge at 24 h is different from zero).

Residual capillary pressure method
To investigate the physical processes that occur during a threshold test with the residual method and understand what are the effects that the properties of caprock and the test conditions (inlet scCO 2 pressure (Table 2), reservoir size (Table 3) and holding time of inlet scCO 2 pressure (Table 4)) have on the residual capillary pressure, 108 different scenarios were simulated.To reproduce the residual method, the upstream and downstream vessels were also simulated.To simulate procedure I, the first stage of rapid increase of the non-wetting pressure at the inlet (30 s) from the initial value to the target one (Table 2), and the second stage with the holding of this pressure were considered.In both stages the condition of closed hydraulic circuit at the outlet was simulated.To reproduce procedure II, the three stages of the procedure were simulated: rapid increase of the non-wetting pressure at the inlet (30 s)  from the initial value to the target (Table 2), holding of this pressure for the target time τ (Table 4), hydraulic circuit closure at the inlet.In all stages the condition of closed hydraulic circuit at the outlet was simulated.For clarity, the Appendix provides information on the hydraulic coupled model, modelled geometry, mesh, and boundary conditions used in the simulations.Consistently with Zweigel et al. (2005) observations, the results of simulations show that the specimen is subject to both drainage and imbibition during the test.The imbibition begins a few hours after the start of the test and it involves only the first 5 mm of the specimen, while the remaining volume experiences only a slight desaturation in all considered scenarios.This can be appreciated by looking at the profiles of scCO 2 saturation at the instant of maximum saturation at the inlet and at the end of the simulation (Fig. 12).By way of example, the sequence of hydraulic states along the capillary pressure curve ('hydraulic paths') taking place at the inlet, the center and the outlet of the specimen during the 21 days of testing according to the two different procedures are shown in Fig. 13.Such results refer to materials A and C, with both the van Genuchten (continuous lines) or the Brooks and Corey (dotted lines) relative permeabilities, tested with an inlet pressure p inlet nw = 33.09MPa, a tank with size V = 1 cm 3 and, for procedure II, a holding time τ = 0 h (closing of the inlet hydraulic circuit immediately after reaching the target scCO 2 pressure).It can be observed that the hydraulic paths are different according to the position within the sample.In all simulations and for all the adopted parameters, all points initially move along the main drainage branch of the capillary pressure curve, whereas imbibition (firstly along a scanning curve and then along the main imbibition curve) takes place at later stages of the test only at the inlet.Moreover, the Brooks and Corey scenarios show a higher water desaturation near the inlet compared to the van Genuchten ones.This evidence is possibly due to the higher relative permeability to the wetting fluid of the Brooks and Corey formulation in the range of water saturations explored in the analyses, i.e. S w > 0.80 (Fig. 8(a)).This effect is larger as the heterogeneity of the pore size distribution increases (material A compared to C), and it is justified by the increasing difference in relative permeabilities between the two formulations with the pore size heterogeneity.For both relative permeability formulations the average (whole sample) and local (inlet) scCO 2 saturation of material C are higher than the corresponding ones of material A (Fig. 13).Furthermore procedure I leads to a higher penetration of the nonwetting fluid into the sample than procedure II with holding time set to zero (τ I ).This depends on the boundary condition imposed on the inlet, since the non-wetting pressure remains constant in time with procedure I, whereas it immediately decays with procedure II and holding time equal to zero (τ I ).
The residual capillary pressure p ∧ c obtained for material A with Brooks and Corey relative permeability is shown in Fig. 14.The dependency of the time evolution of wetting and non-wetting pressures in the vessels on the inlet pressure, the vessel size and the procedure type are also illustrated.
With both procedures the residual capillary pressure p ∧ c increases with the p inlet nw and the size of the reservoir V.It must be considered that p ∧ c increases with the volume of non-wetting fluid which can be stored in the sample, and the undrained conditions at the sample outlet only allow scCO 2 to enter because of the compression of water.Higher pressure and bigger vessel allow a larger volume decrease of the water in the specimen, and thus higher capillary pressures.Procedure I also leads to higher p ∧ c values with respect to procedure II (τ I ) as the pressure at the inlet is constant and equal to a maximum value in the first case while it decreases immediately with time in the second one.Under the first condition, water compression is greater and more scCO 2 can be stored in the sample so, accordingly with the capillary pressure curve, the residual capillary pressure is higher.Similarly, as the holding time τ increases, the determined p ∧ c increases reaching values similar to those of procedure I for holding times of 6 h (τ II ).

Discussion
The threshold capillary pressure p * c identifies the overpressure above which a non-wetting fluid can penetrate and flow into a porous medium initially saturated by a wetting fluid.With the dynamic method, regarding material B with Brooks and Corey relative permeability it was shown in Section 4.1 that p * c can be determined if the effective discharge Q eff w is evaluated at a very short time from when scCO 2 reaches the inlet.This occurs for all scenarios and therefore confirms that the dynamic method is very quick compared to the step-by-step one.The threshold capillary pressure were determined for all the materials and relative permeability relationships, according to the pressure profiles of the two fluids as done in Section 4.1, and they are summarized in Table 6.
The results highlight that the materials with a Brooks and Corey relative permeability have a higher threshold capillary pressure than the materials with a van Genuchten one.Indeed, the van Genuchten relative permeabilities imply higher Q eff w compared to the twin analyses with the   Brooks and Corey ones.The difference in Q eff w between materials with different relative permeability increases with the heterogeneity of the pore size distribution, i.e. it is larger for material A with respect to material C.
In the van Genuchten scenarios Q eff w increases with the pore size heterogeneity.This might be explained observing that, according to such formulation, the relative permeability to the non-wetting fluid at a given saturation significantly increases with the heterogeneity of the PSD (Fig. 8(b)), and it is thus larger for material A than for material C.
On the contrary, in the Brooks and Corey scenarios Q eff w decreases with the pore size heterogeneity.As with the latter formulation the relative permeability to the non-wetting fluid is not so influenced by the PSD, this is likely due to the decrease in the relative permeability to the wetting fluid with the PSD heterogeneity.
The results shown in Table 6 underline that p * c is not only a function of the capillary pressure curve of the caprock but also largely depends on the actual pore interconnectivity.This aspect might implicitly point to a further limitation of determinations of threshold capillary pressure based on the tangent method.This last method turns out to give a good estimate for materials B and C having Brooks and Corey relative permeability but on the other hand, underestimates the threshold capillary pressure of material A with Brooks and Corey relative permeability (17 % error) and overestimates the threshold capillary pressure of materials with van Genuchten relative permeability (error of 258 % for material A, 116 % for material B and 28 % for material C).
The determination of p * c with the dynamic method requires knowing the exact instant of incipient CO 2 penetration in the sample.The use of Electrical Resistivity Tomography (ERT), capable of monitoring the local saturation of the sample, might contribute to such identification, even if it makes the experimental apparatus more complex.An example of the use of ERT can be found in Cosentini et al. (2012), Koestel et al. (2008) and Musso et al. (2023).
In the practice the dynamic method is likely to be applied considering the stabilized effective discharge, obtained at larger times.This 'conventional' approach was also applied to determine the estimates that would be obtained for all the materials, relative permeability and inlet pressure (Table 2).These estimates p ∘ c are plotted in Fig. 15.When the pressure difference across the sample Δp t is below the threshold capillary pressure value, Eq. ( 5) would erroneously predict that p ∘ c = Δp t .In these cases, no CO 2 can flow through the sample and the effective discharge stabilizes to zero.Conversely for Δp t greater than p * c , the stabilized effective discharge is greater than zero and the predicted p ∘ c starts to diverge from Δp t .In the latter conditions the conventional approach leads to p ∘ c value which overestimates the actual threshold capillary pressure, as the stabilized effective discharge is lower than the one which occurs in the instant of incipient scCO 2 penetration.The errors, obtained using the stabilized value of Q eff w in Eq. ( 5) respect to p * c (Table 6), are summarised in Table 7.
The analysis of Table 7 shows that the overestimate of p ∘ c increases with the Δp t .A suggestion for experimental protocols is to check the p ∘ c Δpt ratio obtained experimentally and repeat the test applying a smaller Δp t if such ratio is smaller than 0.70, condition for which the overestimation will be approximately 15 %.
The residual capillary pressure p ∧ c identifies the pressure difference between the non-wetting pressure in the inlet vessel and the wetting pressure in the outlet vessel when the stationary conditions are reached.In all scenarios this occurs in about 7 days (see e.g. the scenarios illustrated in Fig. 14) and it confirms that the residual method is quicker than the step-by-step method and slower than the dynamic one.
The determination of the residual capillary pressure p ∧ c is easy, however it does not represent the threshold capillary pressure since the hydraulic state in the sample is heterogeneous, with the portion close to the inlet standing on the main imbibition curve and the portion close to the outlet still on the main drainage curve (Fig. 13).It is therefore not a property of the material, but is highly dependent on the experimental conditions.The heterogeneity of the hydraulic state also reveals that p ∧ c is not linked to the condition of insular saturation of the non-wetting fluid, i.e. it is not equivalent to the snap-off pressure.
The residual pressure values p ∧ c for all 108 scenarios are summarized in Fig. 16 (procedure I) and Fig. 17 (procedure II), while the error that is made with respect to p * c is reported in Tables 7-9 for each inlet pressure of the non-wetting fluid.
Similarly to what was observed for the dynamic method, the relative permeability plays a significant role in the determination of p ∧ c and with the Brooks and Corey relationship p ∧ c is always larger than the van Genuchten's one.As the van Genuchten formulation predicts a higher relative permeability for the non-wetting fluid, it allows for the penetration of a greater mass of CO 2 into the specimen, which is more than compensated by a higher compression of the water phase, which cannot flow out of the system.This implies a greater increase of p outlet w in the van Genuchten scenarios compared to the Brooks and Corey ones.As the value of p inlet nw is not significantly influenced by the relative permeability law, a higher p ∧ c is obtained with Brooks and Corey permeability.The effect of a different relative permeability increases as the heterogeneity of the pore size distribution increases.Indeed, the two formulations provide quite similar results for material C while they diverge for material A. This effect is accentuated as the inlet scCO 2 pressure increases.
When the heterogeneity of the pore size increases, the residual capillary pressure increases with the Brooks and Corey relative permeability and it decreases with the van Genucthen one.Changes in k vG r,nw  with the pore size distribution prevail with the van Genuchten formulation.The opposite occurs with the Brooks and Corey where k B&C r,nw is not so significantly influenced by the material and therefore, since the initial pressure difference across the sample is greater than the actual threshold capillary pressure(p inlet nw − p outlet w | t=0 > p * c ), scCO 2 penetrates more easily into a material with almost all pores of the same size.
From the results of the sensitivity analyses concerning the residual method, it can be appreciated that p ∧ c increases as the pressure of the non-wetting fluid imposed at the inlet p inlet nw and the size of the vessels V increase.This is due to the lower pressure drop of the non-wetting fluid at the inlet and the lower pressure growth of water at the outlet due to its compression.From Fig. 17, Tables 8-10, it is also observed that by using procedure II the underestimation of the threshold capillary pressure decreases as the holding time of the non-wetting fluid pressure at the inlet (τ) increases.For the materials studied, it is noted that for τ equal to 6 h, the p ∧ c values obtained are similar to those determined with procedure I. On the other hand, the effect of τ increases as the size of the vessels decreases.These two experimental conditions have a competing effect in determining the residual capillary pressure.In the presence of a small vessel, the effect of the compression of the fluid, due to a longer holding time of the scCO 2 pressure at the inlet, becomes more significant compared to when the volumes of the two fluids are greater and therefore the two fluids compress less for the same applied pressure.It shall be observed that generally the residual capillary pressure underestimates the actual threshold capillary pressure value, even if for material A with van Genuchten relative permeability this is not strictly verified.The reason for this different result is related to the high relative permeability to the non-wetting fluid of this material even at high water saturation, and this outcome is likely to hold more for samples with a set of microfractures rather than homogeneous ones.The best estimates were obtained with the scenarios in which the highest scCO 2 pressure at the inlet (p inlet nw = 33.09MPa), the biggest vessels V β and, in procedure II, the greatest holding time τ II were applied.Among these scenarios the underestimation ranges from 6 % for material B and relative permeability according to van Genuchten, to 25 % for material C and relative permeability according to Brooks and Corey.In light of the numerical results, it is suggested to use the residual method using a pressure of the non-wetting fluid equal to three times the approximate estimate of the threshold capillary pressure determinable with the tangent method p .. c from MIP test, big vessels and procedure I.
The longer the holding time of the non-wetting fluid pressure at the inlet, the more the results of procedure II converge to the ones of procedure I.

Conclusions
The results of the sensitivity analysis regarding the dynamic and the residual method confirm that estimates can be obtained in a limited time, ranging from about one day (dynamic method) to about one week (residual method).The results of the two methods are similarly affected by the capillary pressure curve of the sample and by its relative permeability to the wetting and the non-wetting fluid.Depending on the relative permeability to the non-wetting fluid, a greater variability in the threshold capillary pressure is expected for materials with a heterogeneous Pore Size Distribution.However, significant differences between the two method's determinations emerged.
As for the dynamic method, the main theoretical and practical observations are as follows: -the effective discharge decreases during the first few hours after the injection of the non-wetting fluid.The method assumptions hold at the instant of incipient penetration of the non-wetting fluid, which is not immediate, but it requires a previous decrease of the water pressure at the sample inlet; -the instant of incipient penetration of the non-wetting fluid depends on the applied non-wetting pressure and on the material properties.
The corresponding effective discharge is difficult to evaluate experimentally and is not the stabilised one at which conventional determinations can be made; -the use of the stabilised effective discharge in the evaluation of the threshold pressure leads to its overestimation.The overestimation increases with the pressure difference across the sample.To keep the overestimation error within 15 %, the ratio between the estimated threshold capillary pressure and the pressure drop should be greater than 0.70.
Similarly, as for the residual method: -the analyses showed that it generally underestimates the threshold capillary pressure, although this may not be true for samples with a heterogeneous pore size distribution and a high relative permeability to the non-wetting fluid (e.g. this could be the case for a microfractured sample); -the residual capillary pressure is not a property of the material, i.e. it is not equivalent to the snap-off pressure, since the hydraulic state at the end of the test is very different along the sample length; -using larger vessels and applying higher overpressures increases the value of the estimate, especially for samples with heterogeneous PSD and high relative permeability to the non-wetting fluid.The effect is relatively small when the PSD is homogeneous and the relative permeability to the non-wetting fluid is low; -maintaining a constant pressure at the sample inlet gives higher residual capillary pressure values.However, approximately the same   results can be obtained using a pulse decay method if the holding stage is sufficiently long (in the simulations this was 6 hours).

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.

Model geometries, meshes and boundary conditions
1D geometries were adopted for the simulations of both methods.In both cases, the sample had a length L = 40 mm and a diameter D = 20 mm.The mesh maximum length of the elements was 0.1 mm.
The model geometry and the initial and boundary conditions for the dynamic method are presented in Fig. 18 and Table 11.It shall be pointed out that while the mathematical formulation of a two-phase problem requires setting 2 boundary conditions (one for the wetting fluid and another one for the non-wetting fluid) at each model end, the test setup only includes a single drainage line at each sample end (as sketched in Fig. 2).This experimental arrangement does not pose any problem at the sample inlet, as the fluid at this boundary is either the wetting fluid (during the first stage of the test) or the non-wetting fluid (during the second stage of the test).In the first stage, set equal to 6 h, the boundary conditions at the inlet of the specimen were set equal to no flow for the non-wetting fluid (q nw = 0) and constant pressure p inlet w for the wetting one.These conditions were switched to no flow for the wetting fluid (q w = 0) and constant pressure p inlet nw for the non-wetting fluid.The values of p inlet nw adopted in the simulations of the different scenarios can be found in Table 2: the same values were adopted for the p inlet w for the simulations of the first 6 h.On the other hand, it cannot be excluded that at a certain time during the second stage of the test both the wetting and the non-wetting fluid could flow out of the sample.However, imposing a fixed pressure for the wetting (p w = 15 MPa) and the non-wetting fluid (p nw = 0 MPa) at the model outlet, as it is in the drainage tubes until the breakthrough of the non-wetting fluid, would result in a fixed fully water saturation state.Consequently the outflow of the non-wetting fluid would be impeded in the simulations, which is unphysical.To exclude this possibility, a dummy porous element having a length of 10 mm was introduced at the sample outlet.This dummy element was fictitiously set to have a very large porosity (ϕ = 0.99) and also an intrinsic permeability largely higher than the one of the sample (k = 1•10 − 12 m 2 ), while the relative permeabilities to both fluids were set equal to 1.The capillary pressure curve of the dummy element was calibrated so to ensure the easy escape of both fluids (α main drying = 15 MPa − 1 , n = 4, m = 0.75).At the external boundary of the dummy element, the pressure of both fluids was set equal to the initial value of the two different fluids.These extreme hydraulic properties of the dummy element, particularly the very high intrinsic permeability and the relative permeability set equal to 1, allow to reproduce the 'free' flow of both the fluid phases occurring in the drainage line without imposing any unnatural condition on the sample boundary.To reproduce the residual method, the upstream and downstream vessels were also simulated as a 'dummy' porous medium with a porosity equal to 0.99 and a very large permeability compared to the one of the sample (k = 1•10 − 12 m 2 ) (see Fig. 19 and Table 12).The relative permeabilities to both fluids were also set equal to 1.The capillary pressure curves of the vessels were calibrated so to ensure full scCO 2 saturation for the vessel at the inlet and full water saturation for the vessel at the outlet (α main drying = 15 MPa − 1 , n = 4, m = 0.75).These properties allow for a very fast equilibration of the pressure within the vessels filled with fluids.The length of the vessels "l" was calibrated so to allow the reproduction of the physical values of V α and V β provided in Table 3 (l = 3.2 mm to reproduce V α and l = 12.7 mm to reproduce V β ).To simulate procedure I, the boundary conditions at the inlet of the upstream vessel were set equal to no flow for the wetting fluid and constant pressure for the non-wetting one after a rapid increase (30 s) of the V.S. Vespo et al. inlet scCO 2 pressure from the initial value to the target (Table 2).The boundary conditions at the outlet of the downstream reservoir were set equal to no flow for both fluids.To emulate procedure II, the three stages of the procedure were considered.So in the 1 st stage, at the inlet of the upstream vessel a rapid increase (30 s) of the non-wetting fluid pressure, from the initial value to the target (Table 2), and a no flow condition for the wetting fluid was set, while at the outlet of the downstream vessel the conditions were set equal to no flow for both fluids.In the 2 nd stage, the inlet scCO 2 pressure was maintained for the target time τ (Table 4) while the no-flow condition was maintained at the inlet for water and at the outlet for both fluids.Finally, in the 3 rd stage the boundary conditions at the inlet of the upstream vessel and the outlet of the downstream vessel were set equal to no flow for both fluids.Validation of the numerical model that simulates the dynamic method The validation of the numerical model that simulates the dynamic method was done against the experimental results of a test on a caprock from the CCS research project of Ketzin (Germany) in Boulin et al. (2013).Nitrogen was used as a non-wetting fluid in this test.The caprock formation is described as consisting mostly of mudstone, clayey siltstone and anhydrite.The sample porosity was ϕ = 0.15, its diameter was D = 41 mm and its length L = 29.5 mm.The size of the numerical model was updated so to be consistent with the experimental data.The intrinsic permeability was k = 1.43 10 − 20 m 2 .The capillary pressure curve, determined through the interpretation of MIP results in Boulin et al. (2013) where the properties of the air-mercury system (θ Hg,air = 140 • and γ Hg,air = 0.485 N/m) were substituted with those of the nitrogen-water system (θ N2,H2O = 0 • and γ N2,H2O = 0.07 N/m), as suggested in Eq. ( 1).The following parameters of the Brooks and Corey equation were found to adequately fit the capillary pressure curve in

Validation of the numerical model that simulates the residual method
The model capabilities to reproduce the residual method were checked against the experimental results of a test on a Opalinus Clay Shale in The van Genuchten parameters of the capillary pressure curve of the shaly Opalinus Clay for an air-water system are provided in Favero & Laloui (2018).They were updated to the liquid CO 2 -water system, according to the expression in Eq. ( 1).The size of the vessels used in the simulations was 66 cm 3 in upstream (l = 6.9 cm) and 28 cm 3 in downstream (l = 2.9 cm).It shall be noted that this size was imposed according to the personal communication of one of the authors (Ferrari, 2023).Procedure II was adopted in the tests.The initial CO 2 pressure in the inlet vessel was raised from 8 MPa to 15.70 MPa and held for 30 min, after which the upstream drainage line was closed and the pressure started decaying.Consistently with the entrance of CO 2 in the sample, the water pressure in the downstream vessel raised starting from its initial value of 8 MPa.The comparison between the experiment and the simulation is presented in Fig.  V.S. Vespo et al.
non-wetting fluid pressure V.S. Vespo et al. the duration of the dynamic tests on claystone (k = 10 − 20 m 2 ) in Boulin et al. (

Fig. 1 .
Fig. 1.(a) capillary pressure curves adopted in this work and estimation of threshold capillary pressure with the tangent method p .. c (b) comparison between the capillary pressure curve of material B and the one of Opalinus Clay in Favero & Laloui (2018).

Fig. 3 .
Fig. 3. Pore pressure profile along the sample according to Egermann et al. (2006) assumptions in a test performed according to the dynamic method.

Fig. 6 .
Fig. 6.Determination of the residual capillary pressure p ∧ c once the pressure of the non-wetting fluid in the upstream vessel and of the water in the downstream vessel reach steady-state conditions according to procedure I (a) and II (b).

Fig. 7 .
Fig. 7. Pore Size Density curves of the samples explored in the sensitivity analysis.

Fig. 8 .
Fig. 8. Relative permeability to the wetting fluid (a) and the non-wetting one (b) according to the van Genuchten and the Brooks and Corey formulations.

Fig. 9 .
Fig. 9. Time evolution of the Q eff w Qw ratio for different values of inlet pressure (material B, Brooks & Corey permeability)

Fig. 12 .Fig. 13 .
Fig. 12. Profiles of scCO 2 saturation at the instant of maximum saturation at the inlet (blue lines) and at the end of the simulation (black lines) referring to material A with Brooks and Corey relative permeability (inlet scCO 2 pressure equal to 33.09 MPa, vessel size equal to 1 cm 3 (V α ), and holding time equal to zero (τ I ) for procedure II (continuous lines))

Fig. 14 .
Fig. 14.Residual capillary pressure p ∧ c obtained for material A with Brooks and Corey relative permeability: (a) procedure I, (b) procedure II.The dependency of the time evolution of wetting and non-wetting pressures in the vessels on the inlet pressure p inlet nw , the vessel size V and the holding time of the non-wetting pressure τ are

Fig. 15 .
Fig. 15.Abacus of estimated threshold capillary pressure p ∘c for all the materials, the relative permeability laws and inlet pressure (dynamic method).

Fig. 16 .
Fig. 16.Abacus of residual capillary pressure p ∧ c for all the materials, the inlet pressures, the Brooks and permeability, (b) the van Genuchten permeability with procedure I (residual method)

Fig. 17 .
Fig. 17.Abacus of residual capillary pressure p ∧ c for all the materials, the inlet pressures, (a) the Brooks and Corey permeability and holding time equal to zero τ I , (b) the van Genuchten permeability and holding time equal to zero τ I , (c) the Brooks and Corey permeability and holding time equal to 6 h τ II , (d) the van Genuchten permeability and holding time equal to 6 h τ II with procedure II (residual method)

Fig. 18 .
Fig. 18.Model for the simulation of the dynamic tests.
Fig. 20.Comparison between the experimental data of the dynamic test on the Ketzel sample in Boulin et al. (2013) and model predictions.

Fig. 21 .
Fig. 21.Comparison between the experimental data of the residual test on Opalinus clay sample in Minardi et al. (2021) and model predictions.

Table 2
Non-wetting fluid pressure used in the sensitivity analyses.

Table 6
Threshold capillary pressure p * c of all scenarios studied in the sensitivity analyses.

Table 11
Initial and boundary conditions for the dynamic tests.

Table 12
Initial and boundary conditions for the residual tests.