Generic differential equation for fractional flow of steady two-phase flow in porous media
We report on generic relations between fractional flow and pressure in steady two-phase flow in porous media. The main result is a differential equation for fractional flow as a function of phase saturation. We infer this result from two underlying observations of steady flow simulations in two and three dimensions using biperiodic boundary conditions. The resulting equation is solved generally, and the result is tested against simulations and experimental relative permeability results found in the literature.
main.tex
%% ****** Start of file template.aps ****** %
%%
%%
%%   This file is part of the APS files in the REVTeX 4 distribution.
%%   Version 4.0 of REVTeX, August 2001
%%
%%
%%   Copyright (c) 2001 The American Physical Society.
%%
%%   See the REVTeX 4 README file for restrictions and more information.
%%
%
% This is a template for producing manuscripts for use with REVTEX 4.0
% Copy this file to another name and then work on that file.
% That way, you always have this original template file to use.
%
% Group addresses by affiliation; use superscriptaddress for long
% author lists, or if there are many overlapping affiliations.
% For Phys. Rev. appearance, change preprint to twocolumn.
% Choose pra, prb, prc, prd, pre, prl, prstab, or rmp for journal
%  Add 'draft' option to mark overfull boxes with black boxes
%  Add 'showpacs' option to make PACS codes appear
%  Add 'showkeys' option to make keywords appear
%\documentclass[aps,pre,preprint,groupedaddress]{revtex4}
%\documentclass[aps,prl,preprint,superscriptaddress]{revtex4}
\documentclass[aps,prl,twocolumn,groupedaddress]{revtex4}

% You should use BibTeX and apsrev.bst for references
% Choosing a journal automatically selects the correct APS
% BibTeX style file (bst file), so only uncomment the line
% below if necessary.
%\bibliographystyle{apsrev}

\usepackage{graphicx}
\begin{document}

% Use the \preprint command to place your local institutional report
% number in the upper righthand corner of the title page in preprint mode.
% Multiple \preprint commands are allowed.
% Use the 'preprintnumbers' class option to override journal defaults
% to display numbers if necessary
%\preprint{}
%Title of paper
\title{Generic differential equation for fractional flow of steady two-phase flow in porous media}

% repeat the \author .. \affiliation  etc. as needed
% \email, \thanks, \homepage, \altaffiliation all apply to the current
% author. Explanatory text should go in the []'s, actual e-mail
% address or url should go in the {}'s for \email and \homepage.
% Please use the appropriate macro foreach each type of information

% \affiliation command applies to all authors since the last
% \affiliation command. The \affiliation command should follow the
% other information
% \affiliation can be followed by \email, \homepage, \thanks as well.
\author{Henning Arendt Knudsen}
\email[]{arendt@phys.ntnu.no}
%\homepage[]{Your web page}
%\thanks{}%\altaffiliation{}
\noaffiliation{}
\author{Alex Hansen}
\email[]{Alex.Hansen@phys.ntnu.no}
\noaffiliation{}

%Collaboration name if desired (requires use of superscriptaddress
%option in \documentclass). \noaffiliation is required (may also be
%used with the \author command).
%\collaboration can be followed by \email, \homepage, \thanks as well.
%\collaboration{}
%\noaffiliation

\date{\today}


\begin{abstract}
We report on generic relations between fractional flow and pressure in steady two-phase flow in porous media. The main result is a differential equation for fractional flow as a function of phase saturation. We infer this result from two underlying observations of steady flow simulations in two and three dimensions using biperiodic boundary conditions. The resulting equation is solved generally, and the result is tested against simulations and experimental relative permeability results found in the literature.
\end{abstract}

% insert suggested PACS numbers in braces on next line
\pacs{47.55.Mh}
% insert suggested keywords - APS authors don't need to do this
%\keywords{}

%\maketitle must follow title, authors, abstract, \pacs, and \keywords
\maketitle


% body of paper here - Use proper section commands
% References should be done using the \cite, \ref, and \label commands
% Put \label in argument of \section for cross-referencing
%\section{\label{}}
%\subsection{}
%\subsubsection{}

Two-phase flow in porous media is possible to classify as either steady or unsteady flow. While steady flow is the one that is statistically invariant in time, the typical unsteady flow is the displacement of one phase by another. In particular, many authors have studied the motion of the displacement front. This was done by experiment\cite{LTZ88}, network modeling\cite{LTZ88}, lattice Boltzmann methods\cite{R90} and also statistical models: diffusion limited aggregation\cite{WS81} and invasion percolation\cite{WW83}. The study of steady flow is often associated with finding the relative permeability curves, although this concept can be defined and is used both for steady and unsteady flow\cite{D92,TD96}. Comparisons between relative permeability obtained by the two methods exist in the literature\cite{PK87}. However, in the complexity of two-phase flow in porous media, pure displacements and true steady flow can be considered to be two limiting situations. Thus, by nature the two situations are very different, and one should expect differences in the relative permeability curves. Nevertheless, unsteady methods are generally preferred because steady methods take long time and are more expensive. The results that we present here concern steady-state flow. A lattice-Boltzmann approach to steady-state flow is found in\cite{LP01}. Further, a very interesting study using a network approach is found in\cite{VP01}.

