Accurate voidage prediction in fluidisation systems for full-scale drinking water pellet softening reactors using data driven models

In full-scale drinking water production plants in the Netherlands, central softening is widely used for reasons related to public health, client comfort, and economic and environmental benefits. Almost 500 million cubic meters of water is softened annually through seeded crystallisation in fluidised bed reactors. The societal call for a circular economy has put pressure on this treatment process to become more sustainable. By optimising relevant process conditions, the consumption of chemicals can be reduced, and raw materials reused. Optimal process conditions are feasible if the specific crystallisation surface area in the fluidised bed is large enough to support the performance of the seeded crystallisation process. To determine the specific surface area, crucial variables including voidage and particle size must be known. Numerous models can be found in the literature to estimate the voidage in liquid-solid fluidisation processes. Many of these models are based on semi-empirical porous-media-based drag relations like Ergun or semi-empirical terminal-settling based models such as Richardson-Zaki and fitted for monodisperse, almost perfectly round particles. In this study, we present new voidage prediction models based on accurate data obtained from elaborate pilot plant experiments and nonlinear symbolic regression methods. The models were compared with the most popular voidage prediction models using different statistical methods. An explicit model for voidage estimation based on the dimensionless Reynolds and Froude numbers is presented here that can be used for a wide range of particle sizes, fluid velocities and temperatures and that can therefore be directly used in water treatment processes such as drinking water pellet softening. The advantage of this model is that there is no need for applying numerical solutions; therefore, it can be explicitly implemented. The prediction errors for classical models from the literature lie between 2.7 % and 11.4 %. With our new model, the voidage prediction error is reduced to 1.9 %.


Drinking water softening
Water softening involves the removal of calcium, magnesium, and other metal cations from water [1]. Central softening of drinking water is currently frequently applied in several countries (e.g. the Netherlands, Belgium, Germany, France, and the USA) while domestic softening is the most frequently applied way of softening in other countries [2][3][4]. In full-scale drinking water production plants in the Netherlands, central softening is widely used for reasons related to public health, client comfort, and economic and environmental benefits [5][6][7]. Specifically in areas with high water hardness, centralised drinking water softening can reduce the consumption of soap, detergents, and other household chemicals and increase the service life and energy efficiency of household appliances such as coffee machines due to a reduction in calcium carbonate scaling [8,9].
In the Netherlands, in 2020 almost 500 million cubic meters of https://doi.org/10.1016/j.jwpe.2020.101481 Received 6 May 2020; Received in revised form 17 June 2020; Accepted 22 June 2020 water is being softened annually through seeded crystallisation in fluidised bed reactors, a process known as pellet softening [10], primarily removing calcium cations from the water. Pellet softening was developed and introduced in the Netherlands in the 1970s [11]. In other fields, crystallisation in fluidised bed reactors is also becoming increasingly popular. Examples include the removal of hardness from natural hard ground waters [12], the recovery of CO 2 in direct air capture processes [13], the improvement of pellet characteristics for reuse potentials [14], reduced sludge production [15], groundwater softening in circulating pellet fluidised bed reactors usage in thermal power plants [16], and organic micropollutant removal from groundwater [17]. There is also a growing interest in fluidisation of biomass particles [18], fluidised bed reactors used in wastewater treatment [19], and other liquid-solid fluidisation techniques with many applications in engineering [1,20,21]. In pellet softening [22], drinking water is treated in an up-flow fluidised cylindrical bed reactor (with flow velocity of 60−120 m/h). The Amsterdam reactor [23] is currently a widely applied fluidised bed reactor (60−90 m/h). A strong chemical base, usually caustic soda or calcium hydroxide combined with a seeding material, is dosed to obtain supersaturated conditions in the water phase and is subsequently well mixed in the reactor. High pH leads to an alteration of the calcium carbonate equilibrium in which the solubility product is exceeded and CaCO 3 precipitates on the surface of the seeding material, thereby forming pellets. Seeds, typically quartz sand [24], garnet sand [25] or limestone calcite, are introduced in the reactor and fluidised under upward flow. The calcium carbonate crystallisation on the seeds mainly occurs in the lower region of the reactor [11,23]. Due to the crystallisation, the calcite pellets grow and migrate to the bottom of the reactor. When a defined grain size threshold is exceeded, a small part of the pellets is extracted from the bottom of the reactor [26].
The calcite pellets that are produced as a by-product in these plants can be reused as seeding material, after grinding and sieving of the produced pellets [27]. Pellets are the main by-product from the softening process and by identifying reuse applications, they potentially represent a resource rather than a waste product, thus promoting the establishment of a circular economy. The main advantages of this reuse are an expected significant decrease of the ecological footprint of both calcite as raw material and the drinking water treatment plants using pellet softening and an increased valorisation of the pellets [28]. Process optimisation [26] and control [29] has been focussed primarily on pellets grown on garnet sand as seeding material. The transition from garnet sand as a seeding material to reused crushed calcite took place at Waternet from 2014 to 2016.

