A moving-boundary model of reactive settling in wastewater treatment. Part 2: Numerical scheme

A numerical scheme is proposed for the simulation of reactive settling in sequencing batch reactors (SBRs) in wastewater treatment plants. Reactive settling is the process of sedimentation of flocculated particles (biomass; activated sludge) consisting of several material components that react with substrates dissolved in the fluid. An SBR is operated in cycles of consecutive fill, react, settle, draw and idle stages, which means that the volume in the tank varies and the surface moves with time. The process is modelled by a system of spatially one-dimensional, nonlinear, strongly degenerate parabolic convection-diffusion-reaction equations. This system is coupled via conditions of mass conservation to transport equations on a half line whose origin is located at a moving boundary and that models the effluent pipe. A finite-difference scheme is proved to satisfy an invariant-region property (in particular, it is positivity preserving) if executed in a simple splitting way. Simulations are presented with a modified variant of the established activated sludge model no.~1 (ASM1).


Introduction
A sequencing batch reactor (SBR) is a tank (with possibly varying cross-sectional area) used for the purification of wastewater (e.g., [1,2,3]). It has a controlled outlet at the bottom and at the surface of the mixture, a floating device allows for controlled fill or extraction of mixture; see Figure 1. It is operated batch-wise in a sequence of cycles of fill, react, settle, draw and idle stages; see see Figure 2. The mixture consists of flocculated particles of biomass (activated sludge) consisting of several components that react with dissolved substrates (nutrients) in the liquid. (The removal of these substrates is the purpose of an SBR.) The full model of the process consists of a system of nonlinear partial differential equations (PDEs); see [4], to which we also refer for references to previous related works. To the authors' knowledge, models for SBRs with reactions during all stages (including the settle) have rarely been proposed in the literature. Models of reactive settling in continuously operated secondary settling tanks (SSTs) based on PDEs and numerical schemes are presented in [5,6,7]; see also references therein.
In addition to the mathematical and numerical difficulties of an SST model, the moving boundary in an SBR means a special challenge. Although in normal operation only liquid is extracted through the pipe in the draw phase, our model and numerical scheme are capable to handle any concentrations at the mixture surface.
The moving-boundary model has a connected half-axis with transport equations and nonlinear mass-preserving coupling conditions that do not define the coupling concentrations uniquely. Such a problem of nonuniqueness arises already for a scalar conservation law with discontinuous flux,  which has been investigated widely [8,9,10,11,12,13]; in particular, in the context of continuous sedimentation [14,15,16,17,18] where a monotone numerical scheme approximates the correct solutions [14]. It is the purpose of this contribution to present a numerical scheme that can handle the strong type degeneracy of the PDEs as well as the moving boundary where both a source is located and a half-line model attached. The scheme is monotone when the reaction terms are zero. We prove that if the scheme is used in the simple Lie-Trotter-Kato splitting way (e.g. [19]), namely, one explicit time step is taken without reactions and a another step with only reactions, then the numerical solutions have an invariant-region property under a convenient Courant-Friedrichs-Lewy (CFL) condition. In particular, concentrations are nonnegative.
The PDE model is presented in Section 2. Section 3 contains the numerical scheme based on the one proposed in [5] but modified to handle the moving boundary. The fully discrete and explicit scheme is presented in Section 3.4 and its splitting version with the CFL condition and invariantregion property is presented in Section 3.5. The numerical scheme during the full mixing stage can be found in Section 3.6. In Section 4, we show a numerical examples of SBR operation with a constant cross-sectional area (cylindrical vessel) and a commonly used activated sludge model for the biokinetic reactions; however, slightly adjusted to deliver non-negative concentrations only. Some conclusions can be found in Section 5.

