The main purpose of the recently introduced Bürger–Diehl simulation model for secondary settling tanks was to resolve spatial discretization problems when both hindered settling and the phenomena of compression and dispersion are included. Straightforward time integration unfortunately means long computational times. The next step in the development is to introduce and investigate time-integration methods for more efficient simulations, but where other aspects such as implementation complexity and robustness are equally considered. This is done for batch settling simulations. The key findings are partly a new time-discretization method and partly its comparison with other specially tailored and standard methods. Several advantages and disadvantages for each method are given. One conclusion is that the new linearly implicit method is easier to implement than another one (semi-implicit method), but less efficient based on two types of batch sedimentation tests.

## INTRODUCTION

Benchmark simulations of entire wastewater treatment plants (WWTPs) are today performed with one-dimensional simulation models of the secondary settling tank (Gernaey *et al.* 2014; Li & Stenstrom 2014). In the model by Bürger *et al.* (2011, 2013), sometimes referred to as ‘Bürger–Diehl model’, the physical phenomena of hindered settling, volumetric bulk flows, compression of the sludge at high concentrations and dispersion of the suspension near the feed inlet can be included in a flexible way. Each phenomenon is associated with a separate constitutive function with its model parameters, and can be activated or deactivated at the user's discretion. The possibility to include sludge compression is particularly important, since this improves the predictive power considerably (De Clercq *et al.* 2008; Ramin *et al.* 2014; Guyonvarch *et al.* 2015; Torfs *et al.* 2015; Van Loosdrecht *et al.* 2015). However, the inclusion of physical phenomena that result in second-order derivative terms in the model partial differential equation (PDE), e.g., compression, means that straightforward time discretization by easy-to-implement, explicit methods such as the explicit Euler method leads to long computational times.

Bürger *et al.* (2013) presented implementation details of a numerical algorithm, which gives reliable simulations with respect to the underlying physical principles and can be obtained with a user-defined accuracy. For long-time simulations of entire WWTPs, it is important to keep discretization errors small. This calls for a fine resolution in both space (many layers in the settler) and time (short time steps), which implies long computational times. Conversely, fast computations obtained with a low resolution in space and time come at the cost of poor accuracy.

In scientific computing, the numerical error is a measure of how close a numerical solution is to the exact solution or reference solution (obtained by a very fine discretization) of the model, i.e., the governing differential equation. The efficiency of a numerical method is assessed by relating the numerical errors to the computational (central processing unit; CPU) times necessary to obtain the numerical solutions for different discretizations. It is then said that one numerical method is more efficient than another if it allows a numerical solution with a determined numerical error to be obtained in less CPU time, or equivalently, a given budget of CPU time allows one to obtain a more accurate numerical solution by the first method than by the second.

The simulation model by Bürger *et al.* (2013) is based on a method-of-lines formulation of the underlying nonlinear PDE. This means a system of time-dependent ordinary differential equations (ODEs), one for each layer of the settler. Simulations of PDEs are stable and reliable if a so-called CFL (Courant–Friedrichs–Lewy) condition is satisfied. This gives a maximal time step for each given layer thickness . If only hindered settling and bulk flows are included, then the CFL condition means that can be chosen proportional to , i.e., , which results in fast simulations. When compression or dispersion is included and standard ODE solvers used, the CFL condition states that , which means very small when the error should be reduced (small is chosen).

*z*is measured from the suspension surface downwards to the bottom at . The model PDE iswhere the flux function contains the hindered settling velocity function and the compression function is , where

*K*is a constant containing the solid and fluid mass densities and the acceleration of gravity, and is the effective solid stress function, which is zero below a critical concentration and increasing above.

The investigated methods are used for the simulation of two different batch settling tests in a vessel with m. While the present study is limited to batch settling in a column, the method proposed herein should be appropriate to handle secondary settling tanks as well (under suitable modifications that would mainly affect the convective flux, but not the compression term for whose discretization the present method has been tailored). It is therefore essential to demonstrate that the method can robustly handle several scenarios, in particular a scenario that exhibits some spatial variation of the solution. For this reason we consider a conventional Kynch test (KT; Kynch 1952), that is, the settling of an initially homogeneous suspension, and a Diehl test (DT; Diehl 2007), where a body of concentrated suspension is initially located above clear liquid, separated by, for example, a membrane. The DT shows a stronger variation of the solution than the KT. Both tests provide complementary information that can be used to identify large portions of the flux function *f* (Betancourt *et al.* 2014).

