A fractional step scheme with modified characteristic finite differences running in a parallel arithmetic is presented to simulate a nonlinear percolation system of multilayer dynamics of fluids in a porous medium with moving boundary values. With the help of theoretical techniques including the change of regions, piecewise threefold quadratic interpolation, calculus of variations, multiplicative commutation rule of difference operators, multiplicative commutation rule of difference operators, decomposition of high order difference operators, induction hypothesis, and prior estimates, an optimal order in l² norm is displayed to complete the convergence analysis of the numerical algorithm. Some numerical results arising in the actual simulation of migration-accumulation of oil resources by this method are listed in the last section.