Abstract
The Buckley–Leverett equation is a fundamental model of two-phase flow in porous media, widely used in petroleum engineering to describe waterflooding. The analytical solution to this nonlinear hyperbolic problem contains discontinuous shock waves and requires Oleinik's entropy condition to select a physically admissible weak solution. This study investigates whether the Fourier Neural Operator (FNO) can learn the time-evolution map of a fixed Riemann problem solution for the Buckley–Leverett equation at mobility ratio M = 1. Unlike operator learning over a family of initial conditions, we consider surrogate modeling where the input is a fixed initial saturation profile together with a time instant , and the output is the saturation profile at that instant. The training dataset comprises 1000 temporal slices of a single analytical solution obtained by the method of characteristics with the Rankine–Hugoniot condition satisfied, for randomly selected time instants T ∈ [0.01, 1.0]. The trained FNO model demonstrated high approximation accuracy: relative L2 error of 0.0191+-0.0253, coefficient of determination R2 = 0.997 +- 0.0127, and mean absolute error = 0.0023 +- 0.0016. The shock front position was predicted with a displacement of less than one spatial grid step (4 x 10^-4) across all time quartiles. The main limitation was the Gibbs oscillations near the discontinuity, manifested as local oscillations. The results show that the FNO architecture can effectively approximate the structure of the shock-wave solution to a fixed Riemann problem, including the rarefaction zone and the saturation front position.