# Impact of Time-Dependent Wettability Alteration on the Dynamics of Capillary
Pressure

Abay M. Kassa1,2 Sarah E. Gasda2 K. Kumar1 F. A. Radu1

###### Abstract

Wettability is a pore-scale property that has an important impact on
capillarity, residual trapping, and hysteresis in porous media systems. In
many applications, the wettability of the rock surface is assumed to be
constant in time and uniform in space. However, many fluids are capable of
altering the wettability of rock surfaces permanently and dynamically in time.
Experiments have shown wettability alteration can significantly decrease
capillarity in CO${}_{\textrm{2}}$ storage applications. For these systems,
the standard capillary-pressure model that assumes static wettability is
insufficient to describe the physics. In this paper, we develop a new dynamic
capillary-pressure model that takes into account changes in wettability at the
pore-level by adding a dynamic term to the standard capillary pressure
function. We simulate the dynamic system using a bundle-of-tubes (BoT)
approach, where a mechanistic model for time-dependent contact angle change is
introduced at the pore scale. The resulting capillary pressure curves are then
used to quantify the dynamic component of the capillary pressure function.
This study shows the importance of time-dependent wettability for determining
capillary pressure over timescales of months to years. The impact of
wettability has implications for experimental methodology as well as
macroscale simulation of wettability-altering fluids.