Sustainability goals
There is a wide range of definitions of the concept of sustainability. The World Commission on Environment and Development [30] defines sustainability as the 'development that meets the needs of current generations without compromising the ability of future generations to meet their needs and aspirations'. The importance of continuing the development of a worldwide clean and sustainable water supply is increasing [31]. In general, drinking water suppliers have become more and more focussed on appropriate sustainable treatment technologies [32]. Optimal process conditions of water treatment processes contribute to the sustainability goals of public water organisations since fewer chemicals, less energy, and fewer raw materials are needed [33,34].
The city of Amsterdam, the Netherlands, has the ambition to develop itself as a competitive and sustainable European metropolis [35]. The flows of energy, water, and resources within the urban environment offer a large potential to contribute to this ambition through a transition from the linear usage of resources and waste production towards the sustainable management of urban resources with circular flows of resources. Currently, the detrimental contributions of pellet softening are in particular related to the consumption of chemicals like caustic soda and raw materials such as garnet sand, mined in Australia and shipped to the Netherlands, both of which contribute to the carbon footprint and the environmental burden presented by the Dutch drinking water companies [8]. A life cycle assessment study has shown that pellet softening can be improved in terms of eco-efficiency [36]. The role of more sustainable pellet softening is twofold: it concerns the transition from garnet sand to reused calcite pellets [28] as well as the reduced consumption of caustic soda [8] in an optimal fluidised bed reactor. Both aspects imply the need for more accurate information of the fluidised bed conditions. To be able to acquire the process state of the fluidised bed, the most important process variable, i.e. the effective specific surface area (SSA) for the purpose of crystallisation, must be known. Therefore, the effective voidage must be determined. SSA is either based on the total reactor volume or on the water phase in the reactor, to be discussed in Section 3.2. An optimal operational configuration will lead to a sustainable operational approach for pellet softening by using as little chemicals and raw materials as possible and ultimately leading to a more environmentally sustainable drinking water supply [14].
Waternet, the public water utility of Amsterdam and surroundings, is seeking a sustainable scenario for producing drinking water and offering services that fulfil the requirements of clients and regulations, and, at the same time, maintaining a sound environmental performance while keeping costs as low as possible [37].
To meet sustainability goals and to promote the development of a circular economy, Waternet has modified its pellet-softening processes [28], in which garnet grains have been replaced by calcite seeding particles that are based on crushed, dried, sieved, and reused calcium carbonate pellets. The garnet core inside the pellets hinders their potential application in market segments such as the glass, paper, food, cosmetics, and feed industries. The pellet market value and the sustainability of the softening process can be increased through the substitution of the sand grain by a calcite grain of 0.5 mm (100 % calcium carbonate). If the calcite pellets are crushed, dried, and sieved, they can be reused as a seeding material [27]. To reduce the required amount of chemicals bases like caustic soda, a large crystallisation SSA for optimal crystallisation purposes is an essential condition. Smaller grains imply a larger SSA, but grains that are too small have a potential risk to be flushed out of the fluidised bed reactor.

Reactor technology conditions and research objective
The efficiency of the entire pellet softening process can be expressed in terms of the amount of supersaturated calcium in the effluent, which can be determined by using the calcium carbonate crystallisation potential (CCCP). The CCCP determines the amount of chemicals used and the corresponding costs and CO 2 -equivalent footprint of the pellet softening process.
Chemical yield is defined as the amount of desired product produced relative to the amount that would have been formed if there were no by-products and the main reaction went to completion. For pellet softening, a high yield implies the optimal capturing of calcium ions in calcite pellets with a minimum achievable CCCP. The overall chemical yield is strongly dependent on the specific crystallisation surface area in the reactor [38]. For optimal process conditions, a large SSA in the fluidised bed is required. To be able to optimise pellet softening processes, the SSA has to be determined, which is possible in case the voidage and the particle dimensions in the whole fluidised bed are known [39]. Here, hydraulic process variables including superficial fluid velocity and viscosity, particle sizes, and densities are crucial. To obtain optimal process states of fluidised bed processes, effective voidage prediction models can be used in process automation and intelligent control [40,41]. Due to particle exchange in the pellet reactor, mainly caused by the extraction of calcite pellets and the dosage of seeding material, a certain particle size profile occurs over the reactor height. Due to the crushing process of calcite pellets [28], the seeding material consists of irregularly shaped particles. When they grow in size, they become increasingly round. If there is no difference in specific density, larger particles will migrate to the lower region of the reactor bed, and a stratified bed will evolve.
Since pellet softening is a continuous process, boundary conditions must be monitored on a regular basis. As a consequence of the principle of seeded crystallisation, the most decisive boundary condition is to maintain the fluidised state permanently. Therefore, the risk of fixed bed state, caused by calcite pellets that are too large or water flows that are too low, must be avoided at all times. For water suppliers using surface water, the temperature also has consequences for process control strategies. In addition, flushing of the smallest grains out of the reactor, mostly the smallest fractions of seeding material, must be prevented due to its effect on subsequent treatment processes. In fullscale pellet softening installations, the particle size mostly varies between 0.3-2.0 mm, and particle density between 2.5-4.0 kg/L. The current pellet size set-point in Amsterdam reactors is 1-1.2 mm. To retain fluidisation conditions, it is important that the largest pellets, usually those that are larger than the given set-point, are extracted from the reactor. They can then be used as a by-product in other processes or reused as seeding material.
In full-scale reactors, the fluid-particle characteristics, i.e. homogeneous and heterogeneous flow regimes [42], as well as many practical matters such as the fluid distribution through nozzles and the use of a bypass flow in operational state [43] determine the process state and quality of the fluidised bed. In addition, pellet softening can be seen as a combination of chemical and physical processes while the biological activity on the surface of the calcite pellets also affects the friction and the degree of bed expansion [44]. The combination of a large SSA aim, the level of particle profile distribution over the reactor height, the degree of irregularity and size of the growing calcite pellets and finally environmental and process state conditions makes optimal control of a pellet softening process a complex matter. Since the bed voidage is a critical process variable, the aim of this work is to develop a highly accurate voidage prediction model, as a function of fluid and particle properties which can be applied in full-scale pellet softening fluidised bed unit operations.
A computational fluid dynamics (CFD) approach for obtaining an effective voidage prediction model is not preferable [42], due to the above-mentioned non-ideal circumstances in combination with complex flow behaviour, numerous particle interactions [45], and the large amount of particles in full-scale industrial fluidised bed reactors, as is the case for water pellet-softening (up to 10,000,000,000 particles). This results in extremely high computational costs and a lower suitability for process optimisation and plant-wide control. To cope with constantly changing operational conditions in full-scale installations, more straightforward models are needed for optimal and robust process control. In particular, there is a need for an explicit and easily applicable voidage prediction model that effectively takes into account the local and global multiphase flow phenomena occurring in full-scale installations.
The aim of the current study was to develop a straightforward model based on the well-known dimensionless particle Reynolds and densimetric Froude numbers because they represent the inertial, viscous, and gravitational forces in the multiphase system. The new model should be applicable for a wide range of particles, sizes, fluid velocities, and temperatures and thus be of direct use in water treatment processes like drinking water pellet softening. For process control purposes, a direct relationship between particle size and voidage for a given water temperature and volumetric flowrate (default operational and design) is necessary.

