“…Moreover, we use the short forms for the value of the transport properties at intermediate grid points Therefore, we extend the scheme from [20] by incorporating the discrete convective and source terms [32]. The convective (first derivative) term is discretized by using an upwind approach [9,10,17,28,29,32], while the diffusion (second derivative) is approached by freezing coefficients [20], see also [4,8,12,13,18,19,26]. The Forchheimer term is discretized by freezing the non-positive coefficient at the previous time level t n , while the velocity factor is evaluated at the new time level t n+1 .…”