  We present a Markov Chain Monte Carlo algorithm based on the Metropolis
algorithm for simulation of the flow of two immiscible fluids in a porous
medium under macroscopic steady-state conditions using a dynamical pore network
model that tracks the motion of the fluid interfaces. The Monte Carlo algorithm
is based on the configuration probability, where a configuration is defined by
the positions of all fluid interfaces. We show that the configuration
probability is proportional to the inverse of the flow rate. Using a
two-dimensional network, advancing the interfaces using time integration the
computational time scales as the linear system size to the fourth power,
whereas the Monte Carlo computational time scales as the linear size to the
second power. We discuss the strengths and the weaknesses of the algorithm.