Voidage prediction models
In a fluidised bed, the voidage, particle size, and physical properties of the fluid and particles are inter-related. In the literature, numerous multiphase models are given to predict the voidage in fluidised bed reactors, mainly aimed at gas-solid systems [46][47][48][49][50][51] but also for liquidsolid systems [52][53][54][55]39]. Specifically, for liquid-solid systems, the voidage can be predicted using classical models, such as the Richardson-Zaki [56] approach based on terminal settling velocity. The other frequently applied method is based on the idea of a flow through an assumed collection of channels in a bed of particles [57]. In the Ergun approach [58], the voidage is determined by the balance between the pressure gradient over the fluidised bed due to the mass of the pellets and the drag force of the water exerting on the pellets. The pressure gradient is accordingly given by the submerged weight of the pellets. The Carman-Kozeny model [59] is derived from a drag model where viscous and inertial forces are balanced using the modified particle Reynolds number. The impact of voidage is non-linear for both Kozeny-Carman and Ergun models [60] and more pronounced at lower voidage. Since in water treatment the operational field lies in the vicinity of incipient fluidisation and, in addition, turbulent flow regimes are exceptional, a very popular model adapted for the transitional flow regime is given by van Dijk et al. [22]. Foscolo [61] presented a similar tube flow approach, albeit with several model improvements. Foscolo considered a tortuosity factor, effective length of the fluid path, and a forced interpolation, to correct the limit of single-particle interaction to be obtained for dilute systems.
Another well-known relationship proposed by Wen and Yu [62,63] is based on the dependency of the drag coefficient on the dimensionless Reynolds number [64], and on the assumption of a common voidage function for the entire flow regime. With this model, the voidage can explicitly be calculated for given particle Reynolds and Archimedes numbers. Di Felice [52] investigated the voidage function and the dependency of the particle Reynolds number and proposed an improved Wen-Yu overall voidage relation. According to Akgiray [65], who provided an extended evaluation of expansion equations for fluidised solid-liquid systems, there is no general agreement regarding which equation is the most accurate. Akgiray proposed a voidage prediction model based on an improved drag relation. Kramer et al. [42] proposed an improved drag relation taking into consideration the fluidisation stability to cope with heterogeneity phenomena in liquid-solid fluidised beds and increasing the overall voidage prediction accuracy. The traditional drag relation based on the particle Reynolds number was extended with the particle Froude number. One of the most popular and frequently used models for describing homogeneous liquid-solid fluidised suspensions is the model developed by Richardson and Zaki [56]. The superficial fluid velocity and terminal settling velocity, together with an empirical index, enables determination of the fluid voidage in a Table 1 Voidage prediction models from the literature.

Model Equation a
Boundary conditions Eq. nr.  [40,67,26] was also based on the Richardson-Zaki principle and a fitted Schiller-Nauman equation [68] to determine the terminal settling velocities for calcite pellets. To be able to predict voidage in the proximity of minimum fluidisation conditions, either the minimum fluidisation velocity must be known or the Richardson-Zaki index must be very accurate. Therefore, the Richardson-Zaki model was extended [41] with proven hydraulics-based models. The minimum fluidisation velocity is acquired with the Carman-Kozeny model [59] where the terminal settling velocity is acquired using the Brown-Lawler model [69], an improved version of the well-known model developed by Schiller and Naumann [68]. In the literature, many velocity-voidage prediction models can be found, the most popular of which are given in Table 1.