## METHODS

^{3}/kg, and the effective solids stress chosen iswhere m

^{2}/s

^{2}and kg/m

^{3}. To obtain a working numerical method an important preparation is to compute the primitiveof the compression function . In the present case, with the functions and chosen herein, we can obtain in closed algebraic form; if this is not possible (for other choices of these functions), the function can be obtained numerically (Bürger

*et al.*2013). Now (1) can be written asFor batch sedimentation, this PDE should be complemented with initial data and zero-flux boundary conditionsSuitable numerical schemes for the approximate solution of (1) are based on subdividing the depth interval into a number

*N*of layers of equal thickness . A discretization in space leads to the following method-of-lines formulation (Bürger

*et al.*2013), which is a system of

*N*ODEs:where is the concentration of layer

*j*, denotes the numerical convective flux due to hindered settling between layer

*j*and layer , which can be chosen, for instance, as the Godunov flux, and is the compressive numerical flux chosen asNote that (2) implies that the boundary fluxes . Implementation details on how to compute these fluxes are provided by Bürger

*et al.*(2013). The last term in (4) now becomes the usual second-order difference approximation of , i.e.,The method-of-lines formulation (3)–(5) is converted into a fully discrete scheme by time discretization. Although a variety of methods could be used for (3)–(5), there is nothing to gain by applying a standard ODE solver of accuracy higher than first order, e.g., a Runge–Kutta method. On the contrary, this will only cause longer computational times producing the same error as explicit Euler time stepping (Diehl

*et al.*2015).

For the discussion of the various time-stepping methods, which will take the numerical solution from to , it will be important to carefully distinguish between numerical fluxes and concentrations that are evaluated at the old time step and those evaluated at , which we mark by the respective upper index *n* and . The time-integration methods compared herein are the following, where are constants that depend on the choice of the functions *f* and *d*.

### The explicit Euler method

*et al.*(2011, 2013). The CFL condition of the fully discrete scheme can be captured by when is small. The method is easy to implement since all terms in the right-hand side of (3)–(5) are evaluated at , and therefore each concentration is an explicit function of the known ones at the previous time . Thus, the fully discrete version of (4) iswith analogous formulas replacing the boundary updates (3) and (5).

### The semi-implicit method

*et al.*2006) is described in detail by Diehl

*et al.*(2015). The CFL condition is . To advance the solution from to we must solve the following nonlinear system of algebraic equations, where , , are the unknowns:supplemented with analogous formulas replacing the boundary updates (3) and (5). These equations are solved iteratively, for example, by the Newton–Raphson method.

### The linearly implicit method

*et al.*(1979) and is based on first considering the contribution from the compressive flux term (6) from time to . The purpose of the LI method is to avoid the numerical solution of the nonlinear system of Equation (7). The nonlinearities in (7) are found in the evaluations of the function

*D*. The idea is to replace in (7) by (where is a parameter connected to the convergence of the method) so thatand to find a simple update formula for , , to be executed first in each time step. A stable implicit Euler time step implies the formulaThis means that a linear system of equations should be solved for at the next time step. Then (8) can be updated explicitly to obtain .

To advance the solution from to , one proceeds as follows:

For , set .

*et al.*2015). A preliminary analysis on the stability of the scheme (not presented here) implies that the CFL condition of the LI scheme is

### Other time-stepping methods

Any ODE solver could be used for the method-of-lines system (3)–(5). Possibly competitive methods are adaptive step-size methods and we choose here the ODE solver ode15s in Matlab (2014), which is an implicit multi-step method of variable order with step-size control. Of several standard ODE solvers investigated by Diehl *et al.* (2015), this was the second most efficient one (after the SI method) in the investigations with both stand-alone settler simulations and benchmark simulations for the entire activated sludge process.

## RESULTS AND DISCUSSION

*N*layers up to a time