The model
The vertical z-axis of the governing model and the moving coordinate system with the x-halfaxis are shown in Figure 1. The characteristic function γ equals one inside the mixture and zero otherwise, i.e., γ(z, t) = χ {z(t)<z<B} , where χ I is the indicator function which equals one if and only if I is true, and z =z(t) is the surface location. We let A = A(z) denote the cross-sectional area of the tank that may depend on depth z.
The solid phase consists of flocculated particles of k C types with concentrations C (1) , . . . , C (k C ) . The components of the liquid phase are water of concentration W and k S dissolved substrates of concentrations S (1) , . . . , S (k S ) . The total concentrations of solids X and liquid L are All these concentrations depend on z and t.
At the surface of the mixture, z =z(t), we model a floating device connected to a pipe through which one can feed the tank with a given volume rate Q f (t) and given feed concentrations C f (t) and S f (t); see Figure 1. Alternatively, this floating device allows to extract mixture at a given volume rate Q e (t) > 0 through the same pipe; hence, one cannot fill and extract simultaneously. If [0, T ] denotes the total time interval of modelling (and simulation in Section 3), we assume that When t ∈ T e , we model in [4] the extraction flow in the effluent pipe by a moving coordinate system, namely a half line x ≥ 0, where x = 0 is attached to z =z(t); see Figure 1. Since we assume that there are no reactions in the pipe and all components have the same velocity, the conservation law for the pipe is the linear advection equation A e ∂ tC + Q e (t)∂ xC = 0 where A e is the cross-sectional area of the pipe. However, we are only interested in the effluent concentrations in the pipe at x = 0 + , which we denote by C e (t) and S e (t). At the bottom, z = B, one can withdraw mixture at a given volume rate Q u (t) ≥ 0. The underflow region z > B is for simplicity modelled by setting A(z) := A(B), since we are only interested in the underflow concentration C u (t), which is an outcome of the model (analogously for S u (t)). Below the surface, the volume fractions of solids and liquid add to one, i.e., The same holds for the feed concentrations. For known C and S, (1) implies the water concentration This concentration is not part of any reaction and can be computed afterwards. For computational purposes, we define a maximal concentrationX of solids and assume that the density of all solids is the same, namely ρ X >X. Similarly, we assume that the liquid phase has the density ρ L < ρ X , typically the density of water. The velocity of a solid particle v X = q + v is the sum of the average volumetric velocity, or bulk velocity, of the mixture and the excess velocity v, which is given by the following commonly used expression [14,20]: where Here, ∆ρ := ρ X − ρ L , g is the acceleration of gravity, v hs = v hs (X) is the hindered-settling velocity, which is assumed to be decreasing and satisfy v hs (X) = 0, σ e = σ e (X) the effective solids stress, which satisfies σ e (X) = 0 for X ≤ X c and σ e (X) > 0 for X > X c , where X c is a critical concentration above which the particles touch each other and form a network that can bear a certain stress.
The reaction terms for all components are collected in the vectors which model the increase of solid and soluble components, respectively, where σ C and σ S are constant stoichiometric matrices and r(C, S) ≥ 0 is a vector of non-negative reaction rates, which are assumed to be bounded and Lipschitz continuous functions. We set (analogously forR S (C, S)).
In the present work the matrices σ C and σ S and the vector r are given by expressions that represent a slightly modified version of the activated sludge model no. 1, usually abbreviated 'ASM1' in wastewater engineering. The model was the first comprehensive activated sludge model developed by a task group of the International Water Association (IWA) (see [21]). The model describes the principal biochemical reactions that occur within the activated sludge process, which is the most widespread technology for the secondary treatment of municipal wastewater and as commented in [22], constitutes 'the heart' of many wastewater treatment plants. Later versions, known as ASM2, ASM2d, and ASM3, account for additional reactions such as fermentation and chemical or biological phosphorus removal (not considered in the present work). Commercial software packages that include these models are used commonly by wastewater process design engineers for the process design of various activated sludge system configurations [3]. We refer to Appendix A for the description of the modified ASM1 model used herein and to [1,3,22] for further information on activated sludge models.
In order to establish an invariant-region property for the numerical solution, we make some technical assumptions. To ensure that the numerical solution for the solids does not exceed the maximal concentrationX, we assume the following: v(X, ∂ z X, z, t) = 0.
These conditions mean that when the concentration is (near) the maximal one, biomass cannot grow any more and its relative velocity to the liquid phase is zero. To obtain positivity of component k of the concentration vector C, we let denote the sets of indices l that have negative and positive stoichiometric coefficients, respectively, and assume the following (analogously for S): if l ∈ I − C,k , then r (l) (C, S) =r (l) (C, S)C (k) withr (l) bounded.
Assumption (8) implies that . , k C (analogously for S), which means that the system of ODEs has a non-negative solution if the initial data are non-negative [23]. This positivity property is carried over to the numerical splitting scheme suggested here. The volume of the mixture is defined bȳ The function V is invertible since V (z) = −A(z) < 0; in particular, It turns out that the surface location is given by the following explicit expression of the given volumetric flows: Alternatively,z(t) can be obtained from To state the governing model, we define the velocities and then introduce the total mass fluxes as The complete model is the following (δ is the delta function): The water concentration W can always be calculated from (2). No initial data are needed for the outlet concentrations, but for C and S, namely During the react stage of an SBR (see Figure 2), full mixing occurs and the system of PDEs (17a) and (17b) reduces to the following system of ordinary differential equations (ODEs) for the homogeneous concentrations inz(t) < z < B: where all concentrations depend on time only since they are averages (below the surface). As before, W can be obtained afterwards from (2). In the region 0 < z <z(t) all concentrations are zero. Because of (15) and (17c) we have C u (t) = C(t) and C e (t) = C(t)χ {t∈Te} (analogously for S).