Reactor performance indicators
The particle size and voidage over the reactor height must be known to be able to estimate the SSA. The following reactor performance indicators are defined: specific surface area based on the reactor volume and the specific surface area based on the water phase as well as the specific space velocity (SSV). The most frequently presented definition of SSA in the literature [1,20,39,52] of granular beds given as the total surface area of the particle material divided by the bed volume. A s r , represents the available area per m3 reactor volume for crystallisation. For monodisperse spherical particles, the SSA, based on the reactor volume A s r , , is given by: The SSA based on the water phase A s w , provides more adequate information for crystallisation of CaCO 3 on the available total particle surface. This performance indicator A s w , resembles the proper interaction area between water and reactive surface and is defined as: Mercer et al. [15] demonstrated that calcite seed crystals improved the removal of dissolved calcium during precipitative softening and that the optimal seed dose depended on the surface area available for crystal growth.
Van Schagen et al. [29] showed that the pellet size, and consequently the SSA, had a significant influence on performance with respect to the water quality parameter. To keep the super-saturation in the pellet softening reactor at acceptable levels, the SSA must be known to determine the crystallisation rate over the height of the reactor. Here, the voidage, particle size, and temperature are important [43]. Also, according to van der Veen et al. [23], the SSA is strongly dependent on the water temperature.
Time also plays an important role in the mixing zone of the reactor, where the majority of the crystallisation reaction takes place. The caustic soda, the water, and seeding material must all be mixed in a minimum period of time in order to prevent undesirable carry-over of CaCO 3 precipitation. Consequently, the water that passes along the available specific surface area per unit of time must be at a maximum level. The following performance indicator, specific space velocity, is derived from 'space velocity' as used in the field of chemical reactor engineering [70]. Space velocity for homogeneous reactions, i.e. uniform fluids, is defined as the number of reactor volumes of feed at specific conditions which can be treated per unit of time. Regarding multiphase systems, the reactor volume is partly occupied by grains. In that case, the effective water volume must be used similar to the Empty Bed Contact Time (EBCT), defined as the volume of the empty bed divided by the flowrate [1,20]. Considering that crystallisation occurs in multiphase systems, space velocity can be translated into the rate in which the water passes the SSA or how often the water is renewed at the water layer above the particle surface. Due to the presence of grains and the effective residence time, the superficial fluid velocity must be corrected with respect to the voidage. The specific space velocity A c is accordingly defined as the contact area per second per m 2 of transfer surface area:

Empirical data driven voidage prediction models
The fluidised bed voidage depends primarily on the fluid and particle properties. In general, one uses the mean particle size, assuming perfect spheres d p , average particle density p , kinematic viscosity T , and superficial fluid velocity v s . For the sake of simplicity, wall effect corrections, the influence of the fluid distributor, fluid circulations, and the irregular distribution of the particles and other non-ideal phenomena are ignored.
Voidage prediction models are only valid for a fluidised state. For this reason, it is important to determine the incipient fluidisation and maximum flushing point to check the prevailing state. Detailed elucidation and model derivations are presented in the Supplementary Material (Chapter 10, 11 and 12).

Voidage prediction polynomials (VPP)
A straightforward way to model the voidage for individual grain types is to use sets of polynomials as a function of velocity, viscosity, and particle size, respectively, while considering the particle density to be constant, either indirectly or directly related to particle size. There are two ways to apply these polynomials. Both methods are based on an experimentally obtained dataset containing a wide range of data in fluid velocity covering the whole temperature regime and various particle diameters as well as a given particle density. A 'floating polynomial model' is given in the Supplementary Material (Chapter 11). However, to keep the fitting procedure feasible for processing, a second approach using a one fit polynomial equation is proposed by means of Eq. (15): The advantage of this approach lies in the simplicity of having one single model for which fitting parameters can be obtained using nonlinear regression software. The disadvantage, however, lies in the considerable computational time required for finding an acceptable prediction accuracy and the number of fitting parameters = N l m n ( ). The computational time increases substantially with the size of the experimental dataset. Using fit polynomials works well, although it demands strictly respecting the given boundary conditions to avoid overshoot and physically unrealistic voidage predictions.

Dimensionless numbers application (DNA)
The effective voidage can also be predicted using dimensionless Reynolds and Froude numbers based on an implicit drag relationship proposed by Kramer [42]. The Reynolds number deals with the relationship between viscous and inertial forces and determines the degree of laminar or turbulent flow regime [71]. The particle Reynolds number Re p is defined as: The Froude number is defined as the ratio of inertial to gravity forces and is a proxy for the fluidisation quality from smooth homogeneous (particulate) fluidisation to heterogeneous or aggregative (bubbling) fluidisation [72]. The densimetric particle Froude number Fr p is given by [73]: To find the voidage from these dimensionless parameters, a numerical method is needed.
This equation can be rewritten as a straightforward function of velocity, viscosity, particle size and density. More information is presented in the Supplementary Material (Chapter 18 and 19).

