Modelling of Hydraulic Dynamics in Sludge Treatment Reed Beds with Moving Boundary Condition

. The conventional method of simulation using fixed mesh method (FMM) of discretization is a well-known and trusted procedure in modelling hydraulic dynamics. However, new ideas of innovation in modelling should be advanced. The moving mesh method (MMM) has been considered as a novel approach in modelling hydraulic dynamics after depending on the existing simulation model for decades. The MMM is capable of describing the moving boundary condition of an actual wetland system due to water ponding. An idealized model should be able to simulate the actual hydraulic flows through the system with the corresponding porosity. Hence, a combination of MMM and FMM (MM-FMM) of discretization for hydraulic dynamics is studied in this project to model the flux with respect to water ponding scenario in a sludge treatment reed bed and unsaturated transient flow within the bed. Such method has evidently proved to simulate the actual hydraulic flows in contrast to conventional method. The application of MMM limits the maximum flux to keep within its saturated conductivity, thus reduces the effect of flow overprediction. Subsequently, the simulated results for hydraulic head and moisture content can be predicted for actual condition of different cases according to their respective fluxes.


Introduction
Implementing constructed wetland (CW) or sludge treatment reed bed (STRB) for sludge dewatering and mineralization has become common in developed and developing countries. STRB uses physical processes to dewater the sludge. The sludge is fed intermittently on top of the wetlands, allowing the wastewater to percolate downwards through the bed, where a multi-layered granular substrate filter is usually found. Thus, the particulate contaminants, which consist of organic and inorganic matter, are physically retained on the top surface of the bed, subsequently forming a layer of sludge deposit. At the same time, the treated liquid waste is discharged from the bottom of the bed [1]. Although the capital cost for a conventional STRB is often higher than a mechanical dewatering device due to extensive land requirement, they offer less energy, are chemical-free, and produce a relatively high-value biosolid that can be used for agricultural purposes. Despite the promising advantages of STRB, experimental studies are still ongoing to determine the optimum operating parameters, including the loading rates and resting periods [2][3][4].
Nonetheless, experimental studies are always restricted by high capital expenditures and time consumption. Thus, several process-based models have been formulated as simulation tools to provide an alternative to studying the system design. There are limited process-based models for STRB, as most existing models were developed to describe wastewater treatment in CWs. For instance, HYDRUS-CW2D [5], HYDRUS-CWM1 [6], and CFD [7] are the most commonly used models for a similar system. The numerical models used in CWs have described the wastewater treatment processes by integrating several sub-models, including hydraulic, reactive transport, plants, biochemical reaction, and substrate clogging [8]. The hydraulic sub-model described the dynamics of water flow, where Richards' equation (RE) has been well-known for predicting the hydraulic flow in the porous medium, yet it has successfully predicted the hydraulic behaviour in similar systems [9].
However, these models assumed the water ponding level as a fixed layer, resulting in overprediction of the dewatering capacity. This assumption is valid when the water level fluctuation is negligibly small, but it causes a critical problem in simulation when the fluctuation is relatively significant [10]. In this scenario, the top boundary of the water ponding is assumed to be in the fixedmesh form, where the water level is constant throughout the simulation. In addition, the build-up of sludge deposits on the top surface of STRB also changes the simulation mesh [11]. Therefore, researchers have developed the adaptive meshes refinement (AMR) scheme to improve the conventional fixed meshes scheme [12]. Upon the presence of AMR, they have successfully predicted the shock hydrodynamics via moving mesh. Results have shown that the adaptive meshes produce better results than the fixed meshes when only a fraction of the domain in the areas of interest is being investigated. Moreover, it is computationally less expensive as a shorter simulation time is needed when the resolution in these areas is increased instead of performing the mesh refinement over the entire grid [13]. Therefore, this study aims to develop a process-based model with adaptive meshes to describe the hydraulic performance of STRB using conventional RE. The main objective of this study is to formulate a robust model in simulating the hydraulic flows through the STRB with respect to the sludge ponding on top of the bed and sludge deposition.
Furthermore, the application of adaptive meshes refinement (AMR) scheme in the wastewater treatment has yet to be generalized. The conventional method of discretization, as known as fixed mesh method by Picard iteration is indeed trustworthy in predicting hydraulic flows but possesses possible drawbacks, perhaps flow overprediction. A combination of moving mesh method (MMM) and fixed mesh method (FMM) of refinement in simulating the hydraulic flow in the porous medium is capable of illustrating the actual dewatering condition in STRB, where the infiltration of sludge above the reed bed surface is described by the MMM while the remaining gravel layers are always static, thus discretized by the FMM. Conclusively, a combination of moving mesh and fixed mesh method of simulation is a novel approach to the STRB in wastewater treatment.