1 Department of Mathematics, University of Bergen, P. O. Box 7800, 5020
Bergen, Norway.
2 NORCE, Nygårdsgaten 112, 5008 Bergen, Norway.
Corresponding author: Abay M. Kassa (E-mail<EMAIL_ADDRESS>
##  1 Introduction

Wettability plays an important role in many industrial applications, in
particular subsurface porous media applications such as enhanced oil recovery
(EOR) and CO${}_{\textrm{2}}$ storage [1, 2, 3, 4, 5]. The wetting property of
a given multiphase system in porous media is defined by the contact angle,
which describes the affinity of one fluid for the surface of the medium over
another. The contact angle is determined from surface chemistry and associated
forces acting at the molecular scale along the fluid-fluid-solid interface
[1]. In porous media applications, microscale wettability determines the
strength of pore-scale capillary forces and the movement of interfaces between
individual pores. At the macroscale, wettability impacts upscaled quantities
and constitutive functions such as residual saturation, relative permeability,
and capillary pressure, which in turn affect large-scale multiphase flow
behavior.

The standard assumption is that wettability is a static property of the
multiphase system. However, the composition of many fluids is capable of
provoking the surfaces within pores to undergo wettability alteration (WA) via
a change in contact angle. A change in contact angle can alter capillary
forces at the pore scale, and thus affect residual saturations of the system
[6, 7]. This effect has been exploited extensively in the petroleum industry,
where optimal wetting conditions in the reservoir are obtained through a
variety of means that includes chemical treatment, foams, surfactants and low-
salinity water flooding (see for example [8, 9, 10, 11, 12]).

Wettability is also recognized as a critical factor in geological
CO${}_{\textrm{2}}$ sequestration which exerts an important role on caprock
performance [13, 14]. WA can influence the ability to prevent
CO${}_{\textrm{2}}$ leakage from the reservoir by altering the capillary
forces that act to trap CO${}_{\textrm{2}}$ as buoyant phase beneath a low-
permeability caprock. This process is highly dependent on CO${}_{\textrm{2}}$
being a strongly non-wetting fluid, and WA may lead to conditions that allow
for CO${}_{\textrm{2}}$ leakage [15, 16]. In addition, WA can affect residual
saturation, results in an impact on trapping efficiency of injected
CO${}_{\textrm{2}}$ [4]. Therefore, reliable quantification of wettability is
key to understand the macroscale processes in subsurface applications.

Despite the fact that WA is known to impact macroscale capillarity and
relative permeability behavior, few detailed measurements are available to
characterize the alteration of the constitutive function themselves. Plug et
al., [17] have reported brine-CO${}_{\textrm{2}}$ (gas, liquid) drainage-
imbibition experiment and showed capillary instability for supercritical
CO${}_{\textrm{2}}$-brine system, meaning that the capillary pressure
measurements changes steadily over time. The imbibition curve also exhibited a
significant deviation from the expected curve that was not explained by
classical scaling arguments. The authors proposed WA as an explanation but did
not explicitly measure any changes in contact angle. Additionally, recent
experiments measured capillary pressure curves for a silicate sample using a
fluid pairing of supercritical CO${}_{\textrm{2}}$ and brine [18]. Repeated
drainage-imbibition cycles were performed over a period of 6 months, and a
clear reduction in capillary pressure was recorded for each subsequent
drainage cycle. The authors also attributed these deviations from the expected
capillary curve to a change in wettability of the rock sample over time due to
CO${}_{\textrm{2}}$ exposure, similar to aging. This hypothesis was confirmed
through observations of a wetting angle increase from $0^{\circ}$ to
$75^{\circ}$ after 6-months exposure. It is also reported similar
$P_{\textrm{c}}$-$S$ instability and deviation in dolomite/carbonate [19], and
quartz [20, 19] sands for scCO${}_{\textrm{2}}$-brine system. More literature
on WA and $P_{\textrm{c}}$-$S$ measurements can be found in [14].

The above experiments reveal that capillary pressure curves are not static
when significant WA occurs, despite the fact they were performed following the
standard multi-step outflow procedure, i.e. where “equilibrium” is obtained
after each incremental step in pressure. Therefore, the standard equilibrium
capillary pressure models cannot be readily applied without additional
dynamics to capture the impact of WA. Capillary pressure dynamics due to WA is
fundamentally different from the well known non-equilibrium flow dynamics (for
example see [21, 22]) and contact angle hysteresis, i.e. receding and
advancing contact angle [23, 24] because these dynamics disappear when the
system is at equilibrium.

In petroleum reservoir simulation, there have been some theoretical
investigations that assume an instantaneous WA to model the constitutive
relations and immiscible flow simulations [25, 26, 3, 27, 28]. The WA
mechanism is incorporated by interpolation between the final and initial
wetting-state constitutive relations [27, 25]. Lashgari and co-authors [26]
have derived an instantaneous WA model from Gibbs and Langmuir adsorption
models. The proposed WA model is coupled to $P_{\textrm{c}}$-$S$ and relative
permeability curves through residual saturation. These models neglect exposure
time dependency of WA and also the impact of pore-scale heterogeneity in
wetting properties in the macroscale simulations. In the available literature,
only one study [29] has included the effect of exposure time on WA and
constitutive relations for core scale simulation. But this numerical study
does not sufficiently incorporate or upscale pore-scale processes to
macroscale constitutive laws.

To our knowledge, a rigorous mathematical characterization of dynamics in
$P_{\textrm{c}}$-$S$ functions introduced by exposure to a WA agent has not
been previously performed. In this paper, we develop a new model for dynamic
capillary pressure that is derived using direct simulation of
$P_{\textrm{c}}$-$S$ curves from a pore-scale model. WA is introduced at the
pore scale using a mechanistic model for contact angle change as a function of
exposure time to a WA fluid. This mechanistic model is not intended as a
replacement for detailed surface chemistry models, but is used merely as a
convenient mathematical way to represent changing contact angle in time. We
acknowledge that the field of wettability alteration at the micro- and nano-
scales is an active area of research covering a broad range of complexity. We
emphasize that this study does not directly address wettability alteration at
pore surfaces, which is well beyond the scope of this paper.

The pore-scale flow model is represented by a cylindrical bundle-of-tubes
model. This pore-scale approach is simplistic in that it does not represent
interactions between pores or realistic porous geometry (i.e. converging-
diverging throat diameters). Additionally, there is no natural hysteresis or
residual saturation with cylindrical tubes. However, a simple pore-scale model
is sufficient as a first step towards develop time-dependent correlated
expressions for $P_{\textrm{c}}$-$S$. It is also important to investigate the
impact of wettability alteration apart from other complicating processes that
can occur in more realistic pore-network models.

##  2 Approach

The extended $P_{c}$-$S_{w}$ relationship introduces a dynamic component that
captures the changing wettability as measured by the deviation of the dynamic
capillary pressure from the equilibrium (static) capillary pressure. This
relationship can be described as follows:

$P_{c}(\cdot)-P_{c}^{\rm st,i}:=f^{\rm dyn}(\cdot),$ (1)

where $P_{c}^{\rm st,i}$ represents the capillary pressure for the system
given a static initial wetting state, while $f^{\rm dyn}$ represents the
deviation from the static state. The static curve can be characterized by a
number of well-known capillary pressure models, such as van Genuchten or
Brooks-Corey. Static systems studied here are described by the Brooks-Corey
model,

$P_{c}^{\rm st}={c_{w}}\left(\frac{S_{w}-S_{wc}}{1-S_{wc}}\right)^{-a_{w}},$
(2)

where $c_{w}$ is the entry pressure, $1/a_{w}$ is the pore-size distribution
index, whereas $S_{wc}$ is the residual water saturation. The static curve is
equivalent to the equilibrium capillary pressure for systems that do not
undergo wettability alteration.

The objective of this study is to characterize and quantify the dynamic term
$f^{\rm dyn}$, the key term of interest in the $P_{\textrm{c}}$-$S$ model, for
a system that undergoes wettability alteration. We propose an interpolation
model can be applied for this system, where the dynamic term acts to shift the
capillary pressure between two end wetting states in time. This type of model
has been used in petroleum reservoir simulation [27, 25, 3, 26]. To obtain an
interpolation model, the dynamic component in Equation (1) can be scaled by
the difference between two static curves, each representing the initial and
final wetting-state capillary pressure curves, to give a non-dimensional
quantity we call the _dynamic coefficient_ , which is defined as follows,

$\omega\big{(}P_{c}^{\rm st,f}-P_{c}^{\rm st,i}\big{)}=f^{\rm dyn},$ (3)

where $P_{c}^{\rm st,f}$ is the final wetting state capillary pressure. This
can be substituted into Equation (1) to obtain a dynamic interpolation model
for capillary pressure,

$P_{c}=(1-\omega)P_{c}^{\rm st,i}+\omega P_{c}^{\rm st,f},$ (4)

We note that in the model presented above, we have defined the “total”
capillary pressure $P_{c}$ as simply the measured difference in phase
pressures at any point in time. In a reservoir simulation, this would be
capillary pressure in a given grid cell, whereas in a multi-step outflow
experiment, it corresponds to the pressure drop across the experimental cell.
For quasi-static displacement in a bundle of capillary tubes, $P_{c}$ is the
difference in reservoir pressures, i.e., $P_{c}=P_{l}^{res}-P_{r}^{res}$, (see
Figure 1). The exact nature of $\omega$ and its functional dependencies can
only be determined from a full characterization of $P_{\textrm{c}}$-$S$ curves
under different conditions. These curves can be derived from laboratory
experiments, but this approach is costly and time-consuming. Alternatively,
one may take a more theoretical approach by simulating $P_{\textrm{c}}$-$S$
curves using a pore-scale model that includes the impact of WA.

For a system that undergoes WA, a significant change in contact angle could
lead to the wetting phase becoming non-wetting and vice versa. For clarity, we
will continue to use $w$ subscript for the phase that was originally wetting
and $nw$ for the phase that is originally non-wetting regardless of the actual
state of wettability in the system.

###  2.1 Capillary flow in cylindrical tubes

There are various choices of pore-scale models available. The easiest to
implement and analyze is the bundle of tubes (BoT) model which is a collection
of capillary tubes with a distribution of radii. Herein, we describe the model
and the approach for implementation of time-dependent WA within individual
capillary tube.

The BoT model has been applied in many immiscible multiphase flow studies and
is employed to study constitutive relations and fluid dynamics in general [30,
22, 31]. The BoT model is a popular approach due to the simplicity of
implementation and the ability to study the balance of energy and forces
directly at the scale of interfaces. The average behaviour of a BoT model can
then be used to inform better constitutive functions at the macroscale.
However, direct use of simulated capillary or relative permeabilities curves
is usually not advised. This is because a collection of capillary tubes is not
representative of the interacting nature of pores and pore throats. The
geometry is over-simplified with respect to converging and diverging pore
throats. Some of the shortcomings of BoT models can be addressed through
interactive BoT models [32], or tubes of more complex geometry. These more
advanced BoT models are beyond the scope of this paper.

In this paper, we consider cylindrical bundle of tubes having length, L, and
is shown in Figure 1.

Figure 1: Fluid displacement in non-interactive bundle of tubes (BoT). Here,
the left reservoir contains non-wetting fluid that displaces the wetting fluid
to the right and vice versa.

These tubes are designed to connect wetting (right) and non-wetting (left)
reservoirs with pressures labeled as $P_{\rm r}^{\rm res.}$ and $P_{l}^{\rm
res.}$ respectively. Let the reservoir pressures difference be defined as:

$\Delta P=P^{res}_{l}-P^{res}_{r}.$ (5)

Equation (5) indicates that the pressure drop across the bundle is equal to
the measured capillary pressure $P_{c}$ defined previously in Equation (1).

We consider the tubes to be filled with water and are water-wet initially. To
displace the water phase, the pressure drop has to exceed the entry capillary
pressure [22]:

$\Delta P>P_{c}^{\rm mic}(R_{m},\theta_{m}),$ (6)

where $P_{c}^{\rm mic}(R_{m},\theta_{m})$ is given by the Young’s equation

$P_{c}^{\rm
mic}(R_{m},\theta_{m})=\frac{2\sigma_{ow}\cos(\theta_{m})}{R_{m}},~{}m=1,2,\cdots,N,$
(7)

where $R_{m}$ is the $m$-th tube radius, $N$ stands for number of tubes,
$\sigma$ is fluid-fluid interfacial tension (IFT) and $\theta_{m}$ measures
the wetting property of the rock surface (it may be static or dynamic in time
per the $m$-th tube). As long as condition (6) is satisfied, the fluid
movement across the length of a given tube can be approximated by the Washburn
flow model,

$q_{m}=\frac{R_{m}^{2}(\Delta P-P_{c}^{\rm mic}(R_{m},\theta_{m}))}{8(\mu_{\rm
nw}x_{m}^{i}+\mu_{\rm w}(L-x_{m}^{i}))},$ (8)

where $q_{m}$ is the interface velocity, and $\mu_{\rm nw}$ and $\mu_{\rm w}$
are non-wetting and wetting fluid viscosities, respectively. The above
expression for interface velocity is equivalent to the rate of change of the
interface location $x^{i}_{m}$ in tube $m$, which can be defined as
${dx_{m}^{i}}/{dt}$. From equation (8), one can then determine the required
time to reach a specified interface position, and subsequently the time to
drain each individual tube in a multi-step outflow experiment.

###  2.2 Pore-scale time-dependent wettability model

In this study, the WA agent is defined as either the non-wetting fluid itself
or some reactive component therein that is transferred to the wetting fluid
across a fluid interface. We consider an alteration process that continues
until the ultimate wetting state is reached. The alteration is permanent, but
can also be halted at some intermediate wettability state if the WA agent is
removed from the system before the final wetting state is reached. If the
agent is reintroduced at some later point, alteration continues until the
final state.

To this end, we introduce a general functional form of pore-scale WA
mechanism,

$\theta_{m}(\cdot):=\theta_{m,\rm i}+\varphi(\cdot){\rm\Delta}{\rm\Theta},$
(9)

where ${\rm\Delta\Theta}=\theta_{m,\rm f}-\theta_{m,\rm i}$, $\theta_{m,\rm
f}$ and $\theta_{m,\rm i}$ are the ultimate and initial contact angles
respectively. Note that $\varphi$ can be any function form that is properly
designed to capture the WA processes over a period of time. As discussed in
the introduction, the process of contact angle change at the pore-scale
involves surface chemistry acting at micro and nano-scale. The mechanistic
model in Equation (9) is not intended to capture the full complexity at the
pore-scale, but is simply a convenient mathematical form that gives smooth
change in contact angle over time between an initial and final contact angle.
Detailed laboratory work would needed to further justify the use of this
model, which is beyond the scope of this paper.

Though it is well known that WA impacts the reservoir system, theoretical
investigation and laboratory measurement on time-dependent WA is very limited.
However, coupling the instantaneous WA processes to the Langmuir adsorption
model have been successfully implemented in EOR simulation technology [27, 25,
3, 26]. In this paper, the parameter $\varphi(\cdot)$ in equation (9) is
assumed to evolve according to Langmuir-type adsorption as a test case. Thus,
$\varphi$ can be defined as follows:

$\varphi:=\frac{\chi_{m}}{C+\chi_{m}}$ (10)

where $C$ is non-dimensional parameter that controls the speed and extent of
alteration from water-wet to intermediate-wet or hydrophobic system. The non-
linear variable, $\chi_{m}$ is a measure of the local exposure time of each
individual tube $m$, defined as

$\chi_{m}(t,S_{nw_{m}}):=\frac{1}{T}\int_{0}^{t}S_{nw_{m}}d\tau,$ (11)

where $T$ is a pre-specified characteristic time. In this paper, the total
time to complete the possible drainage-imbibition cycles from the initial to
final wetting state is used as a characteristic time.

For convenience, we may also define $\overline{\chi}$ as a measure of exposure
time across the entire bundle (i.e. REV),

$\overline{\chi}(t,S_{nw}):=\frac{1}{T}\int_{0}^{t}S_{nw}d\tau.$ (12)

Given the above model in Equation (10), we may consider two WA approaches,
namely the uniform and non-uniform alteration models.

In the case of uniform WA, pores are assumed to be altered identically in time
due to transfer of the WA agent by dissolution from neighboring pores. In this
case, the local exposure time is equivalent to the average exposure time, such
that $\chi_{m}=\overline{\chi}$ in Equation (10). Alternatively, we may
restrict WA to _only_ those pores that have been flooded by the non-wetting
fluid. In this case, we must keep track of the local exposure time of each
tube within the bundle to evaluate Equation (10). This process gives rise to a
_non-uniform_ or mixed alteration of the pores until such time that all pores
have reached the final wetting state.

###  2.3 Simulation approach

The uniform and non-uniform approaches for contact angle evolution above are
coupled into the BoT model following the algorithm pictured in Figure 2.

Figure 2: Algorithm for wettability alteration and capillary pressure
simulation following drainage-imbibition displacements in a pore-scale model
represented by a bundle of tubes.

For a drainage cycle, the pressure drop is increased incrementally. Upon each
increase, the respective tubes are drained according to the entry pressure
inequality. The Washburn model, Equation (8), is used to calculate the elapsed
time to drain each tube. Afterwards the contact angle (and thus $P_{c}^{\rm
mic}$) is updated according to Equation (9). For the uniform model, the
contact angle change is applied uniformly, while in the non-uniform model, the
change is applied according to local exposure time. The average saturation is
recorded and the cycle is repeated until all tubes have drained.

During imbibition, the process is reversed. For the non-uniform model, the
contact angle continues to change locally until the non-wetting phase is
displaced by the wetting phase. For the uniform model, the alteration occurs
continuously throughout imbibition according to the average exposure time
$\overline{\chi}$.

Repeated drainage and imbibition cycles are run according to this algorithm
until the final wetting state is reached. At the completion of the numerical
experiment, we obtain a set of $P_{\textrm{c}}$-$S$ “data points” that can be
plotted in the usual way.

###  2.4 Correlation approach

Once the capillary pressure and associated dynamics have been simulated, the
resulting curves used to quantify the dynamic coefficient in the interpolation
function in Equation (4). The goal is develop a correlated model that involves
only a single parameter, and this parameter should have a clear relationship
to the pore-scale WA process such that changes in the pore-scale WA model
parameter $C$ can be directly accounted for in the dynamic coefficient.

We note that for the interpolation model, the dynamics are contained in the
coefficient $\omega$ that interpolates between two static capillary pressure
curves. Therefore before correlating this model, we determine the functional
form of the static wetting states using a simple power-law model.

##  3 Results

In this section, we present the results of simulated capillary pressure
experiments using pore-scale models represented by a bundle of tubes following
the approach described above. First, the static capillary pressure is
simulated for both the initial and final wetting states. Then, simulations are
performed on the model subject to WA. Two cases are performed: uniform and
non-uniform alteration. The correlated $P_{\textrm{c}}$-$S$ models for the
static cases are presented, followed by the correlated models for the dynamic
capillary pressure curves.

###  3.1 Bundle of tubes model set-up

The pore scale is described by a BoT model (see Section 2.1). Each tube in the
BoT is assigned a different radius R, with the radii drawn from a truncated
two-parameter Weibull distribution [31]

$f(R)=\frac{\Big{[}\frac{R-R_{\rm min}}{R_{\rm
av}}\Big{]}^{\eta-1}\frac{\eta}{R_{\rm av}}{\rm
exp}\Big{(}-\Big{[}\frac{R-R_{\rm min}}{R_{\rm
av}}\Big{]}^{\eta}\Big{)}}{1-{\rm exp}\Big{(}-\Big{[}\frac{R_{\rm max}-R_{\rm
min}}{R_{\rm av}}\Big{]}^{\eta}\Big{)}}$ (13)

where $R_{\rm max}$, $R_{\rm min}$, and $R_{\rm av}$ are the pore radii of the
largest, smallest, and average pore sizes, respectively, and $\eta$ is a
dimensionless parameter. The rock parameters and fluid properties are listed
in Table 1.

parameters | values | unit | parameters | values | unit
---|---|---|---|---|---
$\sigma_{ow}$ | 0.0072 | N/m | no. radii | 500 | [-]
$R_{\rm min}$ | 10 |  $\mu$m | $R_{\rm max}$ | 40 |  $\mu$m
$\theta_{f}$ | 80 | degree | $\theta_{i}$ | 0.0 | degree
$\mu_{\rm w}$ | 0.0015 | Pa.s | $\mu_{\rm nw}$ | 0.0015 | Pa.s
$r_{\rm av}$ | 23 |  $\mu$m | L | 0.001 | m
$\eta$ | 1.5 | [-] |  |  |
Table 1: Parameters used to simulate quasi-static fluid displacement in BoT.

These parameters are coupled to the BoT model to simulate capillary pressure
curves.

###  3.2 Static capillary pressure for end wetting states

A starting point for the dynamic capillary-pressure models presented in
Section 2 (Equations (1) and (4)) is characterizing the capillary pressure
curves for the end wetting states. Given the same tube geometry and fluid
pairing described above, the capillary pressure functions are simulated under
static conditions for both the initial and final wetting states using the BoT
model.

Only a single drainage experiment is needed in the static case to fully
characterize the capillary pressure curve. This is due to the simplicity of
the BoT model and lack of residual saturation in cylindrical tubes. Therefore,
hysteresis is not possible if the contact angles in the tubes (and other
parameters) are held constant.

The simulated static curves can be correlated with the Brooks-Corey model (2).
The resulting correlations can be found in Figure 3, while fitted parameters
for the Brooks-Corey model can be found in Table 2.

Initial wetting state |  |  Final wetting state |  |
---|---|---|---|---
param. |  value |  unit |  param. |  value |  unit
$c_{w}$ |  360 |  [Pa] |  $c_{w}$ |  56 |  [Pa]
$a_{w}$ |  0.2778 |  [-] |  $a_{w}$ |  0.2778 |  [-]
$\mathrm{R}^{2}$ |  1 |  - |  $\mathrm{R}^{2}$ |  1 |  -
Table 2: Estimated correlation parameter values for initial and final wetting
state capillary pressure curves

Figure 3 compares the Brooks-Corey formula (2) and the capillary pressure
curves associated with static contact angles (initial and final wetting
states).

Figure 3: Simulated $P_{\textrm{c}}$-$S$ data for the initial and final
wetting states compared to a Brooks-Corey model with calibrated parameters
given in Table 2.

The Brooks-Corey correlation gives an excellent match to the simulated
$P_{\textrm{c}}$-$S$ data under static conditions. We observe that the pore-
size distribution index $a_{w}$ for the initial and final wetting states are
the same for both end wetting states, which is expected since the same
distribution of tube radii is used in both cases. While, on the other hand,
the coefficient $c_{w}$ decreases by a factor of 0.85 from the initial to
final wetting state, corresponding to a decrease in macroscale capillary entry
pressure. Leverett-J scaling theory [34] predicts that entry pressure scales
by $\cos\theta$, which agrees nicely with the reduction in $\cos\theta$ by a
factor of 0.83 for a contact angle change from 0 to 80 degrees.

We reiterate that for the static case where no wettability alteration occurs,
the Brooks-Corey model describes _both_ drainage and imbibition for the bundle
of tubes.

###  3.3 Dynamic capillary pressure for uniform wettability alteration

Figure 4: Simulation results for WA-induced dynamics of capillary pressure as
a function of wetting saturation given uniform wettability alteration at the
pore scale. The color of each data point indicates the time elapsed in months.
The static $P_{\textrm{c}}$-$S$ curves for the end wetting states,
$\theta_{i}$ and $\theta_{f}$, are plotted as a reference (the black broken
lines).

Here, we present results for the first dynamic case (Figure 4) in which
wettability is altered uniformly throughout the bundle by considering the
pore-scale parameter $C=0.005$. Two dynamic drainage-imbibition cycles are
simulated using the quasi-static BoT approach over a period of approximately 7
months. During this time, the capillary pressure steadily decreases from the
initial to the final wetting state with each subsequent cycle. Having reached
the final wetting state, any additional drainage-imbibition cycle would follow
along the static curve for the final wetting state. We observe that
wettability-induced dynamics also introduces an apparent hysteresis in the
$P_{\textrm{c}}$-$S$ data. We recall that a simple BoT model cannot exhibit
hysteresis under static wettability conditions.

Figure 5: Uniform contact angle $\theta$ change across the bundle as a
function of drainage-imbibition fluid history paths. The color of each data
point indicates the time elapsed in months.

The dynamics in capillarity are coupled to a continuous change in contact
angle between the two wetting states, which is shown in Figure 5 for the same
two cycles. We observe that $\theta$ undergoes a rapid increase by the end of
the first drainage, followed by a slower evolution towards the final wetting
state. Since there is no spatial variability in wettability at the pore scale
in the uniform case, smaller tubes are altered at the same rate as larger
tubes. This is evident in the associated $P_{\textrm{c}}$-$S$ data where we
observe substantial alteration of the first drainage path at low wetting
saturation. The alteration from the initial wetting condition is fast, leading
to non-trivial and almost non-monotonic dynamics in capillarity at early time.

The saturation history is an important factor, defining the length of exposure
in individual tubes. Figure 6a shows that saturation history is driven by the
time required to displace fluid through each tube in sequence, i.e. from
largest to smallest during drainage and vice versa during imbibition. From
integration of this curve, we obtain the exposure time $\overline{\chi}$
(Figure 6b) which controls the contact angle change shown previously in Figure
5. We recall that for the uniform model, contact angle change at the pore-
scale (Equation (10)) is driven by the average saturation not the local
saturation.

Figure 6a also indicates the time taken for each drainage-imbibition cycle,
with the first cycle takes approximately 1 month to complete. Since the
displacement process is speeding up with increasing contact angle, up to 12
cycles would be required in this manner before the end wetting state is
reached. In the interest of presentation, we have controlled the pressure drop
during the final imbibition experiment so that exposure time was lengthened in
order to reach the final wetting state by the end of the second cycle.

Figure 6: Saturation history for the simulated $P_{\textrm{c}}$-$S$ data over
two drainage-imbibition cycles as shown in Figure 4 for the uniform WA case:
(a) average saturation history of the non-wetting fluid with respect to time
in months, and (b) the integrated exposure time $\overline{\chi}$ over the
same period.

We have observed above that capillary pressure is no longer a unique function
of saturation as it is in the static case for this simple geometry. However,
if $P_{\textrm{c}}$-$S$ is transformed to the temporal domain by plotting
against $\overline{\chi}$, we obtain a unique function with respect to
exposure time. However, Figure 7a shows that the capillary pressure relates
with the temporal domain $\overline{\chi}$ non-monotonically i.e., the
capillary pressure increases and decreases with each drainage-imbibition
cycle. Similarly, the associated contact angle change is plotted along with
exposure time (Figure 7b), which is exactly the function given in Equation
(10) at the pore level. The value of $\theta$ is smoothly increasing in time
according to the pore-scale model.

Figure 7: The capillary pressure data (a) and contact angle (b) as a function
of $\overline{\chi}$ for the uniform WA case. The color of each data point
indicates the time elapsed in months.

####  3.3.1 Dynamic capillary pressure model development

Following the approach discussed in Section 2, we applied Equation (3) to
determine the dynamic coefficient $\omega$. As shown in Figure 8a, we observe
that $\omega$ is a non-unique function of wetting-phase saturation but with
values that are continuously increasing as the dynamic capillary pressure
approaches the final wetting state. This non-uniqueness behavior makes
challenging to propose a functional form for $\omega$-$S_{w}$ relation.

Figure 8: Plot of the dynamic coefficient $\omega$ against wetting phase
saturation (a) and $\overline{\chi}$ (b) for the uniform wettability
alteration case. The data points are color-coded with exposure time in months.

On the other hand, Figure 8b shows that $\omega$ is smoothly increasing,
uniquely related to $\overline{\chi}$, and mimics the functional form of the
pore-scale model in Equation (9) (see Figure 7b). Given this insight, we
propose an adsorption-type model to fit the interpolation coefficient $\omega$
to the average exposure time:

$\omega=\frac{\overline{\chi}}{\beta_{1}+\overline{\chi}},$ (14)

where $\beta_{1}$ is a fitting parameter obtained from the best fit to
simulated data in Figure 8b. For this particular case, the calibrated
parameter is estimated to be $\beta_{1}=0.01$.

The proposed model for $\omega$ in Equation (14) is incorporated into the
dynamic capillary pressure model proposed in Equation (4) and compared against
the simulated dynamic capillary pressure curves. The resulting fit, shown in
Figure 9 indicates that a one-parameter model for the dynamic coefficient
agrees well with the simulated data. The correlation coefficient for this
comparison is $R^{2}=0.9921$.

Figure 9: Comparison of the dynamic capillary pressure model for uniform
wettability alteration with simulated data: as a function of wetting
saturation (a) and $\overline{\chi}$ (b).

We hypothesize that the value of $\beta_{1}$ in this model is dependent on the
dynamics of wettability alteration at the pore-scale. Therefore, we
investigate how the parameter $\beta_{1}$ evolves with changes in the pore-
scale parameter $C$ in Equation (9). We repeat the capillary pressure
simulations for different values of the pore-scale parameter $C$ and determine
the correlated value of $\beta_{1}$ in each case. Figure 10 shows that the
interpolation model parameter is directly proportional to the pore-scale model
parameter, with a proportionality constant of 2. Thus, the relationship
$\beta_{1}=2C$ can be used to predict the macroscale parameter directly from
knowledge of the pore-scale process.

The final dynamic capillary pressure model can now be obtained by
incorporating the above model for $\omega$ with the established relationship
for $\beta_{1}$ into the general form of the interpolation model in Equation
(4):

$P_{c}=\frac{\overline{\chi}}{2C+\overline{\chi}}\Big{(}P_{c}^{\rm
st,f}-P_{c}^{\rm st,i}\Big{)}+P_{c}^{\rm st,i}.$ (15)

In its final form, the dynamic capillary pressure model in Equation (15) is
solely a function of the pore-scale parameter $C$. This parameter must be
determined from laboratory experiments for a given sample exposed to a WA
agent.

Figure 10: The relation between pore-scale wettability model parameter $C$ and
the dynamic coefficient parameter $\beta_{1}$.

####  3.3.2 Dynamic model applicability

We note that the saturation history in Figure 6a used to generate the
$P_{\textrm{c}}$-$S$ curves in Figure 4 can be thought of as one arbitrary
path of an infinite number of possible paths. If a different path had been
chosen, such as a flow reversal at intermediate saturation or a prolonged
exposure time at a given saturation, it would result in entirely different
capillary pressure dynamics.

We generate many different $P_{\textrm{c}}$-$S$ curves by taking numerous
different paths in the saturation-time domain. The resulting simulated data
forms a surface with respect to saturation and exposure time as shown in
Figure 11a. We then apply the calibrated dynamic model from Equation (15) to
the same saturation-time paths used to generate the
$P_{c}$-$S$-$\overline{\chi}$ surface. The difference between the calibrated
model and the simulated data is shown in Figure 11b. A good comparison of the
dynamic model to simulated data demonstrates that model calibration to a
single saturation-time path is robust enough to be applied to any possible
path.

Figure 11: (a) Simulated capillary pressure under uniform WA obtained by
taking multiple paths in the $S_{w}\times\overline{\chi}$ domain. (b)
Difference between the dynamic capillary pressure model prediction (15) with
the data in (a).

###  3.4 Dynamic capillary pressure for non-uniform wettability alteration

We present the second case where WA evolves throughout the bundle non-
uniformly at the pore level with the pore-scale parameter $C=1\times 10^{-5}$.
In this case, we performed four WA induced dynamic drainage-imbibition
displacements and are shown in Figure 12.

Figure 12: Simulated capillary pressure curves for the case of non-uniform WA
over four drainage-imibition cycles with respect to wetting phase saturation.
The static $P_{\textrm{c}}$-$S$ curves for the initial and final wetting
states are shown as a reference. The color code shows the $P_{\textrm{c}}$-$S$
dynamics along approximatly 33 months of exposure time.

The first drainage curve follows exactly along the static curve for the
initial wetting state, which is expected since the entry pressure to each
individual tube is not altered before the tube has drained. Only after the
tubes are drained by the WA agent can the contact angle change in each tube.
Therefore, it is only the imbibition data that show the corresponding
reduction in capillarity. At the end of imbibition, the system becomes quasi-
static in the sense of alteration, which means the wettability of a given tube
stops being further altered until a new drainage cycle has begun. This implies
that each new drainage curve must follow along the previous imbibition curve
and thus, imbibition and drainage paths overlap in Figure 12. This effect is
in contrast to the uniform WA case, where each drainage and imbibition curves
follow different paths.

The capillary pressure instability in Figure 12 is attributed to the local
change of wettability in each tube. The larger tubes drain first and are
imbibed last, resulting in a long local exposure time compared to the smaller
tubes that drain last and imbibe first. Therefore, the initial wetting state
persists in the smaller tubes. Figure 13 shows the resulting heterogeneous
contact angle change in time across the bundle.

Figure 13: Dynamic contact angle evolution as a function of exposure time to
the WA agent per each tube. This contact angle distribution was recorded at
the end of each drainage-imbibition cycle and the color code shows the contact
angle dynamics during 33 months of exposure time.

We observe that the impact of fluid history and exposure time on the wetting
property of the bundle and the associated capillarity is distinct from the
uniform WA case. Here, the dynamic capillary pressure continues to rise to the
initial wetting state capillary pressure curve for low water saturation. This
is due to the reason that the smaller tubes were exposed to the WA agent for
shorter time than the larger tubes. This situation persists until the bundle
approaches a uniform wetting condition corresponding to the final wetting
state.

As we noticed for the uniform WA case, we also observe $P_{\textrm{c}}$-$S$
hysteresis induced by non-uniform WA in Figure 12. Therefore, we can transform
the $P_{\textrm{c}}$-$S$ curve to the temporal domain $\overline{\chi}$ in the
same way as before. The temporal domain $\overline{\chi}$ is shown in Figure
14a, and is obtained by integrating the averaged non-wetting fluid history in
Figure 14b over exposure period of 33 months. For the sake of clarity and the
same reason as the uniform WA case, we controlled the pressure drop after the
third cycle so that the tubes experience prolonged exposure time in order to
complete the alteration processes within few number of cycles.

Figure 14: Non-wetting fluid history paths over the four displacement cycles
(a) and the averaged variable $\overline{\chi}$ which obtained by integration
of the fluid history in (a) over the exposure period (b).

The averaged variable $\overline{\chi}$ calculated from the saturation history
is then used to transform the capillary pressure curve into temporal space. In
Figure 15, we observe a unique relation between the capillary pressure and
$\overline{\chi}$ as hypothesised above. However, unlike the
$P_{\textrm{c}}$-$S$ relation, the $P_{\textrm{c}}$-$\overline{\chi}$ relation
is not monotonic which rises and falls along exposure time.

Figure 15: Simulated capillary pressure curves for the case of non-uniform WA
over four drainage-imibition cycles with respect to temporal domain
$\overline{\chi}$. The color code shows the dynamics of capillarity during the
exposure time of approximatly 33 months.

####  3.4.1 Dynamic capillary pressure model development

Given the capillary pressure results above, the dynamic coefficient $\omega$
is calculated according to Equation (3), and plotted as a function of $S_{nw}$
and $\overline{\chi}$ in Figure 16a and 16b respectively. Unlike the uniform
WA case, seen in Figure 8b, the coefficient $\omega$ is not monotonically
increasing in $\overline{\chi}$, but instead it rises and falls. This makes it
challenging to propose a functional relation between the dynamic coefficient
$\omega$ and $\overline{\chi}$ directly as we did for the uniform WA case.

Figure 16: The dynamic coefficient obtained by applying equation (3): as a
function of saturation (a) and $\overline{\chi}$ (b). Here also the data
points are color-coded with exposure time in months.

On the other hand, $\omega$ in Figure 16a has a well-behaved curvature for
each drainage-imbibition cycle along the saturation history. Further, the
dynamic coefficient as a function of saturation is increasing with increased
exposure time. Given these insights, we proposed a model for the dynamic
coefficient,

$\omega(S_{nw},\overline{\chi})=\frac{S_{w}}{\alpha(\overline{\chi})+S_{w}},$
(16)

where $S_{w}=1-S_{nw}$, and $\alpha$ controls the curvature of the
$\omega$-$S_{nw}$ curve for each drainage-imbibition cycle. Since $\omega$ is
increasing function of exposure time, $\alpha$ should decrease along the
averaged variable $\overline{\chi}$.

The function form of $\omega$ in (16) is then matched with the
$\omega$-$S_{nw}$ data to analyze the dynamics of $\alpha$ along
$\overline{\chi}$. The obtained $\alpha$-$\overline{\chi}$ relation is
decreasing as hypothesized and in particular has a form that can be governed
as,

$\alpha(\overline{\chi})=\beta_{2}/\overline{\chi},$ (17)

where $\beta_{2}$ is non-dimensional fitting parameter. For this particular
simulation the parameter $\beta_{2}$ is estimated to be 0.004 for the four of
drainage-imbibition cycles.

The suggested models (16) and (17) are substituted into the interpolation
model (4) to give:

$P_{c}=\frac{\overline{\chi}S_{w}}{\beta_{2}+\overline{\chi}S_{w}}\Big{(}P_{c}^{\rm
st,f}-P_{c}^{\rm st,i}\Big{)}+P_{c}^{\rm st,i}.$ (18)

The calibrated dynamic capillary pressure model is compared with the simulated
time-dependent capillary pressure data and depicted in Figure 17. We observe
that the proposed dynamic model agrees well with simulated dynamic capillary
pressure curves. Thus, we have obtained a single valued single parameter model
that describes the evolution of dynamic capillarity over multiple drainage-
imbibition cycles rather than using a model consisting of multiple parameters
that change with each cycle (or hysteresis models).

Figure 17: Comparison of the dynamic model (18) with the simulated
$P_{\textrm{c}}$-$S$ data for the case of non-uniform wettability alteration.
The data points are color-coded with exposure time in months

####  3.4.2 Dynamic model applicability

As we point out in the uniform WA case, the saturation paths above are
arbitrary paths out of infinitely many $P_{\textrm{c}}$-$S$-$\overline{\chi}$
paths. Keeping this in mind, we generate a data for dynamic capillary pressure
by considering possibly many paths in $S_{w}\times\overline{\chi}$ domain,
resulting in a capillary pressure surface shown in Figure 18a.

Figure 18: (a) Simulated capillary pressure under non-uniform WA obtained by
taking multiple paths in the $S_{w}\times\overline{\chi}$ domain. (b) is the
difference between the dynamic capillary pressure (18) and the
$P_{\textrm{c}}$-$\overline{\chi}$ data in (a).

Following the paths used to generate the simulated
$P_{\textrm{c}}$-$S$-$\overline{\chi}$ surface, we applied the calibrated
model (18) to predict the simulated capillary pressure data. Figure 18b shows
the actual difference between the simulated data and the model predicted data.
The comparison demonstrates that a single parameter value model calibration to
a single arbitrary path in $S_{w}\times\overline{\chi}$ domain can be applied
to any possible path in the domain. In general, this simple model is efficient
to predict the the impact of non-uniform WA at the pore scale on the capillary
pressure at the macroscale for any given exposure time and saturation path.

Figure 19: The relation between pore-scale wettability model parameter $C$ and
the dynamic coefficient parameter $\beta_{2}$ in Equation (18).

Given the successful WA induced dynamic capillary pressure model, we perform
the same exercise as before where we vary the pore-scale parameter $C$, from
Equation (9), and observe its impact on the macroscale parameter $\beta_{2}$
in the interpolation model above. The relation between those parameters is
obtained and depicted in Figure 19. A power law model,
$\beta_{2}=b_{1}C^{b_{2}}$ is observed between these parameters and the model
parameters are estimated to be $b_{1}=3.3\times 10^{6}$ and $b_{2}=1.8$. The
model for $\beta_{2}$ can be substituted into model (18) to characterize the
dynamics of the capillary pressure directly from the pore scale phenomenon.

To summarize, we have obtained a reliable model to describe the impact of the
pore-scale wettability alteration on capillary pressure dynamics.

##  4 Discussion and Conclusions

We have proposed a new interpolation based model for dynamic capillary
pressure caused by time-dependent wettability alteration. The model has been
correlated to pore-scale simulations of capillary pressure-saturation curves
where a time-dependent model for contact angle change has been applied
directly at the pore scale. The resulting correlation shows that a simple
interpolation-based model is highly effective at capturing WA induced dynamics
in macroscale capillary pressure by a single value parameter. In the interest
of completeness, we explored other types of models to capture capillary
pressure dynamics, including the mixed-wet model of Skjæveland et al [33]. For
brevity, we do not report the results of that separate study herein. We found
that although other models could be calibrated with reasonable accuracy, they
all involved more than one calibration parameter (up to four) that need to be
adjusted in each drainage-imbibition cycle. Therefore, the single-parameter
single-valued interpolation model presented in this study is the preferred
model due to its reliability.

The proposed interpolation model is a macroscale model that allows for change
in capillary pressure as a function of macroscale variables–saturation and
exposure time to a WA agent. We recall that exposure time is simply the
integration of saturation history over time. The model consists of three main
components–two capillary pressure functions at the initial and final wetting
state and a dynamic interpolation coefficient that moves from one state to the
other. The initial and final capillary pressure functions can be determined _a
priori_ from static experiments using inert fluids. In this study, the initial
and final states are represented by classical Brooks-Corey functions. The
dynamic coefficient is thus the only variable correlated to dynamic capillary
pressure simulations. In this study, we have shown that the coefficient can be
easily correlated to saturation and exposure time via a single parameter.

The resulting correlation is dependent on the conditions for wettability
alteration. Two conditions are explored, the first is dissolution of the agent
into the wetting phase that alters the contact angle uniformly across the REV.
The second is contact angle change only through direct contact with an
altering agent as a separate non-wetting phase, resulting in non-uniform
contact angle change throughout the REV. The differences in the two
wettability alteration mechanisms changes the complexity of the resulting
capillary dynamics. In the uniform case, the dynamic coefficient can be
correlated to exposure time through a sorption-type model, which seems to be a
natural result given the contact angle change at the pore scale is also based
on a sorption model. In the non-uniform case, the dynamic coefficient has a
similar sorption form with increased exposure time, but now with the product
of saturation and exposure time as the dynamic variable. The additional
complexity is needed to draw the capillary pressure curve back to the initial
wetting at low saturation (a region of the $P_{\textrm{c}}$-$S$ curve
dominated by smaller pores where the contact angle takes longer to change).

An important result of this study is quantifying the link between the pore-
and macroscale. We showed that by varying the parameter that alters the speed
and extent of contact angle change in each individual pore, we could predict
the resulting impact on dynamic capillary pressure. In fact, in both the
uniform and non-uniform cases, there is a very simple scaling from the pore-
scale and macroscale parameters. In the uniform case, the two parameters are
directly proportional, while in the non-uniform case, the macroscale parameter
scales with the pore-scale parameter via a power law. The implication of this
result is that by knowing the mechanism that controls contact angle at the
pore-scale, which can be obtained by a relatively simple batch experiment, we
can quantify _a priori_ the macroscale dynamics without having to perform
pore-scale simulations. This is an important generalization and valuable for
making use of experimental data to inform macroscale constitutive functions.

The success of the interpolation model in this study is a natural development
from previous studies that incorporates the interpolation model directly into
reservoir simulation of wettability alteration. In those studies, the model
was matched directly to macroscale data in a heuristic manner. The
contribution of this study is to quantify the pore-scale underpinnings to the
interpolation model through a direct and systematic manner, thus confirming
the validity of this type of model for use in macroscale simulation. The
insights gained in this study into capillary pressure dynamics due to
wettability alteration. Further study is needed to calibrate the macroscale
dynamic parameter to actual data.

There are additional aspects that were not addressed in this study. First, is
whether the effect of time-dependent wettability alteration can be quantified
in a similar manner for relative permeability. For brevity, we did not explore
this additional constitutive relationship in this paper. This can be explored
in future work. An additional quantity of interest is residual saturation.
Quantifying dynamics and stability of residual wetting and non-wetting
saturation is an active area of research. The scope of this work involved a
bundle-of-tubes pore-scale model, which due to its simple nature cannot model
trapped phases. Therefore, a more realistic pore-network model would need to
be implemented, which is the subject of ongoing work.

The results of this study have important implications for future laboratory
investigations. Previous experimental studies have shown that wettability
alteration can cause a significant decrease in capillarity that occurs at the
same time scale as the experiment. Our work further highlights this time-
dependent effect and shows that proper quantification of the dynamics requires
capillary pressure data from both inert and reactive fluid pairings. In
addition, there should be separate batch measurements on mineral surfaces to
characterize the mechanistic model for contact angle change at the pore scale,
such as the one proposed herein. Also important is whether the alteration
mechanism occurs by dissolution of a WA agent or by direct contact with a non-
wetting phase. Such data are important for determining the macroscale
parameters needed in the interpolation-based model.

At the reservoir scale, there is a potential for dynamics in capillary
pressure to have a significant impact on fluid flow that should be accounted
for during simulations. For instance, for CO${}_{\textrm{2}}$ storage could
there be a change in wetting that could degrade the integrity of the sealing
caprock and allow for seepage of CO${}_{\textrm{2}}$ out of the storage
reservoir, and thus further investigation is warranted. This will require
implementation of dynamic capillary pressure functions into reservoir
simulation, with possible wettability reversal that could prove challenging
for numerical consistency and robustness. Modified simulation methods can then
be used in the future to quantify the impact of capillary dynamics on
reservoir processes for both hydrocarbon recovery and CO${}_{\textrm{2}}$
storage.

## Acknowledgement

The research is supported by CHI project funded by the CLIMIT program of the
Research Council of Norway.

## References

  * [1] Bonn, D. & Eggers, J. & Hindekeu, J. & Meunier, J. & Rolley, E. 2009 Wetting and spreading. Rev. Mod. Phys. 81, 739-805.
  * [2] Iglauer, S. & Rahman, T., & Sarmadivaleh, M. & Al-Hinai, A. & Fernø, M. A. & Lebedev, M. 2016 Influence of wettability on residual gas trapping and enhanced oil recovery in three-phase flow: A pore-scale analysis by use of microcomputed tomography. SPE J. 21, 1916–1929.
  * [3] auYu, L. & Kleppe, H. & Kaarstad, T. & Skjæveland, S. M. 2008 Modelling of wettability alteration processes in carbonate oil reservoirs. Netw. Heterog. Media 3, 149–183.
  * [4] Iglauer, S. & Pentland, C. H. & Busch, A. 2014 CO${}_{\textrm{2}}$ wettability of seal and reservoir rocks and the implications for carbon geo-sequestration. Water Resour. Res.51, 729–774.
  * [5] Blunt, M. J.2001 Flow in porous media– pore-network models and multiphase flow. Curr. Opin. Colloid Interface Sci. 6, 197–207.
  * [6] Ahmed, A. & Patzek, T. W. 2003 Impact of Wettability Alteration on Two-Phase Flow Characteristics of Sandstones: A Quasi-Static Description. Water Resour. Res. 39, 1-12.
  * [7] Blunt, M. J. 1997 Pore Level Modeling of the Effects of Wettability. SPE J. 2, 494–510.
  * [8] Morrow, N. R. & Lim, H. T. & Ward, J. S. 1986 Effect of crude-oil-induced wettability changes on oil recovery. SPE J., 89-103
  * [9] Buckley, J. S. & Liu, Y. & Monsterleet, S. 1988 Mechanisms of wetting alteration by crude oils. SPE J. 3, 54-61.
  * [10] Jadhunandan, P. P. & Morrow, N. R. 1995 Effect of Wettability on Waterflood Recovery for Crude-Oil/Brine/Rock Systems. SPE Reservoir Engineering 10, 40–46.
  * [11] Haagh, M. E. J. & Siretanu, I. & Duits, M. H. G. & Mugele, F. 2017 Salinity-Dependent Contact Angle Alteration in Oil/Brine/Silicate Systems: the Critical Role of Divalent Cations. Langmuir 33, 3349–3357.
  * [12] Singh, R. & Mohanty, K. 2016 Foams with wettability-altering capabilities for oil-wet carbonates: A synergistic approach. SPE J. 21, 1126-1139.
  * [13] Kim, Y. & Wan, J. & Kneafsey, T. J. & Tokunaga, T. K. 2012 Dewetting of silica surfaces upon reactions with supercritical CO${}_{\textrm{2}}$ and Brine: Pore-scale studies in micromodels.Environ. Sci. Technol.46, 4228–4235.
  * [14] Tokunaga, T. K. & Wan, J.  2013 Capillary pressure and mineral wettability influences on reservoir CO${}_{\textrm{2}}$ capacity. Rev. Mineral. Geochem. 77, 481-503.
  * [15] Chiquet, P. & Broseta, D. & Thibeau, S. 2007a Wettability alteration of caprock minerals by carbon dioxide. Geofluids, 112–122.
  * [16] Chalbaud, C. & Robin, M. & Lombard, J. & Martin, F. & Egermann, P. & Bertin, H. 2009 Interfacial tension measurements and wettability evaluation for geological CO${}_{\textrm{2}}$ storage. Adv Water Resour 32, 98–109.
  * [17] Plug, W. J. & Bruining, J. 2007 Capillary pressure for the sand-CO${}_{\textrm{2}}$-water system under various pressure conditions. Application to CO${}_{\textrm{2}}$ sequestration. Adv Water Resour 30 2339–2353.
  * [18] Wang, S. & Tokunaga, T. K. 2015 Capillary pressure–saturation relations for supercritical CO${}_{\textrm{2}}$ and brine in limestone/dolomite sands: Implications for geologic carbon sequestration in carbonate reservoirs. Environ. Sci. Technol. 49, 7208–7217.
  * [19] Wang, S. & Tokunaga, T. K. & Wan, J. & Dong, W. & Kim, Y. 2016 Capillary pressure-saturation relations in quartz and carbonate sands: Limitations for correlating capillary and wettability influences on air, oil, and supercritical CO${}_{\textrm{2}}$ trapping. Water Resour. Res., 6671–6690.
  * [20] Tokunaga, T. K. & Wan, J. & Jung, J. & Kim, T. W. & Kim, Y. & Dong, W. 2013 Capillary pressure and saturation relations for supercritical CO${}_{\textrm{2}}$ and brine in sand: High-pressure $P_{c}(S_{w})$ controller/meter measurements and capillary scaling predictions. Water Resour. Res., 49, 4566–4579.
  * [21] Hassanizadeh, S. & Celia, M. & Dahle, H.2002 Dynamic effects in the capillary pressure–saturation relationship and its impacts on unsaturated flow. Vadose Zone J 1, 38–57.
  * [22] Dahle, H. K. & Celia, M. A. & Hasanizadeh, S. M. 2005 Bundle-of-Tubes Model for Calculating Dynamic Effects in the Capillary-Pressure Saturation Relationship. Transport porous med 58, 5–22.
  * [23] Krumpfer, J. W. & McCarthy, T. J. 2010 Contact angle hysteresis: a different view and a trivial recipe for low hysteresis hydrophobic surfaces. Faraday Discuss. 146, 103–111.
  * [24] Eral, H. B. & ’t Mannetje, D. J. C. M. & Oh, J. M. 2013 Contact angle hysteresis: a review of fundamentals and applications. Colloid Polym Sci 291, 247–260.
  * [25] Delshad, M. & Najafabadi, N. F. & Anderson, G. A. & Pope, G. A. & Sepehrnoori, K. 2009 Modeling Wettability Alteration by Surfactants in Naturally Fractured Reservoirs. SPE J. 12, 361-370.
  * [26] Lashgari, H. R. & Xu, Y. & Sepehrnoori, K. 2016 Modelling dynamic wettability alteration effect based on contact angle. SPE, 1–17.
  * [27] Andersen, P. Ø. & Evje, S. & Kleppe, H. & Skjæveland, S. M. 2015 A Model for Wettability Alteration in Fractured Reservoirs. SPE J. 20, 1261–1275.
  * [28] Adibhatia, B. & Sun, X. & Mohanty, K. 2005 Numerical Studies of Oil Production from Initially Oil-Wet Fracture Blocks by Surfactant Brine Imbibition. SPE, 1-15.
  * [29] Al-Mutairi, S. M. & Abu-Khamsin, S. A. & Hossain, M. E. 2012 A Novel Approach to Handle Continuous Wettability Alteration during Immiscible CO2 Flooding Process. SPE, 1-12.
  * [30] Bartley, J. T. & Ruth, D. W. 1999 Relative Permeability Analysis of Tube Bundle Models.Transport porous med 36, 161–187.
  * [31] Helland, J. O. & Skjæveland, S. M. 2006 Physically based capillary pressure correlation for mixed-wet reservoirs from a bundle-of-tubes model. SPE, 11, 171–180.
  * [32] Dong, M. & Dullien, F. A. L. & Dai, L. & Li, D. 2005 Immiscible Displacement in the Interacting Capillary Bundle Model Part I. Development of Interacting Capillary Bundle Model. Trends Anal. Chem. 59, 1–18.
  * [33] Skjæveland, S. M. & Siqveland, L. M. & Kjosavik, A. & Thomas, W. L. H. & Virnovsky, G. A. 2000 Capillary pressure correlation for mixed-wet Reservoirs. SPE, 60–67.
  * [34] Xu, W. S. & Luo, P. Y. & Sun, L. & Lin, N. 2016 A Prediction model of the capillary pressure J–function. Plos One 11, 1–9.


∎

11institutetext: Sowmini Devi Veeramachaneni, Vineet Padmanabhan
22institutetext: School of Computer and Information Sciences, University of
Hyderabad, Hyderabad, India
Tel.: +91-8985061848
22email<EMAIL_ADDRESS><EMAIL_ADDRESS>33institutetext: Arun K Pujari 44institutetext: Mahindra University Ecole
Centrale School of Engineering (MEC), Hyderabad, India
44email<EMAIL_ADDRESS>55institutetext: Vikas Kumar 66institutetext:
University of Delhi, Delhi, India
66email<EMAIL_ADDRESS>
# A Hinge-Loss based Codebook Transfer for Cross-Domain Recommendation with
Nonoverlapping Data

Sowmini Devi Veeramachaneni  Arun K Pujari  Vineet Padmanabhan  Vikas Kumar

(Received: date / Accepted: date)

###### Abstract

Recommender systems(RS), especially collaborative filtering(CF) based RS, has
been playing an important role in many e-commerce applications. As the
information being searched over the internet is rapidly increasing, users
often face the difficulty of finding items of his/her own interest and RS
often provides help in such tasks. Recent studies show that, as the item space
increases, and the number of items rated by the users become very less, issues
like sparsity arise. To mitigate the sparsity problem, transfer learning
techniques are being used wherein the data from dense domain(source) is
considered in order to predict the missing entries in the sparse
domain(target). In this paper, we propose a transfer learning approach for
cross-domain recommendation when both domains have no overlap of users and
items. In our approach the transferring of knowledge from source to target
domain is done in a novel way. We make use of co-clustering technique to
obtain the codebook (cluster-level rating pattern) of source domain. By making
use of hinge loss function we transfer the learnt codebook of the source
domain to target. The use of hinge loss as a loss function is novel and has
not been tried before in transfer learning. We demonstrate that our technique
improves the approximation of the target matrix on benchmark datasets.

###### Keywords:

Matrix Factorisation Collaborative Filtering CodebookTransfer LearningCross-
Domain Recommendation

##  1 Introduction

Recent years have witnessed a clear explosion in the amount of e-commerce data
available. As a result, the quantity of information that needs to be searched
for suggesting items of interest to users is becoming more challenging.
Recommender systems 10.5555/2931100 BOBADILLA2013109 10.5555/1941884  can be
used for predicting items that the users will be interested in and has become
a core part of several e-commerce applications. Machine learning techniques
article  have become an integral part in the process of building recommender
systems and in the last decade or so there has been a surge in this area of
research. In general, the existing algorithms for recommender systems can be
broadly categorized into three major categories namely content-based (CB),
collaborative (CF) and hybrid recommendation. In content-based method, the
recommendation is made by first calculating a similarity score between a user
profile and the profile information of an item. Based on the similarity score,
a set of top $k$ items are then recommended to a user. For instance, in movie
recommendation task, the profile information of a user contains his/her
interest over the genre. A movie is then can be recommended by finding out the
similarity between the user interest for the genre and the movie description
(genre). In the case of collaborative filtering, the preferences of the target
user are matched with that of other users having similar tastes and a
recommendation is made. For example, if a set of users likes the same movies
as that of user $A$ (say target user), then there is a chance that user $A$
will like the movies liked by the set of users which the user $A$ has not
seen. The hybrid method combines both the content-based and collaborative-
based models.

The major shortcoming with content based recommendation is the situations
where the item features are not meaningful or situations where there is a need
to capture the change in user interests over time. On the other hand,
collaborative filtering handles the above mentioned situations as it only
requires the preference (implicit or explicit) information for recommendation.
The rating matrix is memorised in the memory-based methods and the missing
entries are predicted based on the relation (ratings) between similar
users/items. Among memory-based methods, neighbourhood-based methods are
prominent wherein users/items that are similar to target user/item are
considered for prediction. The idea is that any two users having similar
ratings for few items will also tend to have similarity in ratings for the
remaining items too. In addition to the assumption, it is also assumed that if
rating similarity is there for two items among some users then other users
will also give the same rating to the items. Identifying similar users as that
of the target user and using the mean of their ratings for prediction is the
crux of _user-based_ CF methods. In the case of _item-based_ CF it is the
other way around. Depending on how the weighted average is computed (Pearson
correlation, Vector Cosine, Mean-squared-difference etc.), neighbourhood
methods may vary.

In the case of model based methods, missing ratings can be predicted more
quickly by making use of a parametric model which is built from the training
data. Low-rank matrix factorisation, which is a latent factor model, is a
successful and an effective model-based CF approach. In the Matrix
Factorisation (MF) approach as proposed in koren2009matrix  and one of the
most widely-used CF technique, a small set of user and item latent factors are
learnt based on the user-item rating matrix. By making use of the matrix
factorisation method as outlined in koren2009matrix  for each user (or, each
item) one latent vector is constructed by extracting the latent factors of
users/items and this vector captures the user’s (or, item’s) characteristics.
To obtain the user ratings for an item, the product of a user latent vector
with that of an item is made use through a process called low-dimensional
embedding. Maximum Margin Matrix Factorisation (MMMF) srebro2004maximum  is
perhaps the most widely used and successful matrix factorisation technique.

Recommender systems often face with the problem of _sparsity_ which results
from less amount of information being available of the data in hand. Machine
learning techniques like transfer learning pan2010survey  has been applied in
collaborative filtering based recommender systems to avoid the problems
related to sparsity. Transfer learning involves the construction of a
predictive model across multiple domains wherein relevant knowledge from one
domain is extracted and transferred to another domain. Transfer learning in
collaborative filtering-based recommender systems is otherwise called as
cross-domain collaborative filtering. When multiple domains are involved it is
necessary for transfer learning to address crucial issues related to the
presence/absence of aligned users/items during knowledge transfer. It is often
the case that the presence of common users/items across domains do not exist
and even if it does exist it is difficult to map their correspondence.
Transfer learning in the absence of common users/items in the two domains is
often hard to comprehend and _CodeBook transfer_ (CBT) li2009can  has been
proposed as a transfer learning approach to address these type of problems.
Codebook construction is based on the idea that the abundance of information
in one domain is to be made use to predict the missing information in another
domain. Since it is assumed that certain kind of rating pattern is involved in
any kind of data, the codebook tries to capture this rating pattern in a
summarised form. The belief/hypothesis is that across domains these rating
patterns are unvarying.

Not only in recommender systems, but also in many machine learning algorithms,
loss functions plays a significant role in empirical risk minimization and
computational complexities loss1 ; loss2 . So, choice of loss function is very
important. Among all the loss functions, hinge loss is more suitable for
discrete classification task. Also, among all the matrix factorisation
techniques, maximum margin matrix factorisation(MMMF) treats the collaborative
prediction problem as a classification task and takes full advantage of the
hinge loss. Until now MMMF has been used in single domain recommender systems
and in our proposed approach, we adapt it to work for cross domain
recommendation by taking advantage of the hinge loss function.

This paper proposes a codebook based transfer method for collaborative
filtering in cross domain recommender systems. Our proposed approach consists
of two steps: in the first, we construct the codebook by making use of the
source data and apply a technique called co-clustering on it. In the second
step, the idea is to transfer the learned codebook in a novel way to the
target domain by making use of hinge loss as the loss function. In previous
works, researchers have made use of squared loss while transferring the source
knowledge to the target domain in order to predict the missing ratings of
target domain. The proposed codebook based transfer method for collaborative
filtering outperforms MMMF (i.e., when MMMF is directly used on target data)
and other major codebook based transfer learning approaches for CF as can be
seen from the experimental results.

The remainder of the paper is structured in the following manner: Section 2
presents a quick review of current work on transfer learning in recommendation
systems and the proposed approach is explained in Section 3. The experimental
results are given in Section 4, and Section 5 provides the conclusion of the
work.

##  2 Related work

Matrix Factorisation has been very popular and successful in predicting the
missing values of the user-item rating matrix even when the data is too
sparse. The advantage of matrix factorisation technique is that for a matrix
which is very sparse and has less information available, the matrix can be
filled up accurately to a satisfactory level. Although this is quite good as
far as MF techniques for recommender systems are concerned, the major drawback
with having more sparsity is the failure in the accuracy of the prediction
levels. To overcome such scenarios researchers have come up with the idea of
using information from diverse domains so that information of one domain could
be utilized in another. Transfer learning pan2010survey ; pan2016survey  is
one such technique proposed, which makes use of cross domain collaborative
filtering in the domain of recommendation systems to alleviate the sparsity
problem. Transfer learning intents to transfer the information(knowledge) from
a rich(dense) source domain to the sparse target domain. For example, suppose
that a particular user has watched many movies and have rated the same. The
same user has very less ratings in another domain related to books but wants a
book to be recommended, then by using his ratings from the movie domain book
can be recommended. Formally, given a dense source rating matrix (source
domain), and a sparse target rating matrix (target domain), the goal is to
predict the missing entries of target domain by utilizing the information
available in the source domain.

The major question that was unanswered is to determine whether both domains
are suitable for knowledge transfer. This is necessary because transfer of
knowledge can take place only if a correspondence can be established among two
domains. Certain assumptions need to be made for applying transfer learning
strategy which includes having either common subset of items/users or
attribute similarity in items/users. Domains can be linked and the transfer
can happen explicitly via inter-domain similarities, common item attributes,
etc. The transfer can happen also implicitly via shared user latent features
or item latent features or by rating patterns which can be transferred between
the domains.

In Chung:2007:IPR:1282100.1282113 , a framework was proposed in which the
pertinent items in the source domain are picked based on their common
attributes with the target domain (user interested domain). The inter domain
links were created utilising the common item attributes, however there is no
overlap of users/items required between the domains. On the other hand the
transfer of knowledge by the shared latent features (of users/items) is
addressed in pan2010transfer . In this work the latent features of users and
items of source domain are learnt and fitted to a target domain by integrating
the features into the factorisation of target rating matrix via
regularization. However, it require either common users or items among the two
domains. In Pan:2011:TLP:2283696.2283784 , the latent
characteristics(features) of source and target are shared in a collective way.
In this paper a method called matrix factorisation is employed wherein
features of both domains are learned simultaneously rather than learning from
source first and later applying in the target domain. The only constrain is
that in both domains the users and items needs to be identical. The main
assumption that is made in this method is that the change of domain does not
affect the rating behaviour of the same user/item or different users/items
with similar attributes. In order to maximise knowledge transfer, in
zhao2013active ; ZHAO201738  the authors presume that across different domains
there is a correspondence between users/items. Here rather than starting from
source domain, matrix factorisation technique is applied first in the target
rating data so as to learn the latent features of the target. After that based
on certain criteria, random users/items gets selected from the target domain
and is matched with the corresponding users/items from source domain. In order
to ensure that the latent factors of the selected users/items of both the
target and source domain are same, the authors make use of some kind of
regularisation function.

The methods as outlined in itcf ; inproceedings ; PAN201684  takes care of
scenarios wherein the feedback data of target and source are of different
types. The users and items in the target and auxiliary(source) domains are
assumed to be the same in itcf . Two sets of source data are taken into
consideration in tbt , one of which shares a common set of users with the
target data and the other of which shares a common set of items with the
target data. The source data’s latent factors are extracted, and similarity
graphs are constructed from these latent factors. Both latent factors and
similairty graphs are then transferred to the target data.

In collaborative filtering approach that is based on codebook transfer(CBT)
li2009can  the idea is to capture the hidden correlations that exist between
the ratings given by groups of users and groups of items. This pattern is
usually denoted by the name codebook or cluster-level rating pattern. In
codebook transfer(CBT) what happens is that initially, the codebook of the
source domain is extracted by analysing the rating matrix and thereafter this
codebook is used for the prediction of the target domain. Regardless of the
domains, the way the users and items cluster together it so happens that
across domains the cluster level rating behaviour remains invariant.

Figure 1: Illustration of the proposed method

One assumption that is made in this connection is regarding the source matrix
being a full matrix. Cases wherein the rating matrix is not full, the average
rating of a particular user is used for filling-in missing entries of that
user. Thereafter the rows (users) and columns (items) of the filled-in rating
matrix are clustered using co-clustering li2009can . The way the codebook is
generated is such that each user group is assigned a row and each item group
has a column, and each user group and item group has a rating. The average
rating of users and items in a particular group is taken into consideration
for giving a rating for user-group/item-group. In order to use the codebook in
the targt domain it is necessary to select the users/items in the target
domain that are more appropriate with the user-group/item-group of the
transferred codebook.

A linear combination of codebooks from multiple domains is shown in
Moreno:2012:TTL:2396761.2396817  wherein using an optimisation technique the
coefficient of linear combinations is learned. The assumption of having a
fully dense rating matrix is relaxed as demonstrated in
He:2018:RTL:3159652.3159675 ; 8233662 .

A different way of generating the codebook has been outlined in
ji2016improving ; DBLP:journals/asc/VeeramachaneniP19 . In this method, by
making use of the technique of matrix factorisation the user and item latent
factors are generated from the source domain and thereby tends to avoid any
pre-processing and co-clustering of the rating matrix. The user and item
latent factors thus generated are used to obtain the user latent facor group
and item latent factor group. Codebook is generated by multiplying the mean
latent vecores of the group.

In li2009can , the authors focus on extracting the group level behavior of
users on items by assuming that, though the users/items are different across
systems, the groups (groups - based on age, interest etc) of them behave
similarly. All the missing ratings in the source rating matrix are filled in a
priori by a preprocessing step and then co-clustering is applied on a filled-
in source rating matrix to directly get the cluster-level rating pattern
(codebook). The cluster-level rating pattern is then transferred to another
domain (target). In the target domain, the user and item membership to
clusters encode in the codebook so that a user (or, an item) is member of one
user-cluster (or, item-cluster) represented in the codebook. In our approach,
we do not use a separate preprocessing stage in the source domain, and we use
hinge loss instead of squared loss when transferring the learned codebook from
the source to the target domain.

##  3 Proposed Method

Given a dense source user-item rating matrix $X\in\mathbb{R}^{m\times n}$ and
a sparse target user-item rating matrix $Y\in\mathbb{R}^{m^{\prime}\times
n^{\prime}}$ where $m$, $m^{\prime}$ are the number of users in source and
target domain, and $n$, $n^{\prime}$ are the number of items in source and
target data, the goal is to predict the missing entries (as numerous users do
not rate many items) in target domain by leveraging the data from the source
domain. Prediction of missing entries should take place in such a way that the
existing ratings must be approximated with less error rate. The illustration
of the proposed method is shown in Figure 1. Initially (Step-1), the users
(rows) and items (columns) of the source rating matrix need to be
simultaneously clustered (co-clustering) to construct the codebook. We need
the cluster indicators of users and items, and in order to get the cluster
indicators we can choose any of the co-clustering algorithms. We employed the
Orthogonal nonnegative matrix tri-factorisation technique (ONMTF) Ding:2006
to get the cluster indicators, which has been shown to be equivalent to a two-
way k-means clustering approach. The source rating matrix $X$ can be tri-
factorised as follows.

$\min_{P,Q,S\geq 0}||[X-PSQ^{T}]\odot
W||_{F}^{2}+\alpha||P\textbf{1}-1||_{F}^{2}+\beta||Q\textbf{1}-1||_{F}^{2}$
(1)

where $W$ is an indicator matrix and the size of $W$ is $m\times n$. The entry
of $W$ is $1$ if the rating exists in $X$ and $0$ for the rest. $||.||_{F}$
represents the Frobenius norm. The dimension of $P$ is $m\times k_{1}$, $S$ is
$k_{1}\times k_{2}$, and $Q$ is $n\times k_{2}$. $P\textbf{1}-1$ ensures the
row sum of $P$ to be one, similarly with $Q$. Maximum value of each row of $P$
and $Q$ becomes the cluster indicator for the user/item in that row. These $P$
and $Q$ are not in a recognizable form in terms of user/item membership
matrices. To represent $P$ and $Q$, we use binary values by setting the
maximum valued entry in each row to be $1$ and the others to be $0$. As a
result, these binary user/item cluster indicators form (membership) matrices,
denoted by $P_{s}$ and $Q_{s}$, for the source rating matrix. From these
matrices, we can construct the codebook $B$ as follows.

$B=[P^{T}_{s}XQ_{s}]\oslash[P^{T}_{s}11^{T}Q_{s}]$ (2)

$\oslash$ indicates element-wise division. Averaging of all the ratings in
each of the user/item co-cluster takes place in Eq. 2

Once the codebook from the source domain is constructed, transfer (Step-2) it
to the target domain ($Y$) by substituting it in Eq. 3 and solve the objective
function ($\mathcal{J}$) in order to get $U$ and $V$ of target data. We use
Hinge loss to learn $U$ and $V$ instead of squared loss. As the data of ours
is discrete, apart from $U$ and $V$, $r-1$ thresholds $\theta_{ic}$ $(1\leq
c\leq r-1)$ for each user $i$ has to be learned to classify the prediction
value into $r$ discrete values. Here, $r$ is the maximum rating (say 5).

$\mathcal{J}(U,V,\Theta)=\sum_{(i,j)\in\Omega}\sum_{c=1}^{r-1}h(\mathcal{T}_{ij}^{c}(\theta_{ic}-UBV^{T}))+\frac{\lambda}{2}(||U||_{F}^{2}+||V||_{F}^{2})$
(3)

where,

$\mathcal{T}_{ij}^{c}=\begin{cases}-1&\text{if $c<y_{ij}$}\\\ +1&\text{if
$c\geq y_{ij}$}\end{cases}$

$h(\cdot)$ is a smoothed hinge-loss function defined as,

$h(d)=\begin{cases}0&\text{if $d$ $\geq$ 1}\\\ \frac{1}{2}(1-d)^{2}&\text{if
$0<d<1$}\\\ \frac{1}{2}-d&\text{otherwise.}\end{cases}$

$\Omega$ is the set of observed entries, $\lambda>0$ is regularization
parameter. Gradient based approach can be used to solve the optimization
function (Eq. 3) and following are the gradients.

$\frac{\partial\mathcal{J}}{\partial{U_{ik_{1}}}}=\lambda
U_{ik_{1}}-\sum_{c=1}^{r-1}\sum_{j|ij\in\Omega}\mathcal{T}_{ij}^{c}.h^{\prime}(\mathcal{T}_{ij}^{c}(\theta_{ic}-UBV^{T}))VB^{T}$
$\frac{\partial\mathcal{J}}{\partial{V_{jk_{2}}}}=\lambda
V_{jk_{2}}-\sum_{c=1}^{r-1}\sum_{i|ij\in\Omega}\mathcal{T}_{ij}^{c}.h^{\prime}(\mathcal{T}_{ij}^{c}(\theta_{ic}-UBV^{T}))UB$
$\frac{\partial\mathcal{J}}{\partial{\Theta_{ir}}}=\sum_{j|ij\in\Omega}\mathcal{T}_{ij}^{c}.h^{\prime}(\mathcal{T}_{ij}^{c}(\theta_{ic}-UBV^{T}))$

where,

$h^{\prime}(d)=\begin{cases}0&\text{if $d$ $\geq$ 1}\\\ d-1&\text{if
$0<d<1$}\\\ -1&\text{otherwise.}\end{cases}$

In gradient descent algorithms, we start with random $U$, $V$, $\Theta$ and
iteratively update them using the following update rules.

$U_{ik_{1}}^{t+1}=U_{ik_{1}}^{t}-c\frac{\partial\mathcal{J}}{\partial{U_{ik_{1}}^{t}}}$
$V_{jk_{2}}^{t+1}=V_{jk_{2}}^{t}-c\frac{\partial\mathcal{J}}{\partial{V_{jk_{2}}^{t}}}$
$\Theta_{ia}^{t+1}=\Theta_{ia}^{t}-c\frac{\partial\mathcal{J}}{\partial{\Theta_{ia}^{t}}}$

$c$ is the trade-off parameter.

Once $U$, $V$ of target data are learnt, product of $U$, $B$, and $V^{T}$
(i.e., $UBV^{T}$) mapped with the threshold matrix ($\Theta$) becomes the
approximation (prediction$\in\\{1,2,...,r\\}$ say r = 5) of the target rating
matrix (Step-3). This is illustrated in Algorithm 1.

1: Input: A source rating matrix $X$ of size $m\times n$ and a sparse target
rating matrix $Y$ of size $m^{\prime}\times n^{\prime}$ with $y_{ij}$ known
for $(i,j)\in\Omega$

2: Output: $\hat{Y}$ ($y_{ij}$ for $(i,j)\notin\Omega$)

3: Tri-factorise $X$ using Eq. 1 inorder to get $P$, $Q$ and $S$.

4: Calculate codebook $B$ using Eq. 2.

5: Transfer the codebook to the target data by substituting in Eq. 3.

6: Solve the optimization (Eq. 3) using gradient descent technique inorder to
get the target $U$, $V$ and $\Theta$.

7: Product of $U$, $B$, $V^{T}$ (i.e., $UBV^{T}$) by mapping with the
threshold matrix gives $\hat{Y}$.

Algorithm 1 Codebook transfer via hinge loss

##  4 Experimental Results and Analysis

To evaluate the performance of our method we have experimented the method with
different datasets. The datasets used are MovieLens
100K111https://grouplens.org/datasets/movielens/, MovieLens 1M1,
Goodbooks222https://github.com/zygmuntz/goodbooks-10k, Douban
Music333https://github.com/hezi73/TRACER/blob/master/douban.rar, Douban Book3.
From the Goodbooks dataset we have considered the first $5000$ users and
$3000$ items. From the Douban Music and Douban Book data, we have considered
the first $2000$ users and $2000$ items. The entries of all the datasets fall
in {0,1,2,3,4,5}, where $0$ indicates the missing value, $1$ indicates that
the item is leastly liked and $5$ indicates that the item is heavily liked.
The statistics of the datasets are given in Table 1. In all the experiments,
we have divided the observed data into 80% and 20%, in which 80% is used for
training, and 20% is used for testing.

Table 1: Datasets statistics

Dataset | # of Users | # of Items | % of Observed entries
---|---|---|---
MovieLens 100K | 945 | 1682 | 6.29
MovieLens 1M | 6040 | 3952 | 3.77
Goodbooks (Subset) | 5000 | 3000 | 1.08
Douban Music (Subset) | 2000 | 2000 | 11.26
Douban Book (Subset) | 2000 | 2000 | 5.51

###  4.1 Metrics used for Evaluation

As far as collaborative filtering algorithms are concerned the computation of
the prediction accuracy is crucial for assessing the performance of the
algorithms. Root Mean Square Error (RMSE) (Eq. (4)) and Mean Absolute Error
(MAE) (Eq. (5)) are the two most widely used measures for computing the
prediction accuracy. The difference between the predicted rating and the
true(actual) rating forms the basis for these measures and it is natural that
better performance means smaller values for RMSE and MAE. The values stated in
the tables are the average of five runs.

$RMSE=\sqrt{\sum\limits_{(i,j)\epsilon\Omega}\frac{{(y_{ij}-\hat{y}_{ij})}^{2}}{|\Omega|}}$
(4)
$MAE={\sum\limits_{(i,j)\epsilon\Omega}\frac{|(y_{ij}-\hat{y}_{ij})|}{|\Omega|}}$
(5)

where $y_{ij}$ is the original rating and $\hat{y}_{ij}$ is the predicted
rating, $|\Omega|$ is the count of test ratings.

Figure 2: Impact of number of clusters ($k_{1}$, $k_{2}$) on RMSE of MovieLens
1M data when MovieLens 100K is considered as source Figure 3: Impact of number
of clusters ($k_{1}$, $k_{2}$) on MAE of MovieLens 1M data when MovieLens 100K
is considered as source

###  4.2 Comparison Techniques

Some of the methods we consider for evaluating our proposed method are
outlined below;

  * •

MMMF srebro2004maximum ; Jas05 : Maximum Margin Matrix Factorisation (MMMF) is
the dominant factorisation technique used in collaborative filtering. MMMF is
usually applied on the input rating matrix consisting of the user-item
ratings. The idea is to find the user and item latent-factor vectors which are
of low rank by making use of the existing ratings. MMMF can be applied on a
single domain only, and hence in our experiments, we applied it on the target
domain directly.

  * •

MINDTL 8233662 : In MINDTL, codebook is constructed by taking into
consideration the data from all the incomplete source domains. Here codebook
for each domain is constructed. Following that, the constructed codebooks are
linearly integrated and transferred to the target domain, and the
absent(missing) values of the target rating matrix gets predicted. As far as
our experimental setup is concerned only a single domain is taken into
consideration.

  * •

TRACER ZHUANG2018287 : In TRACER, data from multiple domains are accounted for
and based on this, ratings (which includes missing ratings) for all the source
matrices are predicted. Thereafter the predicted knowledge is utilized by
transferring it into the target domain. By making use of consensus
regularisation during the knowledge transfer process, all the predicted values
are forced to be similar. In a way it can be said that in TRACER at the same
time learning and transferring happens. We have considered a single domain in
our experiments, therefore there is no necessity for consensus regularisation.
We thank the authors for giving the code444https://github.com/hezi73/TRACER on
the web.

  * •

CBT li2009can : In this approach, the dense part of the source user-item
rating matrix is considered, and the missing values of the rows of the dense
matrix get imputed using the average of the ratings of particular row (user).
The codebook is obtained from the dense user-item matrix by applying the
technique of co-clustering. In our experiments, unlike in li2009can , which
consider only the dense part of the input data, the codebook is constructed by
making use of the whole source data and the learnt codebook is transferred to
the target domain.

Table 2: Values of Root Mean Square Error and Mean Absolute Error of baseline
methods and proposed method on MovieLens 1M data, Goodbooks data, Douban Book
data

Dataset | Method | RMSE | MAE
---|---|---|---
MovieLens 1M | MMMF | 0.9361 | 0.6402
MINDTL | 0.9948 | 0.7965
TRACER | 0.9800 | 0.8039
CBT | 0.9676 | 0.7746
Proposed | 0.9123 | 0.6134
Goodbooks | MMMF | 0.9604 | 0.6524
MINDTL | 1.2818 | 0.9224
TRACER | 0.9617 | 0.7772
CBT | 0.9634 | 0.7886
Proposed | 0.9470 | 0.6381
Douban Book | MMMF | 0.7976 | 0.5414
MINDTL | 1.1970 | 0.9054
TRACER | 0.8047 | 0.6600
CBT | 0.7828 | 0.6130
Proposed | 0.7785 | 0.5022

Table 3: Values of Root Mean Square Error and Mean Absolute Error of baseline
methods and proposed method on MovieLens data by considering same dataset as
source and target

Dataset | Method | RMSE | MAE
---|---|---|---
MovieLens 100K | MMMF | 0.9828 | 0.6808
MINDTL | 1.9026 | 1.6538
TRACER | 1.0027 | 0.8213
CBT | 1.0264 | 0.8239
Proposed | 0.9653 | 0.6600
MovieLens 1M | MMMF | 0.9349 | 0.6389
MINDTL | 1.8545 | 1.5984
TRACER | 0.9892 | 0.8145
CBT | 1.0729 | 0.8725
Proposed | 0.9138 | 0.6175

We have conducted the experiments with varying number of clusters ($k_{1}$,
$k_{2}$ \- 25, 50, 75, 100, 125, 150, 175, 200) and Fig. 2 and Fig. 3 shows
the impact of number of clusters on RMSE and MAE of MovieLens 1M data when
MovieLens 100K data is considered as source. By observing the figures we can
say that the best performance for this scenario (source - MovieLens 100K,
target - MovieLens 1M) is achieved when the number of clusters are in between
$100$ and $150$ and hence we have fixed both $k_{1}$, $k_{2}$ to $125$. In the
similar way, when MovieLens 1M is used as source and Goodbooks is the target,
the number of clusters is fixed to $150$. Similarly, when Douban music data is
used as source and Douban book is the target, we set the number of clusters to
$100$.

Table 2 shows the values of RMSE and MAE of MovieLens 1M dataset, Goodbooks
dataset, Douban Book dataset. For MovieLens 1M data, MovieLens 100K is
considered as source and for Goodbooks data, MovieLens 1M is considered as
source, and for Douban Book data, Douban music is considered as source data.
The first column gives different datasets considered, whereas second column
gives the methods considered for comparison along with the proposed method.
The third and the fourth columns gives the RMSE and MAE values of the
considered methods on datasets. By observing the table, we can say that the
proposed method is performing well on any of the datasets considered.

We have also experimented our method by considering the same datasets as
source and target data. It is nothing but reconstructing the same matrix by
using the extracted knowledge. Table 3 gives the results on MovieLens 100k
data when the same is used as source, and also on MovieLens 1M data when
itself is used as source. By observing the metric values (RMSE, MAE) in the
table, we can claim that our method is outperforming all the comparision
methods on the datasets considered.

##  5 Conclusions and Future work

A novel method based on codebook based transfer learning for cross domain
recommendation has been proposed in this paper. Maximum margin matrix
factorisation technique is used to construct the codebook from the source
user-item rating matrix. By making use of the _hingeloss_ function in a novel
way, the constructed codebook is transferred to the target domain. The
experimental results show that the sparse target matrix is approximated well
by the proposed method. Rather than considering only ratings of items, in the
future it would be worthwhile to consider social tags as well as other types
of data. Applying transfer learning in other real world scenarios other than
recommender systems could be another direction in which this research could be
taken forward.

## References

  * (1) Aggarwal, C. C. Recommender Systems: The Textbook, 1st ed.  Springer Publishing Company, Incorporated, 2016.
  * (2) Bobadilla, J., Ortega, F., Hernando, A., and Gutiérrez, A. Recommender systems survey.  Knowledge-Based Systems 46 (2013), 109 – 132.
  * (3) Chung, R., Sundaram, D., and Srinivasan, A. Integrated personal recommender systems.  In Proceedings of the Ninth International Conference on Electronic Commerce (New York, NY, USA, 2007), ICEC ’07, ACM, pp. 65–74.
  * (4) Ding, C., Li, T., Peng, W., and Park, H. Orthogonal nonnegative matrix tri-factorizations for clustering.  In Proceedings of the 12th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (New York, NY, USA, 2006), KDD ’06, ACM, pp. 126–135.
  * (5) He, M., Zhang, J., Yang, P., and Yao, K. Robust transfer learning for cross-domain collaborative filtering using multiple rating patterns approximation.  In Proceedings of the Eleventh ACM International Conference on Web Search and Data Mining (New York, NY, USA, 2018), WSDM ’18, ACM, pp. 225–233.
  * (6) He, M., Zhang, J., and Zhang, J. Mindtl: Multiple incomplete domains transfer learning for information recommendation.  China Communications 14 (11 2017), 218–236.
  * (7) Ji, K., Sun, R., Li, X., and Shu, W. Improving matrix approximation for recommendation via a clustering-based reconstructive method.  Neurocomput. 173, P3 (Jan. 2016), 912–920.
  * (8) Koren, Y., Bell, R., and Volinsky, C. Matrix factorization techniques for recommender systems.  Computer 42, 8 (Aug. 2009), 30–37.
  * (9) Li, B., Yang, Q., and Xue, X. Can movies and books collaborate?: Cross-domain collaborative filtering for sparsity reduction.  In Proceedings of the 21st International Jont Conference on Artifical Intelligence (San Francisco, CA, USA, 2009), IJCAI’09, Morgan Kaufmann Publishers Inc., pp. 2052–2057.
  * (10) Masnadi-Shirazi, H., and Vasconcelos, N. On the design of loss functions for classification: Theory, robustness to outliers, and savageboost.  In Proceedings of the 21st International Conference on Neural Information Processing Systems (USA, 2008), NIPS’08, Curran Associates Inc., pp. 1049–1056.
  * (11) Moreno, O., Shapira, B., Rokach, L., and Shani, G. Talmud: Transfer learning for multiple domains.  In Proceedings of the 21st ACM International Conference on Information and Knowledge Management (New York, NY, USA, 2012), CIKM ’12, ACM, pp. 425–434.
  * (12) Pan, S. J., and Yang, Q. A survey on transfer learning.  IEEE Trans. on Knowl. and Data Eng. 22, 10 (Oct. 2010), 1345–1359.
  * (13) Pan, W. A survey of transfer learning for collaborative recommendation with auxiliary data.  Neurocomputing 177 (2016), 447 – 453.
  * (14) Pan, W., Liu, N. N., Xiang, E. W., and Yang, Q. Transfer learning to predict missing ratings via heterogeneous user feedbacks.  In Proceedings of the Twenty-Second International Joint Conference on Artificial Intelligence - Volume Three (2011), IJCAI’11, AAAI Press, pp. 2318–2323.
  * (15) Pan, W., and Ming, Z. Interaction-rich transfer learning for collaborative filtering with heterogeneous user feedback.  IEEE Intelligent Systems 29 (11 2014), 48–54.
  * (16) Pan, W., Nan Liu, N., Wei Xiang, E., and Yang, Q. Transfer learning to predict missing ratings via heterogeneous user feedbacks.  pp. 2318–2323.
  * (17) Pan, W., Xia, S., Liu, Z., Peng, X., and Ming, Z. Mixed factorization for collaborative recommendation with heterogeneous explicit feedbacks.  Information Sciences 332 (2016), 84 – 93.
  * (18) Pan, W., Xiang, E. W., Liu, N. N., and Yang, Q. Transfer learning in collaborative filtering for sparsity reduction.  In Proceedings of the Twenty-Fourth AAAI Conference on Artificial Intelligence (2010), AAAI’10, AAAI Press, pp. 230–235.
  * (19) Portugal, I., Alencar, P., and Cowan, D. The use of machine learning algorithms in recommender systems: A systematic review.  Expert Systems with Applications 97 (11 2015).
  * (20) Rennie, J. D. M., and Srebro, N. Fast maximum margin matrix factorization for collaborative prediction.  In Proceedings of the 22Nd International Conference on Machine Learning (New York, NY, USA, 2005), ICML ’05, ACM, pp. 713–719.
  * (21) Ricci, F., Rokach, L., Shapira, B., and Kantor, P. B. Recommender Systems Handbook, 1st ed.  Springer-Verlag, Berlin, Heidelberg, 2010.
  * (22) Rosasco, L., De Vito, E., Caponnetto, A., Piana, M., and Verri, A. Are loss functions all the same?  Neural Comput. 16, 5 (May 2004), 1063–1076.
  * (23) Shi, J., Long, M., Liu, Q., Ding, G., and Wang, J. Twin bridge transfer learning for sparse collaborative filtering.  In Advances in Knowledge Discovery and Data Mining (Berlin, Heidelberg, 2013), J. Pei, V. S. Tseng, L. Cao, H. Motoda, and G. Xu, Eds., Springer Berlin Heidelberg, pp. 496–507.
  * (24) Srebro, N., Rennie, J. D. M., and Jaakola, T. S. Maximum-margin matrix factorization.  In Advances in Neural Information Processing Systems 17 (2005), MIT Press, pp. 1329–1336.
  * (25) Veeramachaneni, S. D., Pujari, A. K., Padmanabhan, V., and Kumar, V. A maximum margin matrix factorization based transfer learning approach for cross-domain recommendation.  Appl. Soft Comput. 85 (2019).
  * (26) Zhao, L., Pan, S. J., Xiang, E. W., Zhong, E., Lu, Z., and Yang, Q. Active transfer learning for cross-system recommendation.  In Proceedings of the Twenty-Seventh AAAI Conference on Artificial Intelligence (2013), AAAI’13, AAAI Press, pp. 1205–1211.
  * (27) Zhao, L., Pan, S. J., and Yang, Q. A unified framework of active transfer learning for cross-system recommendation.  Artificial Intelligence 245 (2017), 38 – 55.
  * (28) Zhuang, F., Zheng, J., Chen, J., Zhang, X., Shi, C., and He, Q. Transfer collaborative filtering from multiple sources via consensus regularization.  Neural Networks 108 (2018), 287 – 295.