Double Reynolds-Froude model (Rep2Frp).
Eq. (18) is a rather straightforward approach which can be improved based on previous research in the field of chemical engineering [74]. The general Reynolds number (Re) determines whether the flow is dominated by inertial or viscous forces, i.e. whether the flow is laminar or turbulent. Whitaker [75] proposed a heat transfer-based equation for the Nusselt number Nu as a function of the Reynolds number Re and Prandtl number Pr, for flow in pipes, around spheres, and through packed beds as follows: Note that Nusselt numbers Nu between 1 and 10 are characteristic of laminar flow, while turbulent flow typically corresponds to Nu in the range 100 -1000 [76]. Therefore, it can be understood qualitatively that Whitaker's expression for Nusselt increases with increasing Reynolds number and captures the influence of flow regimes on heat transfer [75]. This understanding is supported by Bedingfield and Drew [77,39], who showed a theoretical analogy between heat and mass transfer. Similarly, for liquid-solid fluidisation, we will formulate an expression for the voidage based on the sum of two terms that capture the influence of the flow regimes.
Here we use the Froude number instead of the Prandtl number. The reason for implementing the densimetric Froude number comes from the hypothesis [42] that voidage is based on laminar-turbulent flow regimes as well as heterogeneity phenomena in liquid-solid fluidised beds. The voidage can accordingly be predicted using Eq.

Symbolic regression model (SRM)
Based on high quality datasets, highly accurate prediction models can be obtained using symbolic regression techniques as applied in genetic programming. Genetic programming is a random-based technique [78] for automatically learning computer programmes based on artificial evolution. It has been successfully used in many applications [79,80]. The advantage of genetic programming is that there is no need to define the structure of a model a priori: the technique randomly generates a population of several mathematical operators. Symbolic regression is the process of determining the symbolic function, which describes a dataset and effectively develops an empirical model [81]. These types of models have two main features: complexity and accuracy. Generally, given a certain dataset, the process starts with the determination of very simple but inaccurate models. With time, more accurate but also more complex models are obtained. To prevent adverse modelling of measurement errors, data noise or deviation, a model should be taken as a compromise between complexity and accuracy. Voidage is a function of fluid velocity, viscosity, and particle size and density given in a general form according to Eq. (21):

Experimental setup
To calibrate and to validate the prediction models, liquid-solid expansion experiments were needed to obtain reliable datasets containing fluid viscosity and superficial velocity as well as particle size and density. Advanced laboratory and pilot plant apparatus were therefore especially designed for this purpose. In addition, drinking water related grains were carefully prepared and selected.
Detailed information about experimental expansion columns can be found in the Supplementary Material (Chapter 5), particle selection (Chapter 3), particle and fluid characterisation (Chapter 8) and fluidisation expansion experiments (Chapter 9).

Particle selection, properties, experimental setup and fluidisation experiments
Liquid-solid expansion experiments were carried out at three locations: in Waternet's Weesperkarspel drinking water pilot plant located in Amsterdam, the Netherlands; at the University of Applied Sciences Utrecht, the Netherlands; and at Queen Mary University of London, United Kingdom.
In this study, we examined two kinds of particles, calcite pellets (100 % CaCO 3 ), and crushed calcite seeding material grains [28], both applied in drinking water softening [11]. Polydisperse calcite pellets were sieved and separated in order to acquire more uniformly dispersed samples.
The morphological particle properties obtained with a Camsizer [82] show that crushed calcite and the smallest fractionised calcite pellets have irregular shapes. The larger the grains become, the more spherically shaped they appear to be. Photographs of these particles can be seen in the Supplementary Material (Chapter 2 and 4).
The acquired experimental dataset consisted of a matrix with varying temperatures, grain sizes, and flowrates, as was required for a comparison of the theoretical fluidisation models. In total, 61 liquidsolid fluidisation experiments were carried out for a wide range of calcite pellets (0.425 < d p [mm] < 2.8), and a total of 42 experiments were carried out for crushed calcite (0.4 < d p [mm] < 1.12), summarised in Table 2.
Additionally, in total, 89 additional independent liquid-solid For validation purposes, additional liquid-solid fluidisation experiments for a wide range of different particles in fluid water systems were conducted at the three given locations.
Detailed information regarding particle and fluid characterisation, standard operating procedure of the fluidisation expansion experiments, photos of grains, and data tabulation can be found in [41] and in the Supplementary Material.

Voidage prediction
Based on the expansion column experimental datasets from our experiments, several empirical data driven models were derived using symbolic regression. Model parameters were found through non-linear curve fitting.
For implicit models, the voidage prediction accuracy was found with a straightforward Bolzano's numerical intermediate value theorem. All equations, regression coefficients, fitting parameters and plots are given in the Supplementary Material. The experiments cover a large range of variables, as given in Table 2, which is assumed to represent the boundary conditions for the presented models.

Voidage prediction polynomials (VPP)
Two different polynomials were tested. The first one is a triple quadratic polynomial (222), for which the voidage can be explicitly estimated based on velocity, viscosity, and particle size with a constant assumed particle density. The second one is a 4th, 2nd, and 3rd order polynomial (423) versus velocity, viscosity, and particle size, also with a constant particle density. These polynomial expressions are relatively accurate for the voidage prediction in case the given boundary conditions are respected. However, when these boundary conditions are exceeded, these VPP models are unpredictable and rather inaccurate. Another disadvantage is the high computational effort required for the non-linear curve fitting process with the aim to find all the fitting parameters. (Rep1Frp). The fitting parameters for Eq. (18), describing a combination of laminar versus turbulent flow regimes and homogeneous versus heterogeneous flow characteristics, are given in Table 3.

