A fractional step scheme with modified characteristic finite differences run- ning 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 differ- ence operators, multiplicative commutation rule of difference operators, decomposition of high order difference operators, induction hypothesis, and prior estimates, an optimal order in 12 norm is displayed to complete the convergence analysis of the numerical algo- rithm. Some numerical results arising in the actual simulation of migration-accumulation of oil resources by this method are listed in the last section.