Methodology
In a conventional model simulation, a mesh equation and the differential equation are often solved together simultaneously to generate the new nodes and solution. However, a moving-mesh method (MMM) is an alternative solution for the moving-mesh scenario, in which the mesh and solution are varied simultaneously in such a way that a fixed number of nodes remain concentrated in regions of rapidly variated solution [13,14]. MMM is one of the most common adaptive mesh refinement methods. By implementing this method, there is no longer a need to interpolate dependent variables from the previous values in the mesh, which provides flexibility in describing the boudnary changes in the simulation. MMM is widely used with the model using the finite difference (FD) method as the solution technique in a one-dimensional domain.
For the discretization of the fixed mesh method (FMM), the explicit scheme is the simplest way to solve the partial differential equation. In the discretization of the Richards' Equation (RE) in the hydraulic model, the first derivative to time is estimated by the first order forward difference and the second derivative to space is determined by a second forward difference. The Picard method is a wellknown iterative procedure adapted for solving the resulting discretized non-linear equations [15]. This study aims to integrate the fixed-mesh method with the moving-mesh approach to numerically simulate the long-term hydraulic performance of a STRB.

Governing Equations
Richards' equation (RE) is commonly implemented to describe the transient water flow processes in the porous medium under unsaturated condition. RE can be expressed in three forms, namely h-based, -based, and mixed form, where ℎ [ ] stands for pressure head, [ 3 / 3 ] stands for moisture content, and mixed form includes both dependent variables. A mixed form of RE is presented below [16]: where ( In the RE, the van Genuchten-Mualem model is always used to determine the volumetric water content (ℎ) and hydraulic conductivity (ℎ) [17], which are given as follows: where and are the residual and saturated water contents, respectively [ 3 / 3 ]; [ −1 ] is the inverse of the air-entry pressure; is a measure of the pore-size distribution; = 1 − 1/ ; is the saturated hydraulic conductivity [ / ]; and = 1/2 is an empirical pore-connectivity parameter. In summary, , , , , , and are the function of soil hydraulic curve and the values must satisfy the principles of 0 ≤ < , > 0, and > 0.

Moving mesh RE
Conventionally, FMM of discretization is often obtained by the backward Euler approach, which is then transformed into a tridiagonal non-linear set of equations, and eventually solved by the tridiagonal matrix algorithm [18]. However, to implement the moving mesh method in the RE, the main concept is to keep the fractional amount of water between adjacent mesh points constant [19]. Thus, the formulation of the conservation equation is given by: where is constant in time.
In the development of moving mesh RE, Equation (4) is differentiated using Quotient rule with respect to time, and substituted back the constant, : Meanwhile, the integral in Equation (4) is integrated by using Leibniz rule: Since this study is dealing with mixed form of RE, hence: Furthermore, the intention to transform the equation into velocity-based, the equation is rearranged by expanding the integrals: Furthermore, the RE is formulated into a velocity-based equation before being discretized using the MMM with the FD approach. The schematic illustration of a vertical flow constructed wetland (VFCW) under moving boundary condition is displayed in Figure 1. Thus, the velocity based RE is displayed as:

Boundary Conditions
Every STRB has its own filtration capacity to receive loadings, which depends on the permeability of the substrate used. When the hydraulic loading rates exceed the infiltration capacity, the top surface layer of the bed is completely saturated. In this case, the excessive loadings accumulate on top of the bed, resulting in temporary ponding. In order to describe the loading of the bed, both Dirichlet or the first type (constant ℎ or ) and Neumann or the second type (constant flux, ) boundary conditions are used in the model [20]. According to the illustration shown in Figure 1, it can be deduced that the top boundary condition, are set to be switching in between first and second types of boundary conditions due to the intermittent loading mode, whereas the bottom boundary, is a free drainage condition. Hence, the upper boundary conditions are as follows: where is the liquid flux at the medium surface, 0 is the prescribed liquid flux input, 0 is the initial time, is the end of the loading period, and is the total simulation time. In contrast, in consideration of no imposed flux and free drainage at the bottom of the bed, the lower boundary is then set to be a zero-gradient flux boundary condition, which expressed as: Therefore, the discretized velocity-based moving boundary RE with consideration of the mentioned boundary conditions, the final equation is given as: where = 1, 2, 3, … , − 1, ∆ ̅̅̅ = (∆ + ∆ +1 )/2, +1/2 ̅̅̅̅̅̅̅̅ = ( + +1 )/( + ), = ∆ / 2, = ∆ +1 /2, ∆ = − −1 , and ∆ +1 = +1 − .