Double Reynolds-Froude model (Rep2Frp).
The fitting parameters for the improved Reynolds-Froude based Eq. (20) are given in Table 4.
Model adjustments for process automation purposes without verifying boundary conditions are described in the Supplementary Material.

Symbolic regression model (SRM)
Based on the experimental expanded bed dataset, numerous solutions, i.e. multiple equations, were found when the software package Eureqa [83] was used. We present two examples of equations: one for calcite pellets in Eq. (22) and one for crushed calcite in Eq. (23). Fitting parameters can be found in Table 5. The expressions are accurate within the given boundary conditions. However, as is the case with the polynomials, they are unpredictable and therefore inaccurate when the boundary conditions are violated.

Voidage-velocity graphs
In Figs. 1 and 2 we present our experimental voidage-velocity data in comparison with two models, respectively: the Richardson-Zaki equation (8), as one of the most used and well-known models; and the new Reynolds-Froude Rep2Frp Eq. (20) proposed in this study. Graphical results for all voidage prediction models are plotted and shown in the Supplementary Material.
The voidage prediction models (Table 6) given in the literature (N = 11) as well as the models proposed in the current work (N = 5) were compared with the experimentally obtained data using five statistical methods: mean average error, average relative error, normalized root mean square error, logarithmic root mean squared error, and Pearson's correlation coefficient.

Reactor performance indicators
The SSA is a relevant performance indicator concerning pellet softening and can be estimated in case the voidage and particle size are known. The superficial fluid velocity determines for a great deal to what extent the bed expands. Based on Eq. (12), describing the bestknown SSA for the reactor volume, a 3D plot can be created. Fig. 3 represents the SSA A s r , against the pellet size d p and the linear flowrate v s , for a given water temperature. However, the operation window for full-scale pellet softening reactors is smaller (60 < v s [m/h] < 120): this is plotted in Fig. 4. A surface plot, using Rep2Frp model (Eq. (20)), is accordingly plotted in Fig. 5.
In addition, the voidage is plotted against the pellet size d p and the linear flowrate v s in Fig. 6. Another less familiar reactor performance indicator is the contact surface area per unit of time (Eq. (14)) which is plotted Fig. 7.

Graphical exploration
The prediction models presented in this work, i.e. voidage prediction polynomials, dimensionless number applications, and symbolic regression models, were compared with the most popular and familiar models known from the literature ( Table 1). All voidage-velocity plots are included in the Supplementary Material.
Before assessing the models in terms of statistical fit quality, we applied graphical exploration, as proposed by Anscombe [84], to test the goodness of fit and to determine whether the model describes the experimental data adequately. Fig. 2, for instance, shows both dots (experimental data) and lines (prediction model) that provide a visual image of the degree to which the dots and lines coincide, where the discrepancies appear, whether or not the dots are over-dispersed, and whether the boundary conditions are violated and thus lead to unrealistic values.
Based on graphical explorations (Figs. 1, 2, and remaining figures in the Supplementary Material), the Ergun model [58] shows Table 3 Fitting parameters for Eqs. (18)  overprediction, in particular for higher velocities and larger grains. The Carman-Kozeny model [59] shows overprediction for smaller grains at low and intermediate velocities.
Van Dijk's model [22] is originally fitted for a small operation window [43,26,40,85] and therefore outliers can be found in and outside this specific region. The prediction quality of Foscolo's model [61] is relatively good at lower velocities but underpredicts in case the velocities are increased. Akgiray's model [65], on the other hand, works increasingly better for higher velocities. The Wen-Yu model [62] shows a good fit at low velocities and small grains but works less well for higher velocities and larger grains. The model modification proposed by Di Felice [52] did not improve the fit quality. The RIO 2 model [42] appears to work well in terms of fit quality, not including the smallest grains at maximum velocity. The Richardson-Zaki model [56] in general underestimates the experimental data considerably ( Fig. 1) and, to a lesser extent, so does the van Schagen model [26]. The other RZ-based model proposed by Kramer et al. [41] starts to show discrepancies for higher velocities and larger grains. Regarding the models presented in this study, the polynomial curves coincide with the dots except for the smallest grains and show unrealistic curvature caused by the fitting towards the maximum velocities. The empirical symbolic regression model shows appropriate slopes, albeit with a larger offset for increased velocities. The single Reynolds Froude based model is not accurate for small grains at low velocities while the double Reynolds-Froude model shows good agreement with the experimental data.