Spatial discretization and numerical fluxes
We divide the tank into N computational cells each having depth h = B/N . Assume that the midpoint of cell j has the coordinate z j , hence, the cell is the interval [z j−1/2 , z j+1/2 ]. The top cell 1 is thus [z 1/2 , z 3/2 ] = [0, h], and the bottom location is z = z N +1/2 = B. To obtain the underflow concentrations, we add one cell below z = B. To obtain the extraction concentrations, we add one cell [0, ∆x] of the x-coordinate system. To approximate the cell volumes, we define the average cross-sectional areas The unknowns are approximated by functions that are piecewise constant in each cell j, i.e.
, which are collected in the vector C j (t). We definej(t) := z(t)/h , which is the smallest integer larger than or equal toz(t)/h. Then the surface z =z(t) is located in the surface cellj(t).

Time discretization and surface fluxes
We let T denote the total simulation time, t n , n = 0, 1, . . . , N T , the discrete time points and τ := T /N T the time step that should satisfy a certain CFL condition; see below. The value of a variable at time t n is denoted by an upper index, e.g., C n j and it is thus assumed to be constant in time during t n ≤ t < t n+1 . The discrete surface index is defined byj n :=j(t n ) and we let z n :=z(t n ). For the volumetric flows, we define the averages and assume for simplicity that any of the volumetric flows changes sign at most at the discrete time points t n . This implies thatz(t) is monotone on every interval [t n , t n+1 ]. To ensure that the surface does not travel more than one cell width h during τ , the CFL condition has to imply (cf. To show that the cell concentrations C n j n do not exceed the maximal oneX, we also introduce the concentrationC n j n obtained when all the mass in the surface cellj n is located below the surface within the cell; cf. Figure 3 (a). The mass in the cell is C n j n Ajn α n h = C n j n Ajn h, where α n h := zjn +1/2 −z n . We setX n j :=C If the surface stays within one cell between t n and t n+1 ; then (22) is equivalent to (cf. Figure 3 (b) and (e)). During extraction, the surface cannot rise and is thus located somewhere in cellsj n andj n + 1 (cf. Figure 3 (d) and (e)). In light of (15) and (16), we approximate the fluxes, which have to be non-positive, just below the surface in the following way: Φ S,n e,j n +1/2 := −

Derivation of update formulas
We here derive the update formulas for C n j . Analogous formulas hold for S n j when replacing C by S; however, with different definitions of velocities and fluxes. First come cells that lie below the surface z =z(t) at t n and t n+1 . Special treatment is needed for the cells near the surface. All cells strictly above the surface have zero concentrations. Let κ := τ /h.

Cells away from the surface
Using the integrated form of the balance law on a rectangle [z j−1/2 , z j+1/2 ] × [t n , t n+1 ] strictly below the surface; see the gray rectangle in Figure 3 (a), we get the update formula (mass per h) and the analogous one for S n j . For cell N + 1, we get in the similar way the update formula for the underflow concentration; see Section 3.4.

Cells near the surface during fill (t ∈ T f )
To obtain a monotone scheme with an invariant-region property, we determine how the mass in the surface cell and the one below evolves; see the trapezoids in Figure 3. The mass per h at t n is (we use (21)) m C,n j n = Ajn α nC n j n + Ajn +1 C n j n +1 = Ajn C n j n + Ajn +1 C n j n +1 .
During τ , the feed source along the moving surface is Q n f C n f , whereas the outflux is Φ C,n j n +3/2 . Thus, by the balance law on any trapezoid the mass (per h) at t n+1 is where the in-and outflux and source terms are (the concentration in cellj n below the surface isC n j n ). Case (a): Fill casej n =j n+1 − 1, Figure 3 (a): The surface moves downwards and crosses a cell boundary. All the mass ends up in one cell: Figure 3 (b): When the surface does not cross any cell boundary during τ , the mass (28) is distributed among the two cells with respect to their volumes (below the surface): Case (c): Fill casej n =j n+1 + 1, Figure 3 (c): After the balance law is used on the purple trapezoid, the final mass (per h) is distributed among three cells: Cells near the surface during extraction (t ∈ T e ) During extraction, the surface necessarily moves downwards. The initial mass is (27) and the balance law on a red trapezoid (Figures 3 (d) and (e)) gives m C,n+1 j n = m C,n j n + κΨ C,n e,j n , where Ψ C,n e,j n := Φ C,n e,j n +1/2 − Φ C,n j n +3/2 + h α n Ajn R C C n j n α n , S n j n α n + Ajn +1 R n C,j n +1 .

The cell in the extraction pipe
The conservation law for the cell on the x-axis gives the mass equality (Figure 3 (d)) where the cross-sectional area A e of the effluent pipe is of less importance, since we are only interested in C n e and therefore may choose any ∆x; we set A e ∆x := A 1 h.

Explicit fully discrete scheme
Given data at t n and the valuesz n+1 andj n+1 , the update formulas for the particulate concentrations are given here and we distinguish between fill and extraction. We define λ j := κ/A j = τ /(A j h) and The update formulas for the top cells below the surface are different depending on whether fill or extraction occurs.
The effluent concentration is given by

Other concentrations
For the cells j =j n + 2, . . . , N , the update formula is (analogously for S), at every time point t n , where the numerical flux is computed by (19). For the cells above the surface, we have C n+1 j = 0 for j <j n − 1. Finally, one computes ,n j + · · · + S (k S ),n j .

A splitting scheme and invariant-region property
It is desirable that the solution vectors U := (C, S) and U e := (C e , S e ) of the model (17) stay in the set Ω := U ∈ R k C +k S : C ≥ 0, S ≥ 0, C (1) + · · · + C (k C ) ≤X .
To ensure that Ω is an invariant set for the numerical solutions, we split the scheme in Section 3.4 by taking one time step without the reaction terms and then one time step with only the reaction terms. For a cell strictly below the surface and the solid concentrations, the update formula (30) can be written C n+1 j = C n j − λ j [∆Φ C ] n j + τ R C (C n j , S n j ). and similarly for S n j . The splitting principle is to compute the expression in the curled brackets first and then use the result for the second step as follows: For the cells involving the surface, the first-step update formulas, j ∈ {j n − 1,j n ,j n + 1} arě and the analogous formulas with C replaced by S. The second step consists in the formulas and the analogous formulas with C replaced by S. The time step τ has to be bounded by the CFL condition where β 1 and β 2 depend on h, h 2 , the volumetric flows, and the constitutive functions by and where the constants are given by (here, ξ represents v hs , v hs or d) (cf. (8)) Theorem 1. Consider the numerical splitting method in Section 3.5. If U n j := (C n j , S n j ) ∈ Ω for all j =j n ,Ū n j n := (C n j n ,S n j n ) ∈ Ω, U n e := (C n e , S n e ) ∈ Ω and (CFL) holds, then Theorem 1 is proved by the following lemmas and by the fact that (CFL) implies (20), which we have used in the derivation of the scheme. For the proofs write the update formulas as whereΦ C andR n C,j denote the sum of all components, respectively. Analogous considerations lead to formulas forX n j and X n+1 e . Lemma 1. Let R C ≡ 0. If U n j := (C n j , S n j ) ∈ Ω for all j =j n ,Ū Proof. We write the general update formula (38) as X n+1 j = H n X,j (X n j−1 , . . . , X n j+3 ) and let this include the surface concentrationX n j n . For cells away from the moving surface, we refer to [5, Theorem 3.1] from which we also collect v X,n,+ j+1/2 = q n j+1/2 + γ j+1/2 v hs (X n j+1 ) − J C,n j+1/2 To show the monotonicity for the cells near the surface, we set χ + := χ {v X,n j n +1/2 ≥0} , χ − := χ {v X,n j n +1/2 ≤0} , χ − e := χ {v X,n e,j n +1/2 ≤0} , and calculate for j =j n ∂v X,n,± j+1/2 ≥ 0.
For ease of notation, we introduce and let as usual tilde denote the sum of all components of a vector. In the case t n ∈ T num f , all the coefficients ω f,n j forΥ C,n f,j in (34) are non-negative, so it suffices to show that the derivatives of this function are non-negative under (CFL): In the case t n ∈ T num e , we first estimate This estimation shall be added to the derivatives ofΨ C,n f,j to obtain those forΨ C,n e,j . Introducing we see that the only derivative that differs from those above is We write C Ce (X n e , X n j n +1 ) for any fixed k ∈ {1, . . . , k C }. This formula is trivial for t n ∈ T num f , and for t n ∈ T num e , we get To prove the boundedness, the monotonicity in each variable of H n X,j and the assumption (7) are used to obtain, for t n ∈ T num f and j =j n =j n+1 , For t n ∈ T num f and j =j n − 1 =j n+1 , we first see that (22) implies; cf. Figure 3(c), which we use at the end of the following estimate: = H n X,j n −1 0, 0, X n j n , X n j n +1 , X n j n +2 ≤ H n X,j n −1 0, 0, α nX ,X,X In the case j =j n =j n+1 + 1 (cf. Figure 3(c)) we can still use (41) to obtain 0 ≤ θ n+1 λjn Q f X n f = H n X,j n (0, . . . , 0) ≤ X n+1 j n = H n X,j n 0, 0, X n j n , X n j n +1 , X n j n +2 ≤ H n X,j n 0, 0, α nX ,X,X A similar estimation can be made for the case j =j n +1 =j n+1 +2. For the case j =j n +1 =j n+1 ; see Figure 3(a), we note that (22) implies which we use to estimate = H n X,j n +1 0, X n j n , . . . , X n j n +3 ≤ H n X,j n +1 0, α nX ,X,X,X The remaining fill cases are similar; we omit details. For t n ∈ T num e , similar estimations apply; the only difference is that Q f X n f is replaced byΦ C,n e,j+1/2 , which equals zero when the concentrations are zero, to prove the lower bound. For the upper bound, one uses (22) withQ n = −Q n e instead of Q n f . For the effluent, we get 0 = H n Ce (0, 0) ≤ X n+1 e = H n Ce X n e , X n j n +1 ≤ H n Ce X ,X =X − λ 1 (Q n e + Q n u )X ≤X.
Lemma 2. Let R C ≡ 0. If U n j := (C n j , S n j ) ∈ Ω for all j =j n ,Ū n j n := (C n j n ,S n j n ) ∈ Ω, U n e := (C n e , S n e ) ∈ Ω and (CFL) holds, then Proof. This follows directly from Lemma 1, since each component of the update formula for C n j is equal to that of X n j (there is no reaction term). Proof. We start as in the proof of Lemma 1, use the notation and estimations from there, and prove monotonicity of each component of the right-hand side H (k),n S,j n of the update formula for component S (k),n j , which we write as S n j . We also skip the superscript (k) for components of other vectors. We letρ := 1/(ρ X −X) and The numerical fluxes are different and we get Because of the similarities between Ψ C,n f,j and Ψ S,n f,j , we only write the difference here: To prove the monotonicity in the case t n ∈ T num f , we conclude that the estimations are in fact similar to those in the proof of Lemma 2 with the following difference (where Υ S,n f,j is defined analogously to (39) and Υ S,n f,j denotes an arbitrary component of that vector): To prove the monotonicity in the case t n ∈ T num e , we first estimate Following the same procedure as in the proof of Lemma 2 (with Υ S,n e,j is defined analogously to (40)), we now get For the update of the effluent concentrations, we get the same result for S n e as for C n e . The proof of positivity can be done as in the proof of Lemma 2. Proof. Component k of (33) can be written and estimated by means of (5) and (8) as where the latter inequality follows from (CFL). Consider now the formula (36) for the fill case (a) with the non-zero coefficient ω f,n j n +1 := Ajn/Ajn +1 (the others are trivial). Then component k of (36) is estimated by using the rewriting above applied to both reaction terms: where the latter inequality follows from (CFL). For fill case (b), we get in the case ω f,n j n +1 := η n+1 The other cases are similar, for example, in case (c) with ω f,n j n +1 := θ n+1 one gets The cases of extraction are the same and so are the analogous ones for S n j+1 . Lemma 5. Consider the second step of the splitting scheme with formula (33) and the formula obtained by summing all vector components. IfX n+1 j :=Č (1),n+1 + · · · +Č (k C ),n+1 ≤X and (CFL) holds, then the second step of the splitting scheme corresponding to (33) satisfies The analogous statement is true for the other update formulas near the surface.
Proof. IfX n+1 j ≥X − ε, thenR C (Č n+1 j ,Š n+1 j ) = 0 by assumption (6). OtherwiseX n+1 j <X − ε and (CFL) gives that X n+1 j is estimated by The rest of the update formulas are treated in the same way and the resulting formula contains an additional factor ω f,n j ∈ [0, 1] multiplying τ , so the result also holds for that formula.

Numerics for full mixing
During the react stage, mixing occurs due to aeration or the movement of an impeller; see Figure 2. Then there is no relative velocity between the solids and the liquid, but reactions take place, and the time-dependent concentrations for the mixture below the surface are given by the ODEs (18).
Suppose the (PDE or numerical) solution C(z, T 0 ) is known at t = T 0 = t n0 when a period of complete mixing starts. The initial concentrations for the ODEs (18) are defined as the averages (for k = 1, . . . , k C ; analogously for S) aver .
The ODE system (18) can then be time integrated. If an ODE mixing period ends at t = tñ and the PDE model is to be simulated thereafter, then the total mass below the surface is distributed among the cells by

Material
Notation Unit Particulate inert organic matter X I (g COD) m −3 Slowly biodegradable substrate X S (g COD) m −3 Active heterotrophic biomass X B,H (g COD) m −3 Active autotrophic biomass X B,A (g COD) m −3 Particulate products arising from biomass decay X P (g COD) m −3 Particulate biodegradable organic nitrogen Recall that the ODEs (18) are averages of the PDE model (17a) and (17b) when the convective and diffusive terms are zero [4]. Therefore, if time integration is made with the explicit Euler method and (CFL) holds, then also the numerical approximations of the ODE solutions to (18) belong to Ω.

Numerical simulations
The numerical method in Section 3.4 with condition (CFL) is first used for the simulation of an SBR process. In [4], we demonstrated the process by letting the reaction terms model a denitrification process, which occurs when there is no oxygen present. Here, the reactions are a modified ASM1 model without alkalinity; see Table 1 for the six particulate and six soluble state variables, and Appendix A for the reaction terms. Furthermore, we assume that the mixing during the react stage is achieved by aeration, which means that the liquid is saturated with dissolved oxygen; a concentration about S O = 10 g/m 3 . The constitutive functions used for sedimentation and compression are We simulate one sequence of an SBR with the stages specified in Table 2. (This is the same scenario as in the first example of [4] with a denitrification reaction model and without oxygen supply.) Condition (CFL) implies the step length τ = 4.3716 · 10 −5 s. The initial concentrations have been chosen as  where the total solids feed concentration X f (t) varies with time according to Table 2. Figures 4 and 5 show the simulation results for the concentrations within the vessel of the particulate and soluble components, respectively. The numerical scheme resolves all the discontinuities accurately. At time t = 1 h, the tank has been filled and the react stage starts and there is full mixing by aeration. Consequently, the available components in the tank are distributed homogeneously. During the react stage, the growth of biomass is slow but the fast consumption of S S and increase of S NH are visible in Figures 5(b) and (e), respectively. The nitrification process uses oxygen to produce nitrate and nitrite, which can be seen in plot (d). After the react stage, t ≥ 3 h, the oxygen is quickly consumed (plot (c)), but only where there is biomass and S S and S NH are positive.
The convergence of the numerical scheme is demonstrated in Figure 6, where some of the effluent concentrations are plotted. As was illustrated with the example of denitrification in [4] and the present one for the modified ASM1 model, simulations seem to satisfy the invariant-region property although we only have a proof of this for the splitting scheme of Section 3.5.
For the ASM1 example here and a given number of cells N , we calculate the L 1 relative difference between the simulation result of the two variants of the numerical scheme according to the following formula, where C N is the result without splitting and C split N with splitting: The result at time t = T = 6 h is shown in Table 3. The relative difference is approximately halved as N is doubled.

Conclusions
The general model of multi-component reactive settling of flocculated particles given by a quasione-dimensional PDE system with moving boundary in Section 2 was derived in [4]. The unknowns are concentrations of biomass particles and soluble substrates, and the reaction terms of the PDE model can be given by any model for biochemical reactions in wastewater treatment.
The numerical scheme in Section 3.4 is designed to ensure conservation of mass across the surface during fill and draw. Away from the moving boundary, the scheme is the same as in [24], where it is demonstrated that its order of convergence is not more than one. The extra treatment near surface will of course not improve that. An indication of the convergence of the numerical scheme as the mesh size is reduced is demonstrated in Figure 6. The main result of this work is an invariant-region property (Theorem 1) for the numerical solution if the scheme is computed in a Lie-Trotter-Kato splitting way where the first time step is taken without any reactions and then a step with only the reactions. Then all the concentrations are nonnegative and the solids concentrations never exceed the maximal packing one. In particular, the scheme is monotone when the reaction terms are zero. Simulations with or without splitting have shown to produce very similar outputs and this is demonstrated by the diminishing relative error between such simulations in Table 3.
A proof of convergence of the method (as h → 0) to a suitably defined weak or entropy weak solution, as well as the corresponding well-posedness (existence and uniqueness) analysis, are still pending. That said, we point out that available convergence analyses for related strongly degenerate, scalar PDEs with discontinuous flux (cf., e.g., [12,13,14]) rely on the monotonicity of the underlying scheme as well as a uniform bound on the numerical solution, among other properties. Theorem 1 and its proof may be therefore viewed as a partial result to prove convergence of the numerical scheme presented herein.
Given a moving boundary and a fixed spatial discretization for the numerical scheme, local mass balances have been used to obtain correct update formulas for numerical cells near the surface. This results in a scheme with several cases depending on the surface movement. A certain limitation of the explicit numerical scheme, used without or with splitting, is the restrictive CFL condition (where the time step is esentially proportional to the square of the cell size), implying that very small time steps are needed if accurate approximations on a fine spatial mesh are sought. An alternative approach would be to transform the PDE system and have a fixed number of cells below the moving surface. Such a scheme could possibly also be easier to generalize to a high-order scheme or a more efficient one with semi-implicit time discretization. The advantage of the present fixed-cell-size numerical scheme is, however, that the model can more easily be generalized to include further sources or sinks at fixed locations, a desirable feature in applications to wastewater treatment.