Discrete energy balance equation via a symplectic second-order method for two-phase flow in porous media

We propose and analyze a second-order partitioned time-stepping method for a two-phase flow problem in porous media. The algorithm is based on a refactorization of Cauchy's one-leg $\theta$-method. The main part consists of the implicit backward Euler method on $[t^n, t^{n+\theta}]$, while part two uses a linear extrapolation on $[t^{n+\theta},t^{n+1}]$ to obtain the solution at $t^{n+1}$, equivalent to the forward Euler method. In the backward Euler step, the decoupled equations are solved iteratively. We prove that the iterations converge linearly to the solution of the coupled problem, under some conditions on the data. When $\theta = 1/2$, the algorithm is equivalent to the symplectic midpoint method. In the absence of the chain rule for time-discrete setting, we approximate the change in the free energy by the product of a second-order accurate discrete gradient (chemical potential) and the one-step increment of the state variables. Similar to the continuous case, we also prove a discrete Gibbs free energy balance equation, without numerical dissipation. In the numerical tests we compare this implicit midpoint method with the classic backward Euler method, and two implicit-explicit time-lagging schemes. The midpoint method outperforms the other schemes in terms of rates of convergence, long-time behavior and energy approximation, for small and large values of the time step.