Statistical exploration
Numerical voidage prediction accuracies are given in Table 6 for five statistical metrics (MAE, ARE, NRMSE, LRMSE and R 2 ). In general, the calculated errors of all examined models, based on experimental calcite pellets data, is a 12 % or lower. The prediction errors of 11 of our 16 models are lower than 4 %, and 5 out of 11 are lower than 2 %. The following models have the lowest prediction error: Eureqa, Rep2Frp and the RIO 2 model, with the exception of R 2 . Based on an independent dataset, given in the Supplementary Material, the prediction errors are slightly higher (≈1 %) and show the same top-three ranking. The differences in accuracy are mainly based on the use of older, and less accurate, versions of the experimental setup and measurements devices.
Based on Pearson's correlation coefficient R 2 , it is impossible to make a well-argued choice for a preferable model or to determine which model is the most accurate. Almost all R 2 values are higher than 0.99, and in most cases, they are 0.999. The MAE measures the average magnitude of the errors in a set of predictions, regardless of their direction and where all individual differences have equal weight. The NRMSE gives a relatively high weight to large errors, and therefore the RMSE is more useful when large errors are particularly undesirable. In case of LRMSE, the outliers are drastically scaled down, thus nullifying their effect. Since we are looking for a very accurate voidage prediction model, there is a slight preference for NRMSE as it penalises undesirably large errors.

Richardson-Zaki based models
The Richardson-Zaki model, which is still very popular in the literature, shows the highest prediction error. An explanation can be addressed to the terminal settling velocity, at an assumed voidage = 1, whether the observed voidage at maximum or entrainment state was lower: approximately 0.95. The RZ model uses terminal settling as a starting point and combined with an empirical index, which is also a function of terminal settling, the voidage can be predicted. When a voidage lies in the vicinity of incipient fluidisation, as is the case in pellet softening processes, a relatively small inaccuracy in both the starting point and the index n in equation (8) leads to a high prediction error. When log-log scales are used, these effects are clearly visible (Fig. 1). Van Schagen solved this terminal settling issue through finding better fitting parameters for the Schiller-Naumann [68] equation. He extrapolated expansion data for calcite pellets to terminal settling conditions. With this adjustment, the error was reduced by a factor of almost two. Kramer et al. [41] introduced a second physical point. Besides the point of terminal settling, an extra point of minimum fluidisation was added from which the voidage can be interpolated. The prediction error using these model improvements also reduced with almost a factor of three but still shows some overprediction.

Porous-media-based drag relations
Based on statistical metrics only, it is not evident that the Ergun model is less accurate for the higher flow regime. The model has, however, shown a reasonable fit score at low flow, which is confirmed by Marshall [86]. Based on statistics, the Carman-Kozeny model has a slightly higher fit score, but here too it remains unclear whether this model is actually more suitable in the transitional flow regime. The CK model overpredicts at low flow regimes. The models introduced by van Dijk, Foscolo, Wen-Yu, and Di Felice have acceptable fit scores, but this is mainly at low flow regimes. For pellet softening, voidage is overestimated by van Dijk, but underestimated by Foscolo and Di Felice. Akgiray, Wen-Yu, and RIO 2 show conformity with the data.

Data driven and numerical numbers-based models
The polynomial functions underpredict the voidage for the pellet softening operational range, whereas the empirical symbolic regression model is less accurate at the borders, which is also the case for the single Reynolds-Froude based model. The double Reynolds-Froude based model has a high fit quality and is more reliable since the model is based on dimensionless numbers with a hydraulic physical significance. This implies that, in case the boundary conditions were to be violated, the risk of run-away is smaller, compared to purely empirical data driven models.  (22) and (23).

Modelling aspects
Several factors determine, to different levels and certain degrees, the fit prediction quality: these include non-ideal aspects such as particle polydispersity, morphological properties, density differences, and particle interactions. Non-ideal circumstances such as fluid characteristics, heterogeneous flow phenomena, and the influence of chemicals also play a role. Other determining factors include non-ideal matters related to the piece of apparatus used, like flow distributor, column alignment, sensor inaccuracies, and wall effects. The effort to incorporate all these aspects into a model is significant, but it also makes the model more complex. In full-scale, fully automated unit operations, a reliable, explicit, and programmable model is preferred -despite the penalty it might bring in case the fit quality is slightly lower. The models can be fitted accurately when high-quality experimental data are used, based on calcite pellets extracted from the full-scale reactors.

Bed state control
To reach high performance levels in full-scale pellet softening reactors, an optimal crystallisation process is important, but this is strongly dependent on the governing hydraulic state in fluidised beds. Pellet softening in drinking water production processes is a continuous process. Discharging calcite pellets from reactors and subsequent dosing fresh seeding material are likewise continuous repetitive processes. These particle changes imply that the SSA also varies and will drift away from its ideal setpoint.
At the bottom of the reactor, the voidage is kept relatively low to obtain the highest crystallisation SSA; but nevertheless, fixed bed situations must be avoided. The degree of voidage is dependent on the physical properties of the grains and the water viscosity. De facto, voidage, or fluid bed height, is kept constant through controlling the water flow in the reactor and, depending on water temperature, through particle bed management. In pellet-softening reactors, voidage is approximately ≈ 0.55 at the bottom of the reactor and ≈ 0.8 at the top.
Full-scale pellet softening reactors are always installed in groups, and hard water containing a high concentration of calcium ions is often partly bypassed [43]. Moreover, the total water production changes in volume periodically. In a full-scale operational pellet softening reactor, the process state is subject to changes in water flow, temperature, depending on the season of the year, and ongoing variations in particle sizes, shaping the particle profile over the height of the reactor bed. To cope with these changes, constant process state monitoring is crucial. To control the required pellet sizes, as a standard conventional procedure, particle samples are regularly taken from the reactor bed manually and accordingly analysed in a laboratory. In that case, based on water flow, for a given temperature, together with the particle size, the voidage can be estimated. This is demonstrated in Fig. 8. Engineers can check the quality of the SSA but also monitor the risk of a fixed bed state or flushing.
In industry, soft sensors are also implemented [40] with gauges at the bottom of the reactor in which the voidage and particle size is derived from flow, differential pressure, and temperature. This is done by using models. The advantage of this method is the availability of on-line data results. However, a disadvantage is the vulnerability and consequently sensitivity of the gauges due to the exposure to high lime scaling conditions in the lower zone of the reactor, resulting in less accurate predictions.