*T*is calculated by comparing with a reference solution obtained by Euler's method and layers . The relative error is calculated asFor the KT we choose an initial concentration of and , which means that almost steady state has been reached. This test is also used to study the dependence of the numerical solutions produced by the LI method on the choice of the parameter . For each run, the time step is chosen according to the CFL condition (10) (with equality). Note that the right-hand side of (10) is a decreasing function of and it tends to infinity as approaches 1 from above. Figure 1 shows numerical solutions obtained for three different values of . This figure and Figure 2 indicate that for a given number of layers , the values of closest to 1 produce the solutions with the smallest errors, but the highest CPU times. We note that, for instance, for a given value of

*N*the time step for can be chosen twice larger than that for , so the total CPU time for should be about half of that for . To assess which value of is optimal, an efficiency plot (of relative error versus CPU time) is helpful, see Figure 2. The curves are roughly L-shaped, indicating that should be chosen close to the point corresponding to the point of ‘bend’, since smaller or larger values would lead to larger numerical errors (at only slightly smaller CPU times) or to smaller numerical errors at much increased CPU times. Based on these considerations, the choice seems a good compromise for further comparison with other methods.

*N*, the LI method is the fastest, although not necessarily the one with the smallest error. Note that the efficiency curve for the SI method for the DT (the right plot of Figure 3) is composed of four symbols only. In fact, no information is available for the run with since the Newton–Raphson iterations did not converge for that case. While this situation could be easily overcome by the ad–hoc remedy of further reducing (say, to 80% of its maximal value determined by the CFL condition), it alerts us to a more fundamental problem observed with the SI method, namely that the convergence of (iterative) solvers for the nonlinear equations is not ensured a priori.

For accurate simulation results, should be chosen small, whereby the CFL condition implies that the explicit Euler method requires much smaller time steps than the other methods. On the other hand, the SI method needs more computations at every time step and requires the evaluation of the Jacobian matrix of the nonlinear algebraic system of equations. However, even when this system is solved by the Newton–Raphson method, the SI method has turned out to be far more efficient than explicit Euler (Diehl *et al.* 2015). To reduce the implementation complexity and to remove the necessity to solve numerically a system of nonlinear algebraic equations at every iteration in the SI method, whose convergence depends on the Newton–Raphson iterations, the LI method can be used instead. The advantage of solving a linear system only in each iteration, which also means an easier implementation, is paid fort by the price of a larger numerical error.

## CONCLUSIONS

With the aim of finding an efficient time-integration method for the simulation of sedimentation with compression, we have in this work confined ourselves to two batch sedimentation tests for the comparison of some methods. The ideas of the new LI method and its implementation have been provided. When evaluated on the batch settling tests, the following advantages and disadvantages can be concluded.

Explicit Euler method:

+ Implementation easy.

+ Convergence proof exists also for continuous sedimentation (numerical approximate solutions converge to the PDE solution as ) by Bürger

*et al.*(2005). Robust method.− Not efficient.

SI method:

+ Most efficient of the investigated methods. Under the assumption that the Newton–Raphson iterations find a solution each time step, convergence of the numerical method was proved for batch sedimentation by Bürger

*et al.*(2006).− Implementation more complex than Euler and LI. At each time step, a nonlinear system of algebraic equations is solved, e.g., by Newton–Raphson iterations, which require tolerance parameter to be set. There is no guarantee that these iterations converge.

LI method:

+ Implementation easy. The ingredient in addition to the Euler method is basically that a linear system of equations is solved at each time step. Robust method.

± Second most efficient for for batch sedimentation. The efficiency can be adjusted to some extent by a parameter. Fastest method for a given , but least accurate. Convergence proof in preparation.

Matlab's ode15s:

+ Ready-to-use standard time solver.

− Implementation most complex of the investigated methods. Well-established robust ODE solver for stiff problems, but not developed for solutions containing discontinuities.

− Least efficient of the investigated methods.

## ACKNOWLEDGEMENTS

RB is supported by Fondecyt project 1130154; Conicyt project Anillo ACT1118 (ANANUM); Red Doctoral REDOC.CTA, MINEDUC project UCO1202; BASAL project CMM, Universidad de Chile and CI^{2}MA, Universidad de Concepción; and Fondef ID15I10291. RB and CM are also supported by CRHIAM, Proyecto Conicyt Fondap 15130015. SD acknowledges support from The ÅForsk Foundation, Sweden.