Preliminary Results and Discussion
MATLAB ® Simulink R2022b is used in this study for the simulation. As discussed in the previous Section 2.2, water ponding is the main reason for causing the variation of the top boundary layer. Hence, in this study, the sludge ponding is considered a moving boundary system, while the substrate filter is a fixed system. By combining moving and fixed mesh methods (MM-FMM), the hydraulic flux, water content, and pressure head profiles showed a better simulation result than the conventional method of mesh discretization. Figure 2 compares the preliminary results between MM-FMM and FMM of simulation. A mesh increment of 1 cm height was used to simulate the actual wetland system of 60cm height over 600 mins. Temporal discretization of 1/360 minutes was applied to improve the accuracy of the model. Initial profile of pressure head, h was assumed to be identically -12 cm, with the minimum and maximum tolerance of 0.01 and 1.00 cm, respectively.
Initially, a 60-cm mesh profile of volumetric water content obtained from the RE and pressure head are input into the model to initiate the simulation. The top 5 cm meshes of the system indicate the sludge deposits layer, followed by 3 different substrate layers (5,22, 27 cm) with respective gravel sizes (small>medium>large). Hence, only the top layer of meshes is expected to move according to the sludge ponding in MM-FMM simulation, while the other layers remained static. In fact, the first 5 cm meshes of MM-FMM also represent the moving-boundary condition due to the ponding, as it does not undergo any infiltration on top of the substrate filter. However, infiltration is always happening across the wetland bed in FMM.
At the beginning of infiltration, the water diffusing downward with an increasing trend as the bed is well drained and dry, but it gets saturated over time. Under ponding, the infiltration is at its maximum, where the water flows with a constant flux of 0.01 cm/min above the substrate layers in MM-FMM, which is the set value of saturated hydraulic conductivity, as shown in Figure 2a and 2b. In FMM, the maximum flux has exceeded the saturated amount as it does not consider the moving boundary condition due to ponding. Consequently, the simulated fluxes in the bed are overpredicted, thus affecting the model's actual performance.
The pressure head is showing a similar trend as illustrated in Figure 3a and 3b, where the pressure has started to drop at 300 mins in FMM. The pressure is the lowest when it comes to the surface of the substrates due to the porosity of the bed, as well as the influence of the sludge deposits layer. The development of sludge deposit is the main reason for the reduction of hydraulic flux and sudden drop in pressure within the bed. However, the pressure increases vertically going down the bed due to water saturation. The surface pressure only drops to negative when there is no more ponding.
Furthermore, the water content in the wetland bed has a decreasing trend from top to bottom across the simulation profile, as displayed in Figure 4a. The water content is the summation of volumetric/mobile and immobile water contents. It shows the highest value of 0.5 at the top of the sludge deposits layer in MM-FMM, resulting from the ponding scenario. However, the overall water content in the sludge deposits layer of FMM is less than MM-FMM, as some of the water is predicted to have percolated downward through the substrate filter rapidly, as indicated in Figure 4b. The increase in the thickness of the sludge deposits layer is mainly due to the sedimentation of particulate contaminants. The deposition of solid sludge on top of the bed has led to an increment of meshes, where the length of is no longer equivalent to one. Figure 5 presents the top mesh of the , where it gives a significant increment at the beginning due to both sludge deposits and ponding, then increases gradually based on the sludge accumulation, and decreases significantly at the end when there is no more ponding. The second to fifth meshes in the sludge deposits layer are maintained at approximately 1.11 cm due to the water ponding, whereas the other meshes of the substrate's layer are always 1 cm.
The flux simulation by FMM can predict the ideal cases of flow without ponding condition. However, the simulated results are roughly doubled of the measured data in ponded case. The results deviation could be insignificant for higher fluxes, but it shows a huge difference when the flux is small. A comparison among the simulated fluxes for FMM [21], MM-FMM, and measured data is shown in Figure 6. The overprediction of the flow was mainly due to the assumption that the surface above the substrate and sludge deposit layers is always remained stationary. Thus, the acting force of hydraulic pressure decreases from the very first mesh on top of the reed bed surface over simulation time.
In contrast, the simulated flux by MM-FMM assumes the meshes above the reed bed surface are always moving due to boundary changes. The further decreases of hydraulic head due to fall in sludge ponding level, causes the reduction in hydraulic flow over simulation time. Hence, the predicted consecutive flows across the substrate filters are much smaller compared to the FMM.

Conclusion and Future Recommendation
The moving mesh method (MMM) is a novel approach to the simulation of sludge treatment. In contrast to conventional fixed mesh method (FMM), the moving boundary conditions due to the ponding and sludge deposition can be well addressed via the MMM. The combination of MMM and FMM (MM-FMM) provides a better presentation of hydraulic dynamics in the STRB, which considered the actual ponding scenario, yet implements the conventional simulation of unsaturated transient flow in the porous substrate medium. Moreover, the relative pressure, water content, and porosity of the bed in the case of ponding can be ideally estimated as well. The study of model simulation using MMM and MM-FMM of discretization should not be ceased here, as it definitely provides a better result compared to the FMM of simulation.
In recommendation, a thorough research on sludge accumulation due to the continuous sedimentation and deposition on the upper layer of substrates should be conducted. The increment of sludge deposits thickness has definitely reduced the sludge infiltration rate, that has directly decreased the treatment performance. In fact, the permeability and conductivity of the sludge deposits as well as the substrate medium have reduced a significant value. Some of the solid sludge is confirmed to be percolated into the substrate filter, thus led to the reduction in hydraulic flow. Therefore, a continuous and advanced study of sludge deposit incrementation will certainly improve the system performance in wastewater treatment.
Moreover, the crack in sludge deposits has a major impact on the overall system performance. The sludge deposits are crucial to retain the influent raw sludge to increase the overall hydraulic retention time, as to enhance the removal efficiency of the percolation performance. However, the cracked in the sludge deposits was discovered to have caused the influent sludge to "bypass" the wetland bed, which has negatively affected the overall system performance. Therefore, the model simulation in cracked/bypassed cases of sludge deposit should be studied and investigated as well.