Reactor performance indicators
The SSA plotted in Fig. 3 clearly shows that the A s r , decreases when the linear flowrate increases or particles become larger. This is similar for pellet softening conditions (Figs. 4 and 5). This indicates why smaller particles, like seeding material, increase the SSA. Still, these crushed calcite grains however, migrate to the higher zones of the fluidised bed due to stratification and, in case the flow is too high, the risk of flushing out of the reactor emerges. At the same time, due to CaCO 3 crystallisation, particles grow and migrate to the lower zone of the reactor, where the chemical driving force is large and decisive. These larger grains cause the SSA to decrease fast and, in addition, they Voidage is plotted against velocity. Dots represent experimental data and curves represent the actual prediction model.

Table 6
Model overview voyage prediction accuracy calcite pellets arranged in descending order by accuracy: mean average error (MAE), average relative error (ARE), normalized root mean square error (NRMSE), logarithmic root mean squared error (LRMSE), Pearson's correlation coefficient (R 2 ). Model sources are given in Table 1 enlarge the risk of a fixed bed state, which can be seen in Fig. 6 as indicated by the dark blue zone. To maintain a maximum SSA, large calcite pellets should continuously be withdrawn and replaced by smaller crushed calcite seeds, which results in higher operational costs, mainly caused by transportation, and also adversely affects the sustainability goals. Therefore, an optimum fluid velocity must be chosen very precisely. This substantiates the relevance of a very accurate voidage prediction model. In Fig. 7 a maximum specific space velocity A c is reached for calcite pellets, with an average size slightly larger than 1 mm. With Eqs. (12) and (14), an optimal linear flowrate range (75 < v s [m/h] < 85) can be derived for an average water temperature, with a corresponding SSA: A s r , ≈2500 [m 2 /m 3 ]. This is partly in agreement with the current operational window of the Amsterdam reactors at Waternet [38,17,43,87,88]. Another complexity arises when the ratio of calcite pellet size to crushed calcite size is lower than would approximately be the case in a process state where the voidage is also in the vicinity of incipient fluidisation. In case the voidage is too low, caused by either lower flowrate or higher water temperature, dosed crushed calcite remains trapped between the calcite pellets, leading to a non-stratified particle bed.
In full-scale installations with plant-wide control, complex models and the continuous challenge for finding optimal numerical solutions are less desirable. Another disadvantage is that many models are semiempirical and derived for monodispersed perfectly spherical particles.
In de facto all full-scale multiphase flow processes applied in water treatment processes, the particles are irregularly shaped and often highly polydisperse. Due to the complex flow behaviour and large amount of particles (N∼10 10 ), explicit particulate modelling of fullscale industrial fluidised bed system using Computation Fluid Dynamics is also challenging and as yet unachievable.
Therefore, there is a need for an explicit effective model to be able to accurately predict the overall voidage in fluidised beds that effectively takes into account the global multiphase flow phenomena and many other non-ideal matters occurring in full-scale installations. Based on these criteria and findings, an average optimal linear flowrate can be determined (v s ≈ 85 ± 5 [m/h]). Due to the above-mentioned nonideal circumstances, full-scale operational challenges and continuous changes in particle profile, further research is needed to find a more precise optimal process state that also takes into account hydraulic, chemical, and biological phenomena.

Conclusions
The accurate calculation of voidage and specific surface area is of major importance in drinking water treatment processes like pellet softening, because it determines the process conditions and treatment results.
To maintain or provide optimal process conditions in pellet-softening reactors, it is important to accurately determine the fluidised bed voidage. Voidage is a crucial variable for determining the specific   surface area, the minimum fluidisation and flushing conditions as well as the water and particle residence time. The voidage prediction accuracy of 11 models from the literature was compared with five types of data driven and dimensionless numbers-based models. Accurate experimental datasets were used to validate the predictive power. Porous-media-based drag relations models RIO  Carman-Kozeny) have the disadvantage that the voidage must be numerically solved. The polynomials, Wen-Yu, and the Reynolds-Froude based models provide explicit solutions for the voidage. As mentioned above, in full-scale installations with plant-wide control, complex models and the continuous challenge for finding optimal numerical solutions are less desirable. Therefore, there is a need for an explicit effective model to be able to accurately predict the overall voidage in fluidised beds that effectively takes into account the global multiphase flow phenomena and many other non-ideal matters occurring in full-scale installations.
As a consequence and based on statistical metrics as well as graphical exploration, a preferable prediction model for a pellet softening operational range (60 < v s [m/h] < 120) is the Rep2Frp model. The model is accurate, has a physical basis and is easy to use for process automation programmes and therefore suitable for application in fullscale operational pellet-softening reactors.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.