Field-scale impacts of long-term wettability alteration in geological CO$_2$ storage

Constitutive functions that govern macroscale capillary pressure and relative permeability are central in constraining both storage efficiency and sealing properties of CO$_2$ storage systems. Constitutive functions for porous systems are in part determined by wettability, which is a pore-scale phenomenon that influences macroscale displacement. While wettability of saline aquifers and caprocks are assumed to remain water-wet when CO$_2$ is injected, there is recent evidence of contact angle change due to long-term CO$_2$ exposure. Weakening of capillary forces alters the saturation functions dynamically over time. Recently, new dynamic models were developed for saturation functions that capture the impact of wettability alteration (WA) due to long-term CO$_2$ exposure. In this paper, these functions are implemented into a two-phase two-component simulator to study long-term WA dynamics for field-scale CO$_2$ storage. We simulate WA effects on horizontal migration patterns under injection and buoyancy-driven migration in the caprock. We characterize the behavior of each scenario for different flow regimes. Our results show the impact on storage efficiency can be described by the capillary number, while vertical leakage can be scaled by caprock sealing parameters. Scaling models for CO$_2$ migration into the caprock show that long-term WA poses little risk to CO$_2$ containment over relevant timescales.


Introduction
Geological CO 2 storage can be successfully implemented in deep saline aquifers that possess favorable characteristics [IPCC, 2005], e.g., the storage formation has sufficient capacity, injectivity, and containment properties to store the desired quantity of CO 2 at the required injection rate. A competent caprock should be verified that provides a capillary and permeability seal to prevent buoyant CO 2 from migrating upwards and is free of defects that could leak or seep CO 2 out of the storage formation [Busch et al., 2008]. In addition, lateral migration of CO 2 should stay within the defined boundaries of the storage reservoir. In some cases, the boundaries are defined by a structural trap. In open stores with no defined trap, sufficient CO 2 trapping efficiency in residual and dissolved phase is needed to prevent unwanted lateral CO 2 migration [Tucker, 2018].
There is a wide array of geological and operational factors affecting CO 2 storage migration and containment in deep, saline aquifer [Birkholzer et al., 2015]. The physical and chemical rock-fluid interactions play an important role through their impact on the constitutive functions that govern macroscale capillary pressure and relative permeability. These functions, often called flow functions or saturation functions, govern the displacement processes in a multiphase flow setting and are central in constraining both storage capacity of the formation and sealing properties of the caprock [Mori et al., 2015, Negara et al., 2011, Oostrom et al., 2016, Sarkarfarshi et al., 2014.
Capillary pressure and relative permeability functions for the storage reservoir and caprock are strongly influenced by wettability. The term wettability refers to the preference of one fluid over the other for the rock surface [Bonn et al., 2009, Eral et al., 2013 and is a pore-scale phenomenon that has direct influence on macroscale displacement processes [Anderson, 1987, Bobek et al., 1958, Falode and Manuel, 2014.
Wettability has been shown to play a key role in geological CO 2 storage with regard to CO 2 migration, trapping, and sealing capacity. Saline aquifers are brine-bearing formations, and thus often assumed to be strongly water-wet. However, differences in initial wetting state can be caused by a number of factors that affect the surface chemistry of rock minerals, including temperature and pressure conditions, brine salinity, mineralogy, and organic content (see Iglauer et al. [2014] and references therein). Wettability significantly affects the capillarity between CO 2 and brine, which in turn impacts the efficiency of brine displacement. Studies have shown that wettability affects plume migration [Al-Khdheeawi et al., 2017a] and residual trapping [Krevor et al., 2015]. In terms of containment, the capillary seal of the caprock barrier is reliant on a strongly water-wet condition to maintain a high capillary force in combination with smaller pore radii of a tight medium. An effective seal increases storage capacity in structural traps by supporting a thick column of CO 2 [Iglauer et al., 2014], which can be compromised if the seal weakens over time due to changes in wettability [Rezaee et al., 2017].
Although wettability is often considered a static property, there is a possibility for wettability to change dynamically over time. This process of wettability alteration (WA) is a dynamic process, that involves a complex interaction of surface mineralogy, fluid composition, and reservoir conditions that is percipitated by the introduction of a wettability-altering agent [Bonn et al., 2009, Buckley et al., 1988. WA has been long studied in the petroleum sector, and a variety of methods and agents, that include solvents such as supercritical CO 2 , have been applied to enhance WA in oil reservoirs to the benefits of increased recovery [Drexler et al., 2020, Sun et al., 2017. Modeling of WA often assumes an instantaneous change of wettability as a function of solvent concentration [Lashgari et al., 2016]. However, there are some cases where a WA process depends on exposure time that can extend for weeks or months to the wettability-altering agent [Blunt, 2017, Farokhpoor et al., 2013, Fatah et al., 2021, Jafari and Jung, 2016, Saraji et al., 2013.
For CO 2 storage, there has been some intriguing evidence of long-term WA caused by CO 2 exposure that significantly weakens the strength of capillary forces of a water-wet medium, which can apply to both reservoir rocks and caprock material. For example, core-flooding experiments and capillary pressure measurements on storage reservoir samples have reported a instability or slow reduction in capillary pressure-saturation relationships over time [Farokhpoor et al., 2013, Plug and Bruining, 2007, Tokunaga et al., 2013, Wang and Tokunaga, 2015, Wang et al., 2013. In particular, we note in Wang et al. [2013] that wettability was measured through contact angle (CA) change in silicate and carbonate sands from strongly water wet to intermediate wet over the course of weeks and months. Studies on shales exposed to CO 2 over long periods show a steady increase in shale/water CA over days and months [Chiquet et al., 2007a, Fatah et al., 2021, with trends that indicate additional exposure would lead to further CA increase. Similar studies using reservoir rocks also show an evolution from water-to CO 2 -wet conditions due CO 2 exposure in core-flooding experiments [Fauziah et al., 2021, Valle et al., 2018. All these bench-scale experiments indicate that WA introduces long-term dynamics into the capillary pressure and relative permeability functions. Such dynamics could have a significant impact on CO 2 storage at the field scale in ways that may be positive or negative. On the one hand, a caprock that becomes intermediate-wet or hydrophobic is less effective than a water-wet caprock in providing an effective capillary barrier. On the other hand, capillary diffusion may be reduced in the storage reservoir through long-term WA, which can be beneficial to injectivity and capacity by causing the CO 2 to displace the resident fluid efficiently [Kassa et al., 2021a]. A wettability change, in general, affects the space-time evolution of CO 2 plume in the formation in ways that has not been fully understood yet.
There have been several attempts to understand the impact of static wettability on long-term CO 2 storage at field scales. In Al-Khdheeawi et al. [2017a, 2018, 2017b the main conclusions are that CO 2 vertical migration is enhanced with more CO 2 -wet systems, which is due to the lower residual trapping capacity of CO 2 -wet rocks.
The impact is stronger when considering heterogeneous wettability and increased temperature. These studies did not investigate the impact of CO 2 exposure that changes rock wettability over time. Another study investigated wettability alteration caused by increased reservoir temperature under CO 2 injection. The relative permeability curves were set in predetermined zones different temperature, but the curves were not adjusted dynamically in time [Abbaszadeh et al., 2020].
Understanding how CO 2 storage is impacted by dynamic wettability caused by CO 2 exposure requires a different approach than previously used to study static wettability. Two simulation components are needed: (a) flow functions that capture the dynamics due to CO 2 exposure, and (b) analysis and characterization of long-term WA impacts through implementation of dynamic flow functions into field-scale simulation. We have recently addressed the first aspect in Kassa et al. [2020] and Kassa et al. [2021b] where we successfully developed dynamic capillary pressure and relative permeability functions through mechanistic modeling combined with upscaling of wettability dynamics from pore to core scale. The developed models are extended forms of standard constitutive models either by interpolation between two end-state capillary pressure functions [Kassa et al., 2020] or by incorporating WA dynamics directly into the parameters of the existing relative permeability model [Kassa et al., 2021b]. The dynamics in both models are driven by exposure time to CO 2 which is easily calculated from the local saturation history. We found that both models compare well against pore-scale simulation data where WA is modeled in individual pores. More interestingly, we determined that the upscaled model parameters have a clear relationship with the underlying pore-scale WA model parameter. Moreover, the developed constitutive models are robust and lend themselves to implementation within standard flow simulators.
To our knowledge, numerical simulation of the field-scale impacts of WA caused by long-term CO 2 exposure has not been performed previously. In this paper, we implement the aforementioned dynamic relative permeability and capillary pressure functions [Kassa et al., 2020[Kassa et al., , 2021b into flow models for CO 2 storage applications in the open porous media (OPM) framework [Rasmussen et al., 2019]. We investigate the impact of long-term WA process for CO 2 storage with a focus on the migration patterns and containment in a saline aquifer that involves both  Figure 1: Model development involved in preparation for long-term wettability alteration field-scale studies. The boxed content illustrates key features of previous work [Kassa et al., 2020[Kassa et al., , 2021b implemented in this study. Further discussion of the "Modeling" formulations can found in Section 2.1.1.
viscous-driven flow in the storage aquifer and the potential for buoyancy flow through the caprock. Therefore, we concentrate our analysis on one-dimensional horizontal and vertical flow systems in order to understand the impact of dynamic changes in wettability on both flow regimes separately. For the horizontal case, we study the impact of WA on storage efficiency given by the front location, while for the vertical case, we focus on the impact of CO 2 containment under a capillary seal (i.e., caprock). In both cases, we study the sensitivity of wettability dynamics and subsequent field-scale impacts under different conditions (i.e., permeability, entry pressure, injection rate, and porosity) and ultimately seek an appropriate scaling relationship for the quantity of interest in order to generalize our findings. The quantified behavior contributes to increased understanding of the efficiency and integrity of long-term CO 2 storage subject to WA caused by CO 2 exposure.

Methodology
The goal of this study is to investigate the impacts of dynamic changes in wettability caused by long-term CO 2 exposure on field-scale CO 2 plume evolution and containment. In this section, we outline the methods that have been implemented to carry out the desired field-scale simulations. Figure 1 shows the relevant components of the dynamic saturation functions developed in previous work that are linked to the implementation and field-scale simulation study carried out in this work.
As discussed in Section 1, laboratory observations indicate that CO 2 exposure can induce chemical alteration of the native rock. Assuming saline aquifers are initially water-wet, we consider that long-term CO 2 exposure on the order of months can shift wettability to intermediate wet or even CO 2 wet, where each wetting state is defined by the contact angle (CA) between water and solid as shown in the left-most illustrations in Figure 1. As a result, the saturation functions are altered continuously in time, which is shown schematically for both the phase relative permeability and capillary pressure functions (center Figure 1). Alteration continues until the wetting condition has reached a final state under maximum exposure.
When studying this phenomena at the field scale, we expect that dynamics in saturation functions will be inherently linked to local CO 2 exposure, which itself is controlled by migration of CO 2 at the field scale. As a result, different points of the aquifer will have a different local history of CO 2 exposure, which leads to local alteration of saturation functions that varies in space and time throughout the storage site. This spatial/temporal complexity in capillary pressure and relative permeability can in turn impact field-scale CO 2 migration in ways that are not easily predicted a priori. Field-scale studies are needed to quantify the impact of this interplay between dynamics in saturation functions and large-scale flow processes, illustrated in the right-most panels in Figure 1.
The link between CO 2 exposure and time-dependent saturation functions requires a model to evolve capillary pressure and relative permeability functions as the wetting condition continuously changes. As depicted under "Modeling" in Figure 1, we have previously developed new formulations for capillary pressure and relative permeability that incorporate exposure to CO 2 as an additional dynamic variable [Kassa et al., 2020[Kassa et al., , 2021b. These dynamic models are implemented in a standard field-scale simulation approach (described in Section 2.1) that includes a brief overview of the dynamic model development from previous work in Section 2.1.1.

Field-scale governing laws
CO 2 storage in saline aquifers is modeled as a three-dimensional flow system consisting of two immiscible fluids, CO 2 and brine, that are mutually soluble fluids. In this section, we introduce the model equations to describe the TPTC flow in a porous medium. The general form of the TPTC model is presented below, however we introduce some simplifications later during the model implementation (cf. Section 2.2). For the sake of brevity, we consider constant temperature T and salinity in this study, but these effects can be included into the model in a straightforward manner.
Let Ω ∈ R d , d = 1, 2, or 3 be a permeable domain having a Lipschitz continuous boundary ∂Ω and saturated with two-phase fluids that have two components each. The phases, wetting and non-wetting, are indicated by subscripts α ∈ {w, n} and the two components are represented by index k ∈ {CO 2 , water}. The mass conservation law of component k is described by: where X k α is component k mass fraction in phase α, φ denotes the porosity of the porous domain, S α is the phase α saturation, ρ α is the density of phase α, u α is the Darcy flux of phase α, and F k is the source term of component k.
For each phase α ∈ {w, n}, the Darcy flux u α is given by the multiphase extension of Darcy's law: where K : Ω → R d×d is the intrinsic permeability tensor of the rock, k rα is the relative permeability for phase α, µ α is phase viscosity, P α is phase pressure, and g is the gravitational vector. Hereafter, we will use K and K for a tensor and scalar permeability, respectively. The diffusive flux of component k in phase α, j k α , is represented by the Fick's law and has the form of: where D k α is the molecular diffusion coefficient and τ α is tortuosity for phase α which is calculated as [Millington and Quirk, 1961]: The component mass fractions, X k α , satisfies: The phase α saturation takes values between zero and one, and the sum of phase saturation equals one. That is: The pressures, wetting and non-wetting, are connected by the capillary pressure relation as A model for mass transfer between the phases is needed to close the system of equations. In the general form, we consider the fugacity constraints where they express the requirement that non-wetting (gas) and wetting (liquid) fugacities have to be equal for each component [Coats, 1980, Voskov andTchelepi, 2012]. The fugacity coefficients are calculated according to the work of Spycher and Pruess [2005]. We note that other forms of mass transfer may be considered, including equilibrium partitioning. If mass transfer is not modeled, then this closure relation may be omitted from the model formulation.
Equations (1) through (8) with appropriate initial and boundary conditions can be used to describe TPTC flow dynamics in a porous medium. These equations provide a complete description of the physics of isothermal TPTC flow in a porous medium bearing in mind that the relative permeabilities and capillary pressure parameterizations are given.

Previous work: Dynamic saturation functions
Capillary pressure in Equation (7) and phase relative permeability in Equation (2)  The dynamic constitutive models in Kassa et al. [2020] and Kassa et al. [2021b] were developed by applying an upscaling workflow that involves correlating Darcy-scale models to dynamic capillary pressure and relative permeability data generated from pore-scale numerical experiments. The key features to this workflow are: (A1) A mechanistic model for CA change from an initial to final wetting state that is incorporated directly in pore-scale simulations.
(A2) A dynamic model formulation identified at the Darcy scale that incorporates exposure time and can be correlated to the dynamic data with as few parameters as possible.
(A3) A link between the correlated dynamic Darcy-scale parameter(s) and the pore-scale model for CA change.
We note that the above approach does not capture changes in residual saturation for CO 2 and brine as a result of different wetting conditions. Although the connection between residual saturation and wetting condition is well studied, incorporation into a dynamic framework is still a subject of ongoing work. As such, the dynamic models presented are limited, but still provide important insight into CO 2 migration patterns due to changing wetting conditions.
The first step in the previously published upscaling workflow is the pore-scale mechanistic model for CA change.
The chosen model was motivated by experimental evidence that indicates a gradual and permanent change in the fluid-fluid CA when exposed to a WA agent, in this case CO 2 , over a long period of time. CA change is caused by adsorption of the WA agent to the pore surface that subsequently alters surface chemistry and changes the affinity of fluids to the solid. In addition, a longer exposure time induces a greater change in CA. We thus designed a dynamic CA model that is sorption-based and time-dependent and has the form (see Kassa et al. [2020] for development details): where θ i and θ f are the initial and final CA respectively, C is a model parameter that determines the rate of WA over time, with increasing C indicating slower CA change, and χ is a cumulative measure of exposure time to the WA agent at the pore level. The quantity χ was evaluated by integrating the chosen measure of local exposure, i.e., agent concentration or saturation, over time. As χ is a quantity that always increases under exposure or remains constant under the absence of the WA agent, the above model leads to irreversible CA change.
With the above mechanistic WA model incorporated at the pore scale, a set of numerical experiments were performed to simulate laboratory measurements of capillary pressure and relative permeability, i.e., subsequent drainage and imbibition cycles with stepwise changes in saturation over time. The simulated data were obtained by modeling the pore scale as a bundle of cylindrical or triangular tubes. A schematic illustration of capillary pressure and relative permeability data obtained in our previous studies are shown in Figure 1.
Analysis of the simulated data showed that both saturation functions evolve smoothly from an initial to final state as the CA changes dynamically and heterogeneously with increasing exposure time to the non-wetting fluid.
In the reported experiments, CA changed from a strongly water-wet system (θ i = 0 • ) to an intermediate-wet system (θ f = 85 • ), although the opposite WA process could equally have been simulated. Figure 1 shows that capillary gradually decreases over a period of months to years over several subsequent drainage and imbibition processes. Similarly, relative permeability alteration causes an overall decrease in non-wetting phase permeability with increasing exposure time compared with an increase in wetting phase permeability for the same exposure.
This temporal component motivated a new definition for upscaled exposure time at the Darcy-scale, given as: where t ch is a pre-specified characteristic time which is used to scale the history of exposure, and X WA is the chosen measure of exposure of WA agent. Exposure may be either modeled as dissolved mass fraction X WA = X CO2 w or as fluid saturation X WA = S n , depending on whether the WA agent is a dissolved solvent or invading fluid, respectively. We note that χ is a dimensionless non-linear parameter that depends solely on CO 2 concentration or saturation and time, which are all readily available model variables in Darcy-scale simulators.
The above definition of exposure time was then used to formulate the dynamic models for capillarity and relative permeability. In Kassa et al. [2020], we found that alteration of capillary pressure could be modeled as an interpolation between two static states, where the interpolation coefficient is coupled to macroscale exposure time, χ. The resulting formulation is an interpolation-based dynamic capillary pressure model as follows: where P i c and P f c are the capillary pressure functions at the initial and final wetting states, respectively, χ is the exposure time at the Darcy scale, β is a fitting parameter, and S ew is effective saturation which can be defined as: where S rw and S rn are wetting and non-wetting residual saturations, respectively. The end wetting-state capillary pressure functions P i c and P f c are modeled as static functions that can be represented by a Brooks-Corey model: where 1/λ is a parameter that represents pore-size distribution, c i and c f are entry pressures for initial and final wetting-state conditions, respectively.
Similar analysis in Kassa et al. [2021b] showed that relative permeability dynamics could be more efficiently modeled by including the exposure time directly into the function parameters rather than following the interpolation approach used above. Furthermore, we found a reduced LET relative permeability model is preferable due to its flexibility. More specifically, from the original relative permeability models in Lomeland et al. [2005] (six fitting parameters) we have obtained reduced models (two fitting parameters) based on assumptions in the pore geometry, and we have added a dynamic function to capture the wettability dynamics in the relative permeability-saturation curves. This approach results in a dynamic model for relative permeability as follows: where the dynamic function F(χ) is given by E i and E f are empirical data fitting parameters for the initial and final wetting-state condition and η is a dynamic fitting parameter that controls the WA induced dynamics in the relative permeabilities.
The reader is referred to the cited papers for more details on the comparison between the simulated data and the correlated functions. In summary, the correlated models perform well for any general saturation history resulting in an excellent match with the simulated data, thus demonstrating the robustness of the dynamic models.
We emphasize that the above presentation is intended to provide background on the approach used to understand and model the manifestation of pore-scale WA at the field scale. In this study, only the resulting macroscale dynamic flow functions, Equations (11) and (15), and associated parameters are applied further in field-scale simulations. CA change is not explicitly modeled in these simulations, i.e., θ is not a variable at the field scale, but the mechanisms of CA change described in Equation 9 are implicitly captured in the dynamic flow functions. This implicit connection between pore scale and macroscale is described further in the next section.

Dynamic model implementation and parameterization
The above dynamic models in Equations (11) and (15) with associated parameters are substituted directly into the TPTC model in place of the usual static functions, i.e., capillary pressure relation in Equation (7)  There are additional parameters in both Equations (11) and (15) which need to be obtained by characterization of the rock-fluid system of interest. For instance, capillary pressures at end wetting states needed in (11) can be obtained from laboratory analysis, possibly by rescaling curves obtained by nitrogen or mercury drainage experiments. For relative permeability, the initial wetting state needs to be fully characterized using a non-reactive fluid.
The dynamic β and η parameters are obtained by characterization of WA, either by fitting to laboratory or numerical pore-scale experiments performed for this purpose. In our analyses in Kassa et al. [2020] and Kassa et al. [2021b], we found a direct connection between the dynamic parameters and the underlying pore-scale WA model for CA change. In both cases, β and η have a correlated relationship with the pore-scale parameter C in Equation (9). According to Kassa et al. [2020], the parameter β is represented as a power function of the pore-scale parameter C: where b 1 > 0 and b 2 > 0 are fitting parameters and are determined from WA experiments. As C decreases, indicating faster WA, β decreases. If C = 0, wettability changes from the initial to final wetting state instantaneously.
On the other hand in Kassa et al. [2021b], the parameter η is related linearly, but inversely, with parameter C as: where ν 1 > 0 and ν 2 > 0 are empirical fitting parameters. Faster WA (decreasing C) results in increasing η, reaching a maximum value, ν 2 , when C = 0.

Implementation for macroscale simulation
There are numbers of powerful numerical porous media simulators for TPTC flow scenarios. For instance, DuMux [Flemisch et al., 2011], MRST [Lie, 2019], and OPM [Rasmussen et al., 2019] are among the software that are used to solve multi-phase flow models in a porous medium. These simulators solve the flow problem by reducing the continuous system of PDEs to an algebraic system of equations. In this paper, we employed the OPM framework to solve and analyze the TPTC model above.
OPM Flow is a reservoir simulator available in OPM capable of performing TPTC field-scale simulations.
This simulator implements a reduced version of the presented TPTC model, where the components in the phases are computed by equilibrium partitioning and molecular diffusion is not included. This simplifications result in faster simulations for advection-dominated problems (e.g., continuous well injection of CO 2 ). In studies where it is necessary to include molecular diffusion (e.g., CO 2 migration into caprock), then the complete TPTC model implemented in opm-models can be used, which model the mass transfer between phases by fugacity constraints.
We refer to the latter simulator as TPTC simulator.
Originally, both OPM Flow and TPTC simulators were developed for flow systems with standard saturation functions under the assumption of static and uniform wetting condition in space and time. We have implemented the dynamic saturation functions in both simulators. This is done by adding the dynamic saturation functions (Equations (11) and (15)) and a new variableχ (Equation 10) to the simulator, whose value is approximated by the weighted (time step ∆t over characteristic time t ch ) cumulative sum of the measure of exposure of WA agent The OPM version which is used in this study is based on the 2021-04 Release. The two point flux approximation (TPFA) and backward Euler (BE) method are used to discretize the model in space and time, respectively. The resulting system of equations is linearized using the Newton-Raphson method. The simulator uses automaticdifferentiation (AD) to calculate the Jacobian of the system. We refer to Rasmussen et al. [2019] for details on the implementation and model setup of Flow. The links to download the corresponding code to reproduce the numerical studies can be found at the end of the manuscript.

Numerical experiments and results
In this section, we apply the model and implementation in Section 2 to examine the effect of long-term WA on CO 2 plume migration in the reservoir and containment by the caprock. We consider a change in wettability between two wetting states that is caused by exposure to WA agent, which results in alteration of the saturation functions as described earlier. To understand the important factors affecting long-term WA, we consider different parameters that characterize the fluid-rock system, including intrinsic properties, K, φ, c i , flow rate q, in combination with different values of the WA dynamic coefficients, β(C), and η(C).
The impact is investigated by simulating CO 2 migration under WA for two simple systems: (1) one-dimensional horizontal flow system (1D-H) of constant CO 2 injection and (2) one-dimensional vertical flow system (1D-V) of CO 2 placed under a caprock that is initially sealing. These examples give insight into the impact of WA on storage efficiency and storage containment. Finally, scaling relationships are proposed to understand the parameters for which the impact of WA is significant.

End-wetting state saturation functions
To start, we present the saturation functions applied in the reservoir for the initial-and final-wetting states.  [Kassa et al., 2020[Kassa et al., , 2021b. However, any end-state saturation functions can be used if appropriate field or lab data are available.    Table  1 and q = 10 −3 m 3 /s).
We observe in Figure 2a that WA reduces the entry pressure to the non-wetting phase by two orders of magnitude, causing the capillary pressure function to shift downward while maintaining the same curvature. For relative permeability (Figure 2b), WA dramatically decreases CO 2 mobility while having an opposite impact on the water mobility. The reduction of CO 2 relative permeability is explained by CO 2 preferring smaller pores with reduced entry pressure, which reduces the ease at which CO 2 can flow. On the other hand, larger pores become the preferred flow path for water, leading to higher water relative permeability with intermediate-wet/hydrophobic conditions. We note that the relative permeability curves in Figure 2b are consistent with our previous investigation using a bundle-of-tubes approach to upscale the dynamics in saturation functions. These curves exhibit a concavity in the CO 2 relative permeability under water-wet conditions that is not usually observed in the lab. However, we choose these curves in order to explore the impact of a large change in relative permeability curves due to WA even if a realistic system may have a much more narrow alteration.
We further analyze the static functions in order to understand the flow system before adding dynamics. To this end, it is useful to describe the predominant mechanisms affecting horizontal and vertical flow by way of the fractional flow function (neglecting capillary pressure): where λ α are the phase mobilities defined as k rα /µ α , ∆ρ = ρ w − ρ n , q the total flow rate, and ϑ is the flow inclination angle. Note that the fractional flow in Equation (18) Table 1 and depicted in Figures 2c and 2d. We make several observations based on Buckley-Leverett analysis, which is conducted under the assumption of zero capillary pressure for simplicity. When assuming a constant total flow in one dimension, neither end-wetting state will develop a shock as there is no inflection point in the flow functions. The development of a shock would be expected if the chosen relative permeability functions were convex as often observed for CO 2 -brine systems per the discussion above. However, we note that the CO 2 -water front will always advance further for the initial-wet conditions compared to the final-wet conditions, which applies for both horizontal and vertical flow.
As these curves neglect capillarity, their saturation profile under horizontal and vertical flow will be affected by capillary pressure. Capillarity will smooth saturation shocks and retard the displacing fluid front, and this impact will be greater under initial wet conditions given a stronger capillarity. We also note that constant total velocity is also quite restrictive, but nevertheless gives important insight into flow behavior due to changes in relative permeability and provides a good complement to numerical simulations.

1D-horizontal (1D-H) flow system
In this example, we consider a 1600-m homogeneous 1D flow system discretized into 640 equal-sized cells. The system is initially saturated with brine, and CO 2 is injected from the left boundary at a constant rate for one year. The saturation functions are dynamically varied according to the end-wetting states with parameters given previously in Table 1. In this section, we perform all simulations using the dynamic wettability implementation in the OPM Flow simulator. Here, we neglect the effect of molecular diffusion term and consider the WA agent X WA to be the non-wetting saturation S n .
In Kassa et al. [2020Kassa et al. [ , 2021b, we have found mathematical relationships between the pore-scale parameter C and the dynamic parameters β and η given in Equations (16) and (17). Based on Kassa et al. [2020Kassa et al. [ , 2021b, for the simulations we consider the pore-scale parameter C in the interval [10 −5 , 10 −4 ]. The parameters b 1 , b 2 , ν 1 , and ν 2 are estimated by comparing to WA experiments. In this work, we set the values of these parameters to study their dynamical impact in a time scale of order of months (t ch = 10 7 s). For the simulations in the 1D-H study, these values are shown in Table 2. Then, the combination of values β and η are uniquely determined by the value of C. We remark that we use the value of the exponent parameter b 2 as calibrated in Kassa et al. [2021b]. Table 2: Parameters describing the relationship between the pore-scale parameter C associated with CA change and dynamic coefficients β and η for the 1D-H studies (t ch = 10 7 s, [−] 10 −5 1 45 2.5×10 −5 5.2 37.5 10 9 1.8 4.999 × 10 5 50 5 × 10 −5 18.1 25 7.5×10 −5 37.6 12.5 10 −4 63.1 10 −2

1D-H base-case scenario study
We begin with a base-case scenario where CO 2 is injected for one year for a fixed injection rate, permeability, and porosity, namely q = 10 −7 m 3 /s, K = 10 −10 m 2 , and φ = 0.1. First, we examine the base-case scenario under static wetting conditions at the two end states. Figure 3a shows CO 2 migrates far into the domain at low saturation when both functions are in their initial state. Conversely, CO 2 migrates more slowly and fills more of the pore space behind the front when both functions are in their final-wetting states. These results reflect the fractional flow characteristics observed earlier in Figure 2, but now with the added effect of capillarity. We also observe an independent influence of the wetting state of each saturation function. For instance, we observe the impact of different wetting-state relative permeability functions is more significant when the P c function is weaker (final-wetting) than for the initial-wetting P c function.
Next, we investigate the isolated impact of wetting dynamics in capillary pressure by varying only the dynamic WA parameter β(C) according to the values in Table 2 for the base case scenario. Here, the relative permeability model is kept static at the initial-wetting state. Figure 3b shows that a small C is required in order for the dynamics in P c to have an observable impact during the one-year simulated time period, where we recall that smaller C implies faster WA dynamics in P c . In addition, the power-law relationship β(C) results in a non-linear dependence on C. An important observation can be made regarding the comparison between fastest WA case C = 10 −5 and the reference static case using the initial-wetting k r and final-wetting P c . One would expect that the dynamic case should be very similar to the static case since the WA dynamics are fast, but here we see that the CO 2 front migrates more slowly and builds to a higher saturation when dynamics are included. This difference indicates that the development of heterogeneous wettability along the horizontal column created by WA dynamics affects CO 2 flow in a complex way that can not be predicted by static wetting simulations alone.
For the isolated impact of dynamics in relative permeability (Figure 3c), we test the same values of C that result in the values of η listed in Table 2. In these simulations, the capillary pressure model is kept static at the initial-wetting state. Since the model for η(C) is linear with C, the resulting change in CO 2 migration is more gradual with C than observed for dynamics in capillarity alone. We also observe that the fast-dynamics case (C = 10 −5 ) does not match the reference static end-state case (final-wetting k r , initial-wetting P c ). The CO 2 front in the dynamic case approaches the reference static case towards the inlet (where the wettability has mostly reached the end state), but the front is significantly more advanced where the wettability is still in a transition between initial and final states.
Combining the impact in both β and η simultaneously (Figure 3d), the results show that the most C cases with slower wettability dynamics remain close to the reference case with static wettability (initial k r , initial P c ),  Table 2.
which shows the impact of WA dynamics is relatively low given the base-case parameters. The saturation at the inlet begins to build evenly for increasing dynamics reflecting the altered wetting state there, but it is only for the fastest WA dynamics C = 10 −5 where we observe a significant alteration of the location of the CO 2 front. Here again, there is a discrepancy between the fast dynamic case and the final reference static case (final k r , final P c ), showing the complex flow behavior when wettability is varying in space and time along the flow path. Figure 4 shows the spatial evolution of the upscaled exposure time χ and the non-wetting saturation S n for the dynamic case C = 10 −5 in Figure 3d at three different injection times. We observe that after 30 days, the evolution of the dynamic S n profile is practically the same as the static one. This is expected as the wettability alteration effects start to significantly impact the system after a few months, as observed on the saturation profiles at 180 days and one year respectively. These figures confirm the complex flow behavior when fast dynamics are presented.

1D-H capillary scaling
We now compare simulations where rock properties (K, φ), injection rate (q), and the pore-scale parameter (C) are considered to vary independently, as given in Table 3. For the purpose of quantifying the impact of WA dynamics in P c and k rα compared to a static-wet system, we define the scaled front location difference, SFLD, as: where x i and x dy denote the CO 2 front location using the static initial-wetting functions and dynamic-wetting functions, respectively. The values of x i and x dy are considered as the further spatial location from the injection well of the non-wetting saturation value greater than a threshold value (here 10 −4 ). We avoid x i reaching the boundary by considering smaller values of q, K, and larger values of φ as observed in Table 3. Table 3: Tested parameter values used in the 1D-H system. A set of simulations were performed for each individual value of q, K, and φ in which C is varied, where the combined impact on dynamic wettability functions is studied. For each value of C, the corresponding value of β and η is given in Table 2.
[−] 10 −7 10 −10 0.1 10 −5 5×10 −8 5 × 10 −11 0.2 2.5 × 10 −5 2.5×10 −8 2.5 × 10 −11 0.4 5 × 10 −5 7.5 × 10 −5 10 −4 In addition, we employ a macroscale definition of capillary number, Ca, in order to characterize the impact of WA with respect to the flow regime for each parameter combination (Table 3). Following Armstrong et al. [2014], we compute Ca on the injection side using the initial entry pressure c i as: where ∆x is the length of the grid cell at the entry, A the transversal area, and other parameters are defined previously. The transversal area for these simulations is A=1 m 2 . We perform a similar sequence of numerical experiments, first isolating the impact of dynamics in capillarity and relative permeability by allowing dynamics in one function while keeping the opposite function static in the initial-wetting state. Then, we perform experiments of the combined dynamics. The capillary number for these simulations over more than an order of magnitude from 10 −5 to over 10 −4 . For reference, the base-case scenario described above has a capillary number of approximately 5 × 10 −5 .  Table 3) under isolated dynamics in saturation functions: (a) capillary pressure and (b) relative permeability.
We observe the relationship between SFLD and Ca converges onto a single curve for a given C when considering the isolated dynamics in capillary pressure ( Figure 5a) and relative permeability (Figure 5b). This implies that the flow regime influences the impact of WA dynamics. Noting the scale difference on the y-axis for both plots, it is evident that dynamics in P c have a much greater effect than dynamics in k r especially at low capillary number, slowing the CO 2 front migration by as much as 35%. When C is large (fast WA dynamics) where there is only a 7% reduction in front location. However, the influence of capillarity diminishes more quickly with higher Ca and higher C, and the relative permeability dynamics have a slightly larger impact that is sustained for higher capillary number.
For impact of Ca on the combined dynamics shown in Figure 6, we observe an increased impact at low capillary number when dynamics are modeled in both saturation functions. For example, the impact of WA dynamics in SFLD doubles for C = 2.5 × 10 −5 when relative permeability dynamics are added to capillary dynamics, with the increase being greater at higher capillary number. Thus, relative permeability dynamics help to compensate for the disappearing impact of dynamic P c at higher capillary number. However, the impact of the relative permeability dynamics is very minor, which implies that using less extreme relative permeability curves more reflective of  Table 3. laboratory studies will have negligible impact on CO 2 migration under dynamic wettability conditions. Intuitively, the results in Figures 5 and 6 reflect the fact that very viscous flows will not be greatly impacted by wettability dynamics in either saturation function. The minimal effect of capillarity at high capillary number is expected given that capillary effects disappear with higher Ca, and therefore any additional dynamics in P c have a negligible effect. But we also see that the impact of relative permeability dynamics also decreases with higher capillary number albeit more slowly. For CO 2 storage settings, these observations point to the increased importance of wettability dynamics farther from the injection well where capillary forces are likely to be dominant.

1D-vertical (1D-V) flow system
In this section, we consider a 1D vertical cross-section of an aquifer-caprock system. This system is employed to demonstrate the impact of WA on containment of CO 2 beneath an initially sealing caprock. This system is modeled as a one-dimensional flow domain of A = 10 m 2 cross-sectional area by H = 100 m height. The system is discretized into 100 elements along the height. The top and bottom zones have contrasting permeability and capillary entry pressure. The permeability and porosity in the aquifer are homogeneous and fixed at 10 −10 m 2 and 0.2, respectively. The caprock permeability is homogeneous and varied in different simulations (see Table 4), while porosity is fixed at 0.2. The caprock is initially water-wet and altered by exposure to dissolved CO 2 , such that the wettability altering agent is X WA = X CO 2 w for all 1D-V simulations.
The saturation functions in the aquifer are set equal to the final state and are kept static. We employ the same end-wetting state parameters used in the horizontal study (refer to Table 1) for the aquifer. On the other hand, different initial-wet entry pressure and permeability values are tested for the caprock section (see Table 4).
The initial condition, depicted in Figure 7a, is set such that a 50% CO 2 saturation is uniformly distributed in a column of height h = 70 m. The reservoir pressure is considered to be hydrostatic with regard to brine density and depth of the formation. For this case, the temperature and salinity of the reservoir are considered to be 20°C and 0 ppm, respectively. The diffusion coefficient of CO 2 in the brine wetting phase is set to 2 × 10 −9 m 2 /s. All boundaries are closed to flow, and the residual saturation of each phase is S rα = 0.2. CO 2 mole fraction in the brine phase is set initially to 5 × 10 −3 .
The above described condition is initially not at equilibrium due to capillarity and gravity gradients. In the absence of any other driving forces, the CO 2 and brine will redistribute in the aquifer according to buoyancy and capillarity to reach an equilibrium. Simultaneously, CO 2 dissolves into the brine and diffuses into the caprock, altering the capillary entry pressure over time. Since the capillary pressure in the aquifer is small, the equilibration in the aquifer with regard to gravity is in order of days, while the significant WA effects in the caprock starts in the order of years.

Final-wet-pc
Final-wet-kr  It is clear that the capillary seal is sufficient to contain CO 2 under the caprock, and containment is not sensitive to the parameters of the relative permeability functions. We also observe the resulting redistribution of CO 2 in the aquifer from the initial condition in Figure 7a results in a near complete gravity segregation of CO 2 and brine, leaving residual CO 2 (with saturation of 0.2) below the accumulated CO 2 column.
In comparison, when the end-wetting capillary pressure curve is directly applied in the caprock (Figures 7d-e), the capillary seal is no longer sufficient to contain CO 2 . The difference between Figures 7d and 7e shows the extent to which the wetting properties of the relative permeability affects CO 2 migration in the caprock over time. This result shows if the wettability state of an initially water-wet caprock is changed by exposure to dissolved CO 2 , then CO 2 can eventually migrate into the caprock. The dynamics of the wettability change coupled with diffusion of dissolved CO 2 upward into the caprock will ultimately determine the behavior of CO 2 in this system over time.
Given the dimensions and initial conditions of this closed system, then the maximum amount of CO 2 migrating into the caprock can be estimated. Neglecting the dissolved CO 2 in the brine and assuming the non-wetting phase is mostly CO 2 , then the initial mass of CO 2 in the system is M i CO2 /A ≈ hρ n S i n φ = 5 × 10 3 kg/m 2 . Considering the residual CO 2 in the aquifer, this results in a maximum migration of CO 2 into the caprock of ca. M max CO2 /A ≈ 3 × 10 3 kg/m 2 . Figure 8 shows the accumulated CO 2 in the caprock over time for the combinations of final capillary pressure function with the initial-and final-relative permeability functions for three different caprock permeability values. As expected, we observe that the CO 2 migrates faster into the caprock for the final-wetting functions, the lower the rock permeability the slower the CO 2 migrates into the caprock, and the amount of CO 2 is limited by the initial and residual saturation (M max CO2 /A = 2733 kg/m 2 ). We now introduce dynamics into the wetting state of the caprock. In this vertical system, the WA agent is modeled as dissolved mass fraction, i.e., X WA = X CO 2 w . We recall that the alteration process, or the measure of exposure time, depends on the magnitude of the WA agent. The dissolved CO 2 mass fraction in water is of order 10 −2 while the non-wetting saturation is of order 10 −1 . To observe the WA effects within years in the 1D-V studies, we choose relatively higher values for model parameters b 1 , b 2 , ν 1 , and ν 2 (Table 4) in comparison to the ones we considered in the 1D-H simulations (Table 2).
Similar to the studies in Section 3.2, we perform a sensitivity analysis to quantify the impact of dynamic saturation functions on CO 2 migration into the caprock by varying rock properties (permeability and initial entry pressure) and dynamic parameters (β(C) and η(C). Here, we examine the impact of WA on the integrity of caprock by • Case 1: considering WA dynamics only in the capillary pressure function and • Case 2: considering WA dynamics in both the capillary pressure and relative permeability functions. Table 4: Parameters describing the relationship between the pore-scale parameter C associated with contact angle change and dynamic coefficients β and η for the 1Dvfs (t ch = 10 7 s, [−] 10 −5 10 −2 4.5×10 −1 10 7 1.8 4.999 × 10 3 5 × 10 −1 2.5 × 10 −5 5.2×10 −2 3.75×10 −1 5×10 −5 1.81×10 −1 2.5×10 −1 Table 5 presents the parameters combinations for the quantification of the two cases mentioned above. Table 5: Tested parameter values used in the 1D vertical simulation study. A set of simulations were performed for each individual value of c i and K in which C is varied, where the isolated (only capillary pressure) and combined impact on dynamic wettability functions is studied. For each value of C, the corresponding value of β and η is given in Table 4.
1 × 10 −16 10 −5 5 × 10 4 7.5 × 10 −17 2.5 × 10 −5 10 5 5 × 10 −17 5 × 10 −5 An example of the impact of dynamics in wettability on the vertical CO 2 distribution is shown in Figure 9 for the parameter combination of c i =10 4 Pa, K = 10 −16 m 2 , and C = 5 × 10 −5 under two different dynamic cases. For Case 1 (Figure 9a), the caprock capillary entry pressure is reduced by exposure to dissolved CO 2 , allowing mobile CO 2 to migrate upwards into the low permeability domain. For Case 2 (Figure 9b), an added effect occurs in the caprock. CO 2 migrates upwards into the caprock at a higher saturation when relative permeability is altered, which is in agreement with the shift in fractional flow function for vertical flow from initial to final wetting condition ( Figure 2d). Figure 9c shows the evolution of WA induced CO 2 saturation over time in the first grid cell within the caprock domain. For this set of parameters, we observe breakthrough into the caprock after approximately 10 years. Once CO 2 begins to migrate vertically, both cases are characterized by a steep increase followed by a more gradual accumulation of CO 2 in the selected grid cell.
A set of simulations were performed for each individual value of c i and K in which C is varied (Table 5). These numerical results are shown in Figure 10a for Case 1 (isolated P c dynamics) and Figure 10b for Case 2 (combined dynamics). Comparing the two cases, we observe the build-up of CO 2 follows qualitatively the same response to changes in parameter values, but where the amount of CO 2 in the caprock is less when dynamics are limited to capillary pressure function (Case 1, Figure 10a) than when dynamics are presented in both saturation functions (Case 2, Figure 10b).
As expected, the larger the value of caprock entry pressure, the longer it takes for the CO 2 to start entering the caprock. This starting migration time increases with smaller values of C (slower dynamics on the saturation functions). On the other hand, smaller values of rock permeability result in slower CO 2 migration into the caprock (see also Figure 8), which is also affected by the value of the pore-scale parameter C and entry pressure c i as observed in Figure 10.
Given the above observations, it may be possible to obtain a general relationship for the impact of WA on CO 2 containment. We first divide the variables by reference values to make them dimensionless. The reference curves for the scaling are the ones giving by the static saturation functions (purple curves in Figures 8 and 10). We have tested different functions for the scaling and selected the ones we describe next. For the time variable, we suggest  a translation as a function of the entry pressure and pore-scale parameter: where a 0 and a 1 are fitting parameters. For the CO 2 mass in the caprock, we suggest the following scaling function: where a 2 ,. . . , a 6 are fitting parameters.
From Figure 8 we observe a linear behavior between CO 2 mass and time, until reaching a maximum value of mass. This motivates the following relationship between scaled mass and time: where ψ 0 is a fitting parameter. Table 6 shows the values of the reference and fitting parameters for both cases. Figure 11 shows the different scaled simulation results. Thus, we have obtained simple models to describe the impact of WA change on CO 2 containment given caprock parameters K and c i , and dynamic WA parameter C.

Discussion
We have implemented dynamic capillary pressure and relative permeability functions to simulate and quantify the impact of WA on CO 2 storage and containment. Horizontal and vertical field-scale test cases were used to demonstrate the effect of WA in CO 2 front migration. These simplified scenarios allowed for testing of a wide range of different rock-fluid and wettability parameters in order to understand the impact of WA in relation to other multiphase flow processes. Our findings show that long-term WA can be characterized according to capillary number for horizontal flow and to caprock integrity parameters for vertical migration. These results can be used to understand the response of CO 2 storage mechanisms to wettability alteration by long-term exposure to supercritical and dissolved CO 2 .

CO 2 storage efficiency
The horizontal case study that focused on the impact of WA on storage efficiency shows the importance of dynamics in both saturation functions, P c and k rα , during CO 2 injection. Wettability change sharpens and delays the front movement, and thus a greater CO 2 storage efficiency results with increasing exposure to injected CO 2 . A defining feature of long-term WA is time-dependent heterogeneous wettability, resulting in a wettability gradient outwards from the injection well. Longer exposure time towards the inlet significantly alters wettability over time, while the leading edge of the CO 2 front remains near water-wet. The interplay between dynamically changing saturation functions and heterogeneous wettability leads to an increasingly complex evolution of CO 2 migration with increased injected volumes, even for a simple 1D system. Altered wetting conditions can spontaneously draw CO 2 back towards the inlet by capillary action. These interesting dynamics that may occur over months and years and can have important implications for storage efficiency in systems where WA is expected. More investigation is needed to fully understand the complex behavior caused by long-term WA in non-idealized storage systems.
The horizontal case study results also provide additional insight to the role of WA on storage efficiency as a function of capillary number. The advancing CO 2 front, which can be considered a proxy for storage efficiency, can be significantly impacted at lower capillary numbers, but where the rate of dynamics in the saturation functions plays a important role. At higher capillary numbers, storage efficiency is much less impacted by long-term WA, regardless of the underlying parameters controlling the wettability dynamics. This result gives important insight into the flow regimes where long-term WA is a relevant process (and similarly, where WA dynamics could be neglected). CO 2 storage systems often have a range of flow regimes across the domain at any given time, with highly viscous flows near the injection well and capillary-dominated flow further afield.
Another important finding is that the observed impact of wettability dynamics at the field-scale are controlled mainly by dynamics in the capillary pressure function and are mostly insensitive to the dynamics in relative permeability. We recall that the relative permeability dynamics we modeled represented a quite severe transition in curvature, and therefore the expected dynamics if using lab-derived functions will be even smaller. Therefore, it is reasonable to conclude that one can neglect dynamics in relative permeability and only focus on dynamics in capillary pressure for simulating field-scale impacts of long-term changes in wetting condition due to CO 2 exposure.
We emphasize that this conclusion only applies to dynamics in the curvature of the relative permeability functions, and the change in end-point residual saturation due to CO 2 exposure is still a relevant aspect to incorporate at the field-scale. The possibility for including this process in future work is discussed more in Section 4.4.

CO 2 containment
The second vertical case with CO 2 migration into the overburden due to long-term WA shows how the same dynamics in P c and k r impact storage containment. In contrast to the horizontal case, the vertical WA process is dependent on a certain sequence of effects: (1) diffusion of dissolved CO 2 into the caprock; (2) reduction of caprock entry pressure by CO 2 exposure that allows supercritical CO 2 to break the capillary seal; (3) vertical migration of CO 2 in caprock. We note that the impact of altered relative permeability only acts once CO 2 is mobile in the caprock. This sequence of processes implies that even if WA induces a loss of containment, it is first and foremost a slow process. CO 2 continually encounters unexposed caprock along its vertical migration path that must be altered by the same slow process in order for vertical flow to continue.
Although we observe that long-term WA by CO 2 exposure leads to loss of containment, our simulations show that CO 2 advances very slowly, but steadily, over time. Our scaling analysis shows that CO 2 migration follows a characteristic evolution that is amenable to scaling by the underlying parameters, which gives insight into a unified model for long-term CO 2 migration into a caprock due to long-term WA. The scaling model is a function of the rate of WA dynamics, caprock permeability, and initial capillary entry pressure. The model is an important generalization and valuable for making use of experimental/surveyed data to predict the integrity of caprocks exposed to CO 2 over long timescales, i.e., one can quantify a priori the potential for unwanted CO 2 migration without having to perform field-scale simulations. More importantly, this scaling relationship shows that long-term WA poses very little risk to CO 2 containment. Only with a dramatic reduction in the initial strength of the capillary seal and relatively high caprock permeability will CO 2 migrate vertically in non-negligible quantities over several decades. For the Sleipner CO 2 storage project where caprock permeability is on the order of 100 nanoDarcy and a plume that covers an area of 10 km 2 , this would translate to ca. 5 tons over 200 years (assuming WA on both saturation functions, c i =3 MPa, and C = 10 −5 ).

Numerical implementation of dynamic flow functions
All numerical studies presented in this paper are conducted using the open-source framework OPM. One of the advantages of using OPM is that utilizes modern hardware architecture (i.e., multiple cores and vectorization). We extended the TCTP and flow simulators in OPM to include dynamics in the saturation functions (P c and k rα ).
Studies on the 1D-H system for large Ca by increasing the injection flow rate would require a large domain to avoid the saturation front to reach the boundary, while keeping the grid size small to reduce numerical errors in the computation of the scaled front location difference (SFLD). Smaller values of K also result in larger Ca while keeping the saturation front closer to the injection well. However, using small grids close to the injection well result in convergence issues. This is expected since the general purpose of OPM Flow is for large field-scale simulations.
Then, one could modify further the code (specially the well models) to conduct studies with smaller grids. The time scale for the numerical studies on the 1D-V system is of order of decades, which is computationally expensive for simulations using full TPTC models since dynamics in the saturation functions in addition to computations of phase compositions restrict the size of the time step. Further investigation on solution strategies for this system is required to reduce the computational cost (e.g., by not updatingχ after each time step but at certain times).

Extension to realistic systems
This study has applied the dynamic WA model to relatively simple 1D systems under controlled conditions and homogeneous media. Additional effects become relevant if we are to extend this understanding to realistic 2D and 3D storage reservoirs. For example, CO 2 migration in the reservoir is more complex when affected by gravity and reservoir heterogeneity. In the horizontal 1D systems, we observed a longer CO 2 exposure at the inlet versus further out in the reservoir. With the introduction of gravity, this dynamic could lead to slower upward migration for CO 2 injected down-dip or at the bottom of a thick reservoir. Heterogeneity either in permeability, porosity, or initial wettability will add complexity to dynamics in wettability that needs to be further studied. We note that differences expected in realistic systems apply mainly for migration within the reservoir itself. Vertical migration in a low-permeability caprock is primarily a 1D process even in 3D systems. The main difference in moving to 2D/3D is the possibility for varying CO 2 column height under the caprock interface, which plays a role in the buoyancy forces acting on the rising CO 2 .
Another important aspect is the impact of wetting dynamics on CO 2 trapping efficiency as a residual phase.
We have not accounted for dynamics in residual saturations in our simulations as the examples used in this study are mainly focused on drainage processes. However, real systems will experience imbibition, especially for CO 2 storage dependent on migration-assisted trapping. Intermediate-wetting conditions will reduce the capillary trapping capacity of reservoir rocks, which leads to a greater portion of CO 2 remaining in a free phase. The change in residual trapping from water-to CO 2 -wetting conditions by CO 2 exposure will be strongly coupled to the impacts observed due to dynamic alteration of the flow functions studied here. Although residual trapping has been characterized under static wettability, more work is needed to develop models for dynamic changes in residual trapping due to CO 2 exposure and couple them to the dynamic flow functions implemented in this study. Future studies may address this by adapting the approach proposed by Kassa et al. [2020Kassa et al. [ , 2021b; however, with additional complexity to the porous medium (i.e., pore-network model instead of parallel capillary tubes).

Additional pore-scale impacts
In this study, we have seen a marked difference in horizontal and vertical CO 2 migration depending on the rate of wettability alteration as controlled by the pore-scale parameter C. We recall that the implemented dynamic flow functions have been upscaled from pore to core scale. Our results show how pore-scale dynamics are manifested at the field scale and the importance of quantifying C through laboratory experiments. To date, only a few experiments have measured contact angle change through exposure to CO 2 over long timescales, but none have described the temporal evolution of contact angle needed to calibrate pore-scale CA change models such as Equation (9). More experimental evidence is needed to determine if continual CO 2 exposure could fully transform a water-wet rock to strongly CO 2 -wet over the long timescales relevant for CO 2 storage.
There are other pore-scale effects that can be of interest for further study, particularly with regard to geochemical impacts. First, CO 2 -induced reactions can dissolve and precipitate minerals, thus altering the pore topology of reservoir rocks and caprocks [Landa-Marbán et al., 2021]. Changes in pore topology could alter the flow functions in a similar manner as changes in contact angle, but characterization of the connection between topological alteration and capillary pressure/relative permeability curves requires further investigation. Geochemical alteration can also lead to an increase in permeability, which could enhance vertical CO 2 migration in addition to WA. On the other hand, an increased porosity can counteract the permeability increase. Further study is needed to quantify and characterize the combined impact of geochemical and wettability alteration on CO 2 migration.

Conclusion
In this work, we present a macroscale TPTC mathematical model to study long-term WA effects in CO 2 storage applications. Particularly, this model includes macroscale dynamic capillary pressure and relative permeability functions derived from pore-scale WA models. We use OPM to implement the model and perform the numerical studies. Two simple 1-D systems are considered to investigate the WA effects on different flow regimes.
Based on this work our conclusions are as follows: • Horizontal and vertical field-scale test cases can be used to demonstrate the effect of long-term WA in CO 2 front migration.
• Long-term WA can be characterized according to capillary number for horizontal flow and to caprock integrity parameters for vertical migration.
• Field-scale impacts of long-term WA are mostly controlled by dynamics in the capillary pressure function due to CO 2 exposure, while similar dynamics in the curvature of the relative permeability functions is a secondary factor.
• Scaling models to quantify CO 2 migration into the caprock show that long-term WA poses little risk to CO 2 containment. c i , c f Entry pressures (initial and final) θ i , θ f Contact angle (initial and final)

Notation
x i , x dy Non-wetting phase front location (initial and dynamical saturation functions)