The body of knowledge in the field is immense\cite{D92,S95}. Pore-network modeling is popular and a good overview of the state-of-art is given in\cite{B01}. Basic mechanisms of displacement on pore level are largely known, and much qualitative knowledge exists about properties on larger scale. However, there is a good deal of work left to do in order to bridge the gap between these two worlds. Numerical simulations is a unique tool to that respect being very well controllable and more precise. In particular, we find that while experimental data points often are few and scattered, we have the possibility to generate sufficient data points to obtain good statistics. Further, in simulation series we can control the \emph{ensemble} completely. Commonly used ensembles for relative permeability or fractional flow curves as functions of saturation, are constant global pressure, constant total flux or constant capillary number(Ca). We argue that the choice of ensemble can be of major importance and is often underestimated. Our results are presented in the formalism of fractional flow and globally applied pressure. We find that when using the constant Ca ensemble, there exist a generic differential equation for the fractional flow. In turn a general solution is found.

The results of the paper is based on a network simulator for immiscible two-phase flow. This line of modeling which is based on Washburn's equation\cite{W21} dates back to the work of several groups in the mid 1980's\cite{KL85,CP91,CP96}. Our model is a continuation of the model developed by Aker \emph{et al.}\cite{AMHB98,AMH98}. Although a thorough presentation can be found in \cite{KAH02}, we provide a brief r{\'e}sum{\'e} of main aspects for clarity.

The porous media are represented by networks of tubes, forming square lattices in two dimensions (2D) and cubic lattices in 3D, both tilted with respect to the imposed pressure gradient and thus the overall direction of flow. We refer to the the lattice points were four(2D) and six(3D) tubes meet as nodes. Volume in the model is contained in the tubes and not in the nodes, although effective node volumes are used in the modeling of transport through the nodes. For further details, see\cite{KAH02}. Randomness is incorporated by distorting the nodes on random within a circle(2D) and a sphere(3D) around their respective lattice positions. This gives a distribution of tube lengths in the system. Further, the radii are drawn from a flat distribution so that the radius of a given tube is $r\in(0.1l,0.4l)\ $where $l\ $�is the length of that tube.

The model is filled with two phases that flow within the system of tubes. The flow in each tube obeys the Washburn\cite{W21} equation, $q=-(\sigma k/\mu)(\Delta p-\sum{p_c})/l$. With respect to momentum transfer these tubes are cylindrical with cross-sectional area $\sigma\ $�and length $l$. The permeability is $k=r^2/8\ $which is known for Hagen-Poiseulle flow. Further, $\mu\ $�is the viscosity of the phase present in the tube. If both phases are present the volume average of their viscosities is used. The volumetric flow rate is denoted by $q\ $and the pressure difference between the ends of the tube by $\Delta p$. The summation is the sum over all capillary pressures $p_c\ $�within the tube. With respect to capillary pressure the tubes are hour-glass shaped meaning that a meniscus at position $x\in(0,l)\ $in the tube has capillary pressure: $p_c=(2\gamma/r)[1-\cos{(2\pi x/l)}]$, where $\gamma\ $is the interfacial tension between the two phases. This is a modified version of the Young-Laplace law\cite{D92,AMHB98}.

Biperiodic(2D) and triperiodic(3D) boundary conditions are used so that the flow is by construction steady flow. That is to say, the systems are closed so both phases retain their initial volume fractions, i.e. their saturations. These boundary conditions are in 2D equivalent to saying that the flow is restrained to be on the surface of a torus. The flow is driven by a globally applied pressure gradient. Usually invasion processes are driven by setting up a pressure fall between two borders, inlet and outlet. Since, by construction, the outlet is directly joined with the inlet in our system, we give instead a so-called global pressure drop when passing this line or cut through the system\cite{roux,KAH02}. Effectively we put restraints on the pressure gradient that is experienced throughout the network. Integration of the pressure gradient along an arbitrary closed path one lap around the system, making sure to pass the `inlet-outlet'-cut once, should add up to the same global pressure fall.

For a given value of the global pressure fall, and given the instantaneous distribution of the phases in the system, the pressure field is calculated numerically, and the flow distribution at that moment is known. Given the flow field the whole system is updated according to the Euler scheme. All interfaces, menisci, which are the essential entities keeping track of the distribution of the phases, are moved according to the fluid velocity in the tube where they are situated. After having updated the hole system, the new distribution of the phases might give new effective viscosities in the tubes and different capillary pressures across the menisci, leading to a recalculation of the coefficients in the set of equations for the pressure. The tricky part is how to model the transport of the menisci across the nodes, and this is done by applying a set of rules that assures volume conservation of both phases and that are acceptably close to the real world. We refer to the previous study for details on this\cite{KAH02}.

All simulations being the basis for our conclusions are done in series where the saturation is varied as the independent variable and the capillary number as well as the viscosity ratio $M=\mu_{\rm nw}/\mu_{\rm w}\ $is held fixed. We define the capillary number as ${\rm Ca}=Q_{\rm tot}\mu_{\rm eff}/(\Sigma\gamma)\ $�where $\Sigma\ $�is the cross-sectional area of the entire network, $Q_{\rm tot}\ $is the total flux through the network and $\mu_{\rm eff}=\mu_{\rm nw}F_{\rm nw}+\mu_{\rm w}F_{\rm w}$. Here, $F_{\rm nw}\ $�and $F_{\rm w}\ $are nonwetting and wetting fractional flow respectively. In invasion studies one often employs either the wetting or the nonwetting viscosity and refer to the wetting or nonwetting capillary number, respectively. However, in steady flow it is appropriate to take the volume average. The resulting dependent variables are the fractional flow of the two phases and the globally applied pressure. These two variables are in a one-to-one relationship with the two relative permeabilities, but we will stick to the former notation here, for clarity and convenience.

From the simulations two observations are made from which we infer the following differential equation for the fractional flow as a function of saturation;
\begin{equation}
\label{eq:diff}
\frac{d^2F}{dS^2}-(2aF+b)\frac{dF}{dS}=0.
\end{equation}
The saturation is the nonwetting saturation $S=S_{\rm nw}\ $and the fractional flow is the nonwetting fractional flow $F=F_{\rm nw}(S_{\rm nw})\ $here and in illustrations unless where explicitly stated otherwise. However, the relation is equivalent in form for the wetting counterpart $F_{\rm w}(S_{\rm w})$, only the values of $a\ $and $b\ $change.

The two observations on which Eq. (\ref{eq:diff}) is based are of a phenomenological character. The physical problem is intractable from an analytical point of view. By means of simulations the evolution and properties are monitored in detail. The first observation is that the derivative of the fractional flow is related to global pressure by
\begin{equation}
\label{eq:dfds}
P(S)=A\frac{dF(S)}{dS}+B,
\end{equation}
where $P\ $�denotes normalized pressure, i.e., the actual pressure for the fixed flux at a given saturation divided by the pressure for single phase flow at the same flux. This result has been reported on in great detail for the case of viscosity matching phases\cite{KH02}. Those results were solely based on simulations in two dimensions. However, the result is extendable to 3D, see Fig.\ref{fig:dfdsfit}, as well as viscosity contrasts.

\begin{figure}
\includegraphics[width=8cm]{dfdsfit.eps}
\caption{\label{fig:dfdsfit}
The figure shows that the relation between pressure and the derivative of the fractional is of the form in Eq. (\ref{eq:dfds}). The curves are for viscosity matching phases for four different capillary numbers in 2D and two capillary numbers in 3D. The nonwetting saturation $S_{\rm nw}\ $serves as parameter for the curves, and the intervals in which the data is plotted are as follows; (2D) ${\rm Ca}=3.2\times10^{-2}:S_{\rm nw}\in(0.14,0.93)$,${\rm Ca}=1.0\times10^{-2}:S_{\rm nw}\in(0.14,0.83)$,${\rm Ca}=3.2\times10^{-3}:S_{\rm nw}\in(0.14,0.83)$,${\rm Ca}=1.0\times10^{-3}:S_{\rm nw}\in(0.20,0.73)$, (3D) ${\rm Ca}=1.0\times10^{-2}:S_{\rm nw}\in(0.10,0.90)$,${\rm Ca}=1.0\times10^{-3}:S_{\rm nw}\in(0.20,0.75)$. The two curves in 3D are lowered by one unit for clarity.
}
\end{figure}

\begin{figure}
\begin{tabular}{ll}
\includegraphics[width=4cm]{pffit3.eps}
&
\includegraphics[width=4cm]{pffit4.eps}
\\
(a) ${\rm Ca}=3.2\times10^{-3}$ &(b) $ {\rm Ca}=1.0\times10^{-3}$
\\
\includegraphics[width=4cm]{pffit5.eps}
&
\includegraphics[width=4cm]{pffit6.eps}
\\
(c) ${\rm Ca}=3.2\times10^{-4}$ &(d) $ {\rm Ca}=1.0\times10^{-4}$
\end{tabular}
\caption{\label{fig:pffit}
The figure shows fitted curves to $P(F_{\rm nw})\ $according to Eq. (\ref{eq:qvadr}). The results are from the same series of simulations as the curves in 2D in Fig. \ref{fig:dfdsfit}, only showing results for two lower capillary numbers. The shown region of $F_{\rm nw}\ $is the one where the fit is good. For these capillary number this region is large, for the two higher values of ${\rm Ca}\ $�shown if Fig. \ref{fig:dfdsfit}, the region of good fit is; ${\rm Ca}=3.2\times10^{-2}:F_{\rm nw}\in(0.15,0.45)$,${\rm Ca}=1.0\times10^{-2}:F_{\rm nw}\in(0.10,0.50)$. The data is from five realizations of the porous network giving slightly shifted curves as we can see. The scattering of the data becomes larger for smaller capillary numbers due to increased hysteresis and history effects.
}
\end{figure}

The second observation is that for a broad range of saturation values the pressure forms a quadratic function of the fractional flow. In a mathematical formulation this result becomes
\begin{equation}
\label{eq:qvadr}
P(F)=a'F^2+b'F+c'.
\end{equation}
The range of validity of this equation is illustrated by the quadratic fits in Figs. \ref{fig:pffit} and \ref{fig:pfLfit}. In Fig. \ref{fig:pffit} the viscosity ratio is unity, and we observe how Eq. (\ref{eq:qvadr}) holds for each fixed value of ${\rm Ca}$. Likewise, we show in Fig. \ref{fig:pfLfit} quadratic fits to $P(F)$-curves for one given capillary number ${\rm Ca}=3.2\times10^{-3}$, varying the viscosity ratio from $M=30\ $to $M=1/30$. The point is that for most of the range that we are interested in, it is sufficient to expand $P(F)\ $to second order in $F\ $around the maximum point. Alternatively, an expansion of $P(S)\ $in $S\ $�turns out to require higher order terms, and is hence less useful. Other combinations of capillary number and viscosity ratio have been performed and they are found to agree with this result. The shifting of the curves and detailed dependence on ${\rm Ca}\ $and $M\ $�are to be discussed elsewhere\cite{KHunpub02}.

\begin{figure}
\includegraphics[width=8cm]{pfLfit.eps}
\caption{\label{fig:pfLfit}
We show $P(F_{\rm nw})\ $for a fixed value of the capillary number, namely ${\rm Ca}=3.2\times10^{-3}$, which corresponds to Fig. \ref{fig:pffit}(a). We study one sample assigning five different viscosity ratios $M\ $to the phases. We observe how it is possible to fit all the curves to the form in Eq. (\ref{eq:qvadr}) for $M\in(1/30,30)$.
}
\end{figure}


\begin{figure}
\includegraphics[width=8cm]{fsLg3fit.eps}
\caption{\label{fig:fsLg3fit}
The figure shows how four samples of fractional flow curves from our simulations can be fitted to the general solution of Eq. (\ref{eq:diff}) as it is given in Eq. (\ref{eq:solution}). The dimensionality, viscosity ratio and capillary number of these samples are as indicated in the legend. For clarity the curves are shifted so that the two curves in 2D belong the lower x-axis, and the two curves in 3D belong to the upper x-axis. The 3D curves are from the same simulation series as the 3D curves in Fig.\ref{fig:dfdsfit}.
}
\end{figure}

\begin{figure}[!h]
\includegraphics[width=8cm]{PKandSYfit.eps}
\caption{\label{fig:PKandSYfit}
Data from the literature: Sharma and Yen(SY)\cite{SY83}, Peters and Khataniar(PK)\cite{PK87},Braun and Blackwell(BB)\cite{BB81,D92}. The water saturation is also the wetting saturation $S_{\rm w}$, but note that it is partly the different normalization of this entity, that was employed in the respective papers, which decide the placement of the curves. The SY-curves are typical samples of the set of functional forms that are employed when discussing how various physical parameters influence on fractional flow. Regarding PK and BB-curves, see Fig.\ref{fig:PKdfP}. All the curves are fitted to the solution of Eq.(\ref{eq:diff}) as it is given in Eq.(\ref{eq:solution}).
}
\end{figure}

\begin{figure}[!h]
\includegraphics[width=8cm]{PKdfP.eps}
\caption{\label{fig:PKdfP}
The curves correspond to the respective curves of Fig.\ref{fig:PKandSYfit}. For PK the reported value of the viscosity ratio $M=107.2\ $is used for the calculation from relative permeabilities to fractional flow and pressure. However, for BB we have chosen two different reasonable values for the viscosity ratio, to illuminate its effect. 
}
\end{figure}

Recall that $S\ $�is the independent variable upon which both $F\ $�and $P\ $depend. The two Eqs.(\ref{eq:dfds}) and (\ref{eq:qvadr}) for pressure $P\ $must be equal for all saturations $S$. Combining the two gives a first order differential equation that differentiated once again with respect to $S\ $�becomes Eq.(\ref{eq:diff}). Note that in Eq. (\ref{eq:diff}) $a=a'/A\ $�and $b=b'/A$. The second order version being more elegant is, however, autonomous and just as simple to solve using standard methods\cite{BO78}. The general solution is
\begin{equation}
\label{eq:solution}
F(S)=-\frac{1}{2a} \frac{(b-k)+(b+k)\exp{[k(S-S_0)]}}{1+\exp{[k(S-S_0)]}},
\end{equation}
where
\begin{equation}
k^2=b^2-4ak_0 ,
\end{equation}
and $k_0\ $and $S_0\ $are constants of integration. This result is generic and general. In Fig.\ref{fig:fsLg3fit} we provide the fractional flow by simulations for four different sets of parameters. The sets are fitted by the form in Eq.(\ref{eq:solution}). We observe that the fits are excellent. Even though the two underlying observations are valid only in a certain central region (say 70-90\%) of the curves, the solution is satisfactory in the entire two-phase region.

Relative permeability is used both for steady-state flow and unsteady flow. The literature covers several methods to produce relative permeability from experimental data\cite{D92,TD96}. It is a priori not entirely clear to which extent a given study is comparable to our simulations. We have chosen a few samples from literature which are provided in Figs.\ref{fig:PKandSYfit} and \ref{fig:PKdfP}. In literature, when focus is on how some physical or chemical properties influence on fractional flow, this is often done by employing a class of functional forms for the fractional flow. The curves indicated by SY are samples of that\cite{SY83}. It is reassuring that these curves can be fitted by Eq.\ref{eq:solution}. The other three curves are calculated from the relative permeabilities $k_{\rm rnw}\ $and $k_{\rm rw}\ $and the viscosity ratio $M\ $using the formalism $P = (k_{\rm rnw}+Mk_{\rm rw})^{-1}\ $and $F_{\rm nw}�= Mk_{\rm rw} P$, which arise when assuming that the total flux is held constant under the series and $P\ $is normalized pressure, for details see\cite{D92,KH02}.

The PK-curve is a steady-state curve which is comparable to our simulations, even though the ensemble is constant total flux. The validity of the Eq.(\ref{eq:solution}) follows from Fig.\ref{fig:PKandSYfit}, likewise Eq.(\ref{eq:dfds}) follows from Fig.\ref{fig:PKdfP}. Here we use the actual viscosity ratio that was used in the experiments for our calculations. It is well established that relative permeability is a function of a large number of parameters including the viscosity ratio\cite{VP01}. However, to some extent one gets the impression that the curves are regarded as rock properties and are valid for at least a range of viscosity ratios. This might be very system dependent. The curves marked BB\cite{BB81,D92} in Figs.\ref{fig:PKandSYfit} and \ref{fig:PKdfP} are generated from relative permeability data in this way by choosing two values for viscosity ratio: $M=20\ $and $M=40$. Within this range the results are robust.

In general fractional flow curves are obtained for constant total flux. Another possible ensemble is constant pressure drop over a sample. In our simulations we have chosen a third ensemble, namely constant capillary number. By this method we keep the ratio between viscous and capillary forces fixed. The precision and control of parameters we obtain in the curves exceed what is normal in experiments. This is highly advantageous when looking for general relationships between variables. We believe that constant Ca is the appropriate ensemble for the study of steady flow properties. The comparisons with literature in this study indicate the correctness of the result, however, carefully designed laboratory experiments should be made using this ensemble to check our results.

In conclusion it is our claim that this general result is valid for all steady-state two-phase flow system that can be modeled by our numerical model. That is to say immiscible flow where film flow can be neglected. This is more likely to be case at moderate to high capillary numbers. Studies of imbibition and drainage by small scale experiments show that in particular imbibition like steps may be dominated by film flow at low capillary numbers, but not at high capillary numbers\cite{CK85}.


% If in two-column mode, this environment will change to single-column
% format so that long equations can be displayed. Use
% sparingly.
%\begin{widetext}
% put long equation here
%\end{widetext}

% figures should be put into the text as floats.
% Use the graphics or graphicx packages (distributed with LaTeX2e)
% and the \includegraphics macro defined in those packages.
% See the LaTeX Graphics Companion by Michel Goosens, Sebastian Rahtz,
% and Frank Mittelbach for instance.
%
% Here is an example of the general form of a figure:
% Fill in the caption in the braces of the \caption{} command. Put the label
% that you will use with \ref{} command in the braces of the \label{} command.
% Use the figure* environment if the figure should span across the
% entire page. There is no need to do explicit centering.

% \begin{figure}
% \includegraphics{}%
% \caption{\label{}}
% \end{figure}

% Surround figure environment with turnpage environment for landscape
% figure
% \begin{turnpage}
% \begin{figure}
% \includegraphics{}%
% \caption{\label{}}
% \end{figure}
% \end{turnpage}

% tables should appear as floats within the text
%
% Here is an example of the general form of a table:
% Fill in the caption in the braces of the \caption{} command. Put the label
% that you will use with \ref{} command in the braces of the \label{} command.
% Insert the column specifiers (l, r, c, d, etc.) in the empty braces of the
% \begin{tabular}{} command.
% The ruledtabular enviroment adds doubled rules to table and sets a
% reasonable default table settings.
% Use the table* environment to get a full-width table in two-column
% Add \usepackage{longtable} and the longtable (or longtable*}
% environment for nicely formatted long tables. Or use the the [H]
% placement option to break a long table (with less control than 
% in longtable).
% \begin{table}%[H] add [H] placement to break table across pages
% \caption{\label{}}
% \begin{ruledtabular}
% \begin{tabular}{}
% Lines of table here ending with \\
% \end{tabular}
% \end{ruledtabular}
% \end{table}

% Surround table environment with turnpage environment for landscape
% table
% \begin{turnpage}
% \begin{table}
% \caption{\label{}}
% \begin{ruledtabular}
% \begin{tabular}{}
% \end{tabular}
% \end{ruledtabular}
% \end{table}
% \end{turnpage}

% Specify following sections are appendices. Use \appendix* if there
% only one appendix.
%\appendix
%\section{}

% If you have acknowledgments, this puts in the proper section head.
\begin{acknowledgments}
H.A.K. acknowledges support from VISTA, a collaboration between Statoil and The Norwegian Acad. of Science and Letters. This work has received support from NTNU through a grant of computing time on the good supercomputing facilities at NTNU.
\end{acknowledgments}

% Create the reference section using BibTeX:
\bibliography{referanser}

\end{document}
%
% ****** End of file template.aps ******

main.bbl
\begin{thebibliography}{25}
\expandafter\ifx\csname natexlab\endcsname\relax\def\natexlab#1{#1}\fi
\expandafter\ifx\csname bibnamefont\endcsname\relax
  \def\bibnamefont#1{#1}\fi
\expandafter\ifx\csname bibfnamefont\endcsname\relax
  \def\bibfnamefont#1{#1}\fi
\expandafter\ifx\csname citenamefont\endcsname\relax
  \def\citenamefont#1{#1}\fi
\expandafter\ifx\csname url\endcsname\relax
  \def\url#1{\texttt{#1}}\fi
\expandafter\ifx\csname urlprefix\endcsname\relax\def\urlprefix{URL }\fi
\providecommand{\bibinfo}[2]{#2}
\providecommand{\eprint}[2][]{\url{#2}}

\bibitem[{\citenamefont{Lenormand et~al.}(1988)\citenamefont{Lenormand,
  Touboul, and Zarcone}}]{LTZ88}
\bibinfo{author}{\bibfnamefont{R.}~\bibnamefont{Lenormand}},
  \bibinfo{author}{\bibfnamefont{E.}~\bibnamefont{Touboul}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{C.}~\bibnamefont{Zarcone}},
  \bibinfo{journal}{J. Fluid Mech.} \textbf{\bibinfo{volume}{189}},
  \bibinfo{pages}{165} (\bibinfo{year}{1988}).

\bibitem[{\citenamefont{Rothman}(1990)}]{R90}
\bibinfo{author}{\bibfnamefont{D.~H.} \bibnamefont{Rothman}},
  \bibinfo{journal}{J Geophys. Res.} \textbf{\bibinfo{volume}{95}},
  \bibinfo{pages}{8663} (\bibinfo{year}{1990}).

\bibitem[{\citenamefont{Witten and Sander}(1981)}]{WS81}
\bibinfo{author}{\bibfnamefont{T.~A.} \bibnamefont{Witten}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{L.~M.} \bibnamefont{Sander}},
  \bibinfo{journal}{Phys. Rev. Lett.} \textbf{\bibinfo{volume}{47}},
  \bibinfo{pages}{1400} (\bibinfo{year}{1981}).

\bibitem[{\citenamefont{Wilkinson and Willemsen}(1983)}]{WW83}
\bibinfo{author}{\bibfnamefont{D.}~\bibnamefont{Wilkinson}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{J.~F.} \bibnamefont{Willemsen}},
  \bibinfo{journal}{J. Phys. A} \textbf{\bibinfo{volume}{16}},
  \bibinfo{pages}{3365} (\bibinfo{year}{1983}).

\bibitem[{\citenamefont{Dullien}(1992)}]{D92}
\bibinfo{author}{\bibfnamefont{F.~A.~L.} \bibnamefont{Dullien}},
  \emph{\bibinfo{title}{Porous Media: Fluid Transport and Pore Structure}}
  (\bibinfo{publisher}{Academic Press, San Diego}, \bibinfo{year}{1992}).

\bibitem[{\citenamefont{Tiab and Donaldson}(1996)}]{TD96}
\bibinfo{author}{\bibfnamefont{D.}~\bibnamefont{Tiab}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{E.~C.} \bibnamefont{Donaldson}},
  \emph{\bibinfo{title}{Petrophysics}} (\bibinfo{publisher}{Gulf Publishing
  Company, Houston}, \bibinfo{year}{1996}).

\bibitem[{\citenamefont{Peters and Khataniar}(1987)}]{PK87}
\bibinfo{author}{\bibfnamefont{E.~J.} \bibnamefont{Peters}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Khataniar}},
  \bibinfo{journal}{SPE Formation Evaluation} pp. \bibinfo{pages}{469--474}
  (\bibinfo{year}{1987}).

\bibitem[{\citenamefont{Langaas and Papatzacos}(2001)}]{LP01}
\bibinfo{author}{\bibfnamefont{K.}~\bibnamefont{Langaas}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{P.}~\bibnamefont{Papatzacos}},
  \bibinfo{journal}{Transport in Porous Media} \textbf{\bibinfo{volume}{45}},
  \bibinfo{pages}{241} (\bibinfo{year}{2001}).

\bibitem[{\citenamefont{Valavanides and Payatakes}(2001)}]{VP01}
\bibinfo{author}{\bibfnamefont{M.~S.} \bibnamefont{Valavanides}}
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{A.~C.}
  \bibnamefont{Payatakes}}, \bibinfo{journal}{Advances in Water Resources}
  \textbf{\bibinfo{volume}{24}}, \bibinfo{pages}{385} (\bibinfo{year}{2001}).

\bibitem[{\citenamefont{Sahimi}(1995)}]{S95}
\bibinfo{author}{\bibfnamefont{M.}~\bibnamefont{Sahimi}},
  \emph{\bibinfo{title}{Flow and Transport in Porous Media and Fractured Rock}}
  (\bibinfo{publisher}{VCH Verlagsgesellschaft mbH, Weinheim},
  \bibinfo{year}{1995}).

\bibitem[{\citenamefont{Blunt}(2001)}]{B01}
\bibinfo{author}{\bibfnamefont{M.~J.} \bibnamefont{Blunt}},
  \bibinfo{journal}{Current Opinion in Colloid \& Interface Science}
  \textbf{\bibinfo{volume}{6}}, \bibinfo{pages}{197} (\bibinfo{year}{2001}).

\bibitem[{\citenamefont{Washburn}(1921)}]{W21}
\bibinfo{author}{\bibfnamefont{E.~W.} \bibnamefont{Washburn}},
  \bibinfo{journal}{Phys. Rev.} \textbf{\bibinfo{volume}{17}},
  \bibinfo{pages}{273} (\bibinfo{year}{1921}).

\bibitem[{\citenamefont{Constantinides and Payatakes}(1991)}]{CP91}
\bibinfo{author}{\bibfnamefont{G.~N.} \bibnamefont{Constantinides}}
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{A.~C.}
  \bibnamefont{Payatakes}}, \bibinfo{journal}{J. Colloid Interface Sci.}
  \textbf{\bibinfo{volume}{141}}, \bibinfo{pages}{486} (\bibinfo{year}{1991}).

\bibitem[{\citenamefont{Constantinides and Payatakes}(1996)}]{CP96}
\bibinfo{author}{\bibfnamefont{G.~N.} \bibnamefont{Constantinides}}
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{A.~C.}
  \bibnamefont{Payatakes}}, \bibinfo{journal}{AIChE Journal}
  \textbf{\bibinfo{volume}{42}}, \bibinfo{pages}{369} (\bibinfo{year}{1996}).

\bibitem[{\citenamefont{Koplik and Lasseter}(1985)}]{KL85}
\bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Koplik}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{T.~J.} \bibnamefont{Lasseter}},
  \bibinfo{journal}{SPE J} \textbf{\bibinfo{volume}{25}}, \bibinfo{pages}{89}
  (\bibinfo{year}{1985}).

\bibitem[{\citenamefont{Aker et~al.}(1998{\natexlab{a}})\citenamefont{Aker,
  M{\aa}l{\o}y, and Hansen}}]{AMH98}
\bibinfo{author}{\bibfnamefont{E.}~\bibnamefont{Aker}},
  \bibinfo{author}{\bibfnamefont{K.~J.} \bibnamefont{M{\aa}l{\o}y}},
  \bibnamefont{and} \bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Hansen}},
  \bibinfo{journal}{Phys. Rev. E} \textbf{\bibinfo{volume}{58}},
  \bibinfo{pages}{2217} (\bibinfo{year}{1998}{\natexlab{a}}).

\bibitem[{\citenamefont{Aker et~al.}(1998{\natexlab{b}})\citenamefont{Aker,
  M{\aa}l{\o}y, Hansen, and Batrouni}}]{AMHB98}
\bibinfo{author}{\bibfnamefont{E.}~\bibnamefont{Aker}},
  \bibinfo{author}{\bibfnamefont{K.~J.} \bibnamefont{M{\aa}l{\o}y}},
  \bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Hansen}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{G.~G.} \bibnamefont{Batrouni}},
  \bibinfo{journal}{Transport in Porous Media} \textbf{\bibinfo{volume}{32}},
  \bibinfo{pages}{163} (\bibinfo{year}{1998}{\natexlab{b}}).

\bibitem[{\citenamefont{Knudsen et~al.}(2002)\citenamefont{Knudsen, Aker, and
  Hansen}}]{KAH02}
\bibinfo{author}{\bibfnamefont{H.~A.} \bibnamefont{Knudsen}},
  \bibinfo{author}{\bibfnamefont{E.}~\bibnamefont{Aker}}, \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Hansen}},
  \bibinfo{journal}{Transport in Porous Media} \textbf{\bibinfo{volume}{47}},
  \bibinfo{pages}{99} (\bibinfo{year}{2002}).

\bibitem[{\citenamefont{Roux}(unpublished)}]{roux}
\bibinfo{author}{\bibfnamefont{S.}~\bibnamefont{Roux}}
  (\bibinfo{year}{unpublished}).

\bibitem[{\citenamefont{Knudsen and Hansen}(2002{\natexlab{a}})}]{KH02}
\bibinfo{author}{\bibfnamefont{H.~A.} \bibnamefont{Knudsen}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Hansen}},
  \bibinfo{journal}{Phys. Rev. E}  (\bibinfo{year}{2002}{\natexlab{a}}).

\bibitem[{\citenamefont{Knudsen and Hansen}(2002{\natexlab{b}})}]{KHunpub02}
\bibinfo{author}{\bibfnamefont{H.~A.} \bibnamefont{Knudsen}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{A.}~\bibnamefont{Hansen}},
  \bibinfo{journal}{in preparation for Phys. Rev. E.}
  (\bibinfo{year}{2002}{\natexlab{b}}).

\bibitem[{\citenamefont{Sharma and Yen}(1983)}]{SY83}
\bibinfo{author}{\bibfnamefont{M.~M.} \bibnamefont{Sharma}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{T.~F.} \bibnamefont{Yen}},
  \bibinfo{journal}{SPE Journal} pp. \bibinfo{pages}{125--134}
  (\bibinfo{year}{1983}).

\bibitem[{\citenamefont{Braun and Blackwell}(1981)}]{BB81}
\bibinfo{author}{\bibfnamefont{E.~M.} \bibnamefont{Braun}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{R.~J.} \bibnamefont{Blackwell}},
  \bibinfo{journal}{SPE-AIME 56th Conference, Texas}  (\bibinfo{year}{1981}).

\bibitem[{\citenamefont{Bender and Orszag}(1978)}]{BO78}
\bibinfo{author}{\bibfnamefont{C.~M.} \bibnamefont{Bender}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{S.~A.} \bibnamefont{Orszag}},
  \emph{\bibinfo{title}{Advanced mathematical methods for scientists and
  engineers}} (\bibinfo{publisher}{McGraw-Hill}, \bibinfo{year}{1978}).

\bibitem[{\citenamefont{Chen and Koplik}(1985)}]{CK85}
\bibinfo{author}{\bibfnamefont{J.-D.} \bibnamefont{Chen}} \bibnamefont{and}
  \bibinfo{author}{\bibfnamefont{J.}~\bibnamefont{Koplik}},
  \bibinfo{journal}{J. Colloid Interface Sci.} \textbf{\bibinfo{volume}{108}},
  \bibinfo{pages}{304} (\bibinfo{year}{1985}).

\end{thebibliography}
