Mathematical modelling of the internal circulation anaerobic reactor by Anaerobic Digestion Model No. 1, simultaneously combined with hydrodynamics

In this study, the hydrodynamic characteristic of a lab-scale internal circulation (IC) anaerobic reactor was investigated. We found that the gradual Increasing-Size Continuous Stirred-Tank Reactors (ISC) model is desirable to describe the hydraulic character of the reactor. As a generalized anaerobic digestion model, Anaerobic Digestion Model No.1 (ADM1) was combined simultaneously with the ISC model to simulate the effluent of the IC reactor. Both the stable running and overloading shock tests were carried out to validate the simulation. The mathematical simulation results agreed well with the experimental observation. This proposed model may be valuable to simulate the performance of the IC reactor effectively and to supply a useful tool to for designing and operating other anaerobic reactors.

gradients, and these inhomogeneities is bound to lead to uncertainties on the reaction result. Amini et al. 23 ever used two-and three-phase computational fluid dynamics (CFD) simulated the full-scale Membrane Bioreactor, by analysed fluid-flow pattern, shear stress and cross-flow velocity. Based on smoothed particle hydrodynamics model which was coupled with Activated Sludge Model No. 1, Meister et al. 24 developed a wastewater treatment model, which successfully computed spatial distribution of compounds in full-scale plant basins and the quantified the evolution of compounds. Hulle et al. 20 pointed out that a correct description of the mixing behaviour of an anaerobic digester is necessary to accurately model and predict the reactor performance. Consequently, the simplified mixing pattern might lead to mistaken calibration of the bio-kinetic model 25 . For this reason, it is meaningful to consider the hydraulic characteristics of the reactor when ADM1 is implemented.
This study aimed to effectively simulate the effluent performance of the anaerobic reactor. A lab-scale IC reactor was set up, the hydraulic characteristics of which was investigated by the lithium-ion tracing technique. Three kinds of tank-in-series models were used to depict the hydraulic character of the reactor, and we found that the Increasing-Size Continuous Stirred-Tank Reactors (ISC) model was the most suitable. Although ADM1 was used in many previous anaerobic digestion studies, we are the first to simultaneously consider the hydraulic dynamic of the reactor when implementing ADM1. By combining it with the ISC model, the ADM1 mathematical model simulated the effluent COD (COD eff ) both in steady and in over-loading states, the result of which fitted the experimental observation quite well. This mathematical model is expected to provide a beneficial method to definitely and effectively simulate the anaerobic reactor, as well as to provide several useful details on the reactor performance.  experimental setup and operation. The lab-scale IC reactor consisted of a synthetic wastewater tank, a peristaltic pump, an IC reactor and a wet gas flow meter, as shown in Fig. 1. The reactor was made of Plexiglas and had a cylindrical configuration. The internal diameter was 200 mm, and the height of the reaction zone was 1125 mm, with a total working volume of 26 L. The ratio of the under-expanded sludge bed compartment to the upper precipitation compartment was 4:1. The reactor was equipped with a water jacket, and all of the experiments were performed under a mesospheric temperature (37 ± 2 °C).

Methods
The IC reactor was inoculated with methanogenic granule sludge from a full-scale IC anaerobic reactor that treats pulp and papermaking wastewater (Guangzhou City, China). Then, 8 L seed sludge was inoculated and sparged with pure nitrogen gas for a half hour before use. The reactor was cultivated with synthetic wastewater for one month, with a 24 h hydraulic residence time (HRT). When the COD removal rate steadily reached 80%, the reactor was considered to be in stable running condition. experimental data collection. Residence time distribution (RTD) research based on the tracing technique is a common method of studying the mixing pattern of bio-chemical reactors 26 . In this study, eight tracing runs were carried out to investigate the RTD of the IC reactor at 4 different HRTs. The R1-R4 were conducted under normal operation with sludge, and the R5-R8 were conducted under the IC reactor without sludge (shown in Table 1).
It has been reported that with lithium as the tracer, lithium-ion (Li + ) does not adsorb onto the sludge and is not harmful to the biomass 14 . By using an injector, a solution containing 2.397 g Li 2 SO 4 •H 2 O (260 mg Li + , Shanghai Aladdin, China) was injected instantaneously to produce 10 mg/L mean concentration Li + in the reactor. The pulse injection (less than 2 second) 24 was conducted at t = 0, which was exploited to obtain the exit-age function E(t) of lithium-ion (represented as the substrate) 27 . Samples were collected from the effluent pipe every 0.5 hour for 4 h and 8 h HRTs, and every 1 hour for 16 h and 24 h HRTs. Samples were centrifuged and filtered before analysis 16 . The lithium-ion concentration was detected by an inductive coupled plasma emission spectrometer (ICP, PerkinElmer-Optima 8300, US) at a wavelength of 670.8 nm, according to the Standard Methods 28 .
Effluent water COD was determined using the potassium dichromate method 29 . The concentration of volatile fatty acids (VFA) was measured by a gas chromatograph external standard method. To inhibit the adsorption of organic acids by using a chromatographic column, all of the effluent samples were pre-treated by using a high-speed refrigerated centrifuge and a microporous membrane filtrate. The operating condition of the chromatographic detection was set as: the injector and hydrogen flame ion detector (FID) temperatures were 200 °C and 300 °C, respectively; nitrogen was used as the carrier gas under a 5 mL/min flow rate; and a DB-FFAP (30 m × 0.32 mm × 0.5 µm, Agilent, China) capillary column was used. The configuration of the column operating temperature was originally 80 °C holding for 2 minutes, then increased to 200 °C at the rate of 10 °C/min and held for 5 minutes before running detection. theoretical interpretation. Normalized RTD curves. The hydraulic characteristics of the IC reactor were determined based on RTD curves. The curves were obtained from tracer studies. The RTD exit-age distribution function E(t), the mean residence time t and the distribution variance of the RTD σ t 2 are calculated as follows 12 : www.nature.com/scientificreports www.nature.com/scientificreports/ where t is time, h; and C(t) is the tracer concentration at time t, mg/L. To compare the flow patterns at different HRTs, Eqs (1-3) were normalized by the mean residence time t. The normalized time θ, normalized RTD exit-age distribution function E(θ) and the dimensionless variance σ θ 2 are described as follows: Dispersion number & Péclet number. Theoretically, when the back-mixing in the reactor is weak, the axial dispersion plug flow (PF) model can be applied to simulate the hydrodynamic of the reactor. The dispersion number D/υL (dimensionless) is used to characterize the back-mixing intensity of the reactor system. The Péclet number (Pe, dimensionless) is the inverse of the dispersion number D/υL. Pe can be estimated as follows: An ideal totally mixed reactor will show Pe = 1, while the ideal PF will show Pe = ∞ at the other extreme 14 .
Tank-in-series (TIS) models. As a non-ideal flow, the hydrodynamics of the anaerobic reactor can be well described by the multi-CSTRs (continuous stirred tank reactor) tank-in-series model 25,29 . Figure 2 shows these models. The equal-sized CSTRs model (ESC) is commonly used, which is described by the Erlang distribution: The number N value (N ESC ) of ESC in Eq. (8) can be counted as follows: If N ESC is a non-integer, N ESC needs to be rounded to the nearest integer number. Furthermore, an extended equal-sized CSTRs (EESC) model can be used by the subset of the gamma distribution 27,30 , while the N value (N EESC ) is a non-integer in the equation, as follows: The previous research studies have illustrated that the distribution of the substrates and intermediates along the axis of the reactor were far different, and the increasing-size CSTRs model (ISC) would perform better at simulating the performance of the reactor 30,31 . The ISC model can be expressed by an equation set as below 19 : where r i is the volume fraction coefficient and N is the number of tanks.
ADM1 implement. Based on the hydraulic calculation result, AMD1 was used to simulate the effluent of the lab-scale IC reactor by AQUASIM2.0 (published by EAWAG) 32 . The biological degradation processes are described by substrate uptake Monod-type kinetics equations, while other extracellular processes (e.g., disintegration and hydrolysis) and biomass decay are described by first-order kinetics equations 33 , Despite all parameters in the equations theoretically affecting the outcome of the model, the sensitivity degree of the parameters was far different. In particular, most of the parameters have low sensitivity. A sensitivity analysis was used to identify the dominant parameters, which was conducted by employing the absolute-relative sensitivity function supported by AQUASIM2.0. The function is calculated as follows 32 : where y is an arbitrary variable that is calculated, and p is a model parameter. The absolute-relative sensitivity function measures the absolute change in y for a 100% change in p. Their units do not depend on the unit of the parameter, which makes quantitative comparisons of the effect of different parameters p on the same variable y possible.
The model parameters that are represented by constant variables can be estimated by minimizing the sum of the squares of the weighted deviations between the measurements and calculating the results of the model as follows: where y meas,i is the i-th measurement; σ meas,i is its standard deviation; and y i is the calculated value of the model. Using the algorithm provided the by AQUASIM2.0, both the sensitivity analysis and the model parameters were calculated.

Results and Discussion
HRt distribution in the trace experiment. Using the data collected from the R1-R8 tracing experiment, we obtained the HRT distributive information C(t)-t curves (Fig. 3), which were converted to E(θ)-θ curves (Fig. 4) 12,21 . By comparing the C(t)-t curves of the tracing experiment at different runs, we found that the C(t)-t Li + peak value intensity of normal operation was lower than that of the no sludge operation, which confirmed the findings of a previous tracing study 12 . Nevertheless, when the C(t)-t curves were normalized, and the E(θ)-θ peak value intensity of the with-sludge operation was obviously higher than that of counterparts without sludge, which reflected that the RTD in normal operation was intensive and the flow pattern tended to be PF 14 .
From the C(t)-t curves or E(θ)-θ curves, the peak value occurrence time (PVOT) of the normal operation was later than that without the sludge operation on a short HRT (comparing R4 and R8). Along with the extended HRT, the PVOT of the two operations were close (comparing R3 and R7, R2 and R6, and R1 and R5). These results illustrated that: the reactor structure fundamentally determined the RTD of the substrates; the sludge can detain the substrates, which showed apparently at the short HRT; nevertheless, the sludge detaining effect would degrade along with the extension of HRT, which might be attributed to the perturbation effect of generating bubbles in the reactor.
Back-mixing degree and degradation efficiency. Back-mixing refers to the mixing degree of substrates in the reactor, which exerts an important effect on pollutant removal 12 . The back-mixing degree can be described by N or D/υL 34 . Table 2 shows the RTD analyses results (D/υL, Pe and TIS model tank number N value) under all experiment conditions. From Eq. (9), the calculated tank number N value is able to identify the mixing pattern. If N tends to be 1, the reactor approximates to the CSTR (totally mixed); on the contrary, when N tends towards ∞, the reactor approximates to the PF reactor (without mixing) 34 . According to Tomlinson et al. 35 , D/ υL = 0.02 (Pe = 50) is defined as the intermediate degree of dispersion, and D/υL ≥ 0.2 (Pe ≤ 5) is defined as large dispersion.
The Run5-8 were carried out without sludge, and the N values were 2.28, 2.32, 2.15, and 1.93 under the HRT, from 24 to 4, respectively, which revealed that the mixing pattern tends to be completely mixed, and the back-mixing degree was increased with the decreasing HRT. The dispersion number D/υL values (0.31, 0.31, 0.35, and 0.42, respectively) were greater than 0.02. Therefore, the dispersion of the reactor was large and the back-mixing was strong. The Run1-4 were executed under normal conditions (with sludge), and the N values were 2.25, 2.41, 2.66 and 2.99, respectively, under HRT from 24 h to 4 h, which implied that the mixing pattern also tends to be strongly mixed, but the back-mixing degree decreased with the decreasing HRT. The dispersion numbers D/υL (0.32, 0.29, 0.25 and 0.21, respectively) were still greater than 0.02 but were apparently less than the value of the counterparts without sludge. The result illustrated that the configuration of the reactor mainly www.nature.com/scientificreports www.nature.com/scientificreports/ determines the mixing pattern. Within the biomass, the granule sludge and the generation of bubbles produce a coaction to stimulate the mixing pattern to be PF. This effect was obvious under a short HRT.
Theoretically, the pollutant removal rate in the absolute PF pattern is higher than the rate in the absolute CSTR pattern. The calculated tank number N values under the normal experiment were from 2.25 to 2.99 (HRT from 24 to 4) and were no more than 3, which indicated that this IC reactor has a moderate back-mixing pattern between the PF and the completely mixed pattern 34 . With the decreasing HRT, the mixing pattern tends to be a PF pattern. Generally, when the reactor is dominated under Monod or first-order biological kinetics, a PF pattern will be more effective 36 . Consequently, if the anaerobic bioreactor has an absolute PF pattern, the production of VFAs will accumulate at the bottom layer and the pH will descend there; at the same time, the methane production occurs at only the top layer and the pH value will ascend 12 . It should be noted that methanogenic bacteria may be inhibited by a low pH, when the pH of the reaction is less than 6.5 37,38 . This adverse result can break down the methane production reaction.
In contrast, moderate back-mixing can produce a continuous stirring between the liquid-solid phases, which not only enhances the transfer action on pollutants and granule sludge but also increases the balance of substrates, the pH, and nutrients in the reactor. Therefore, the moderate back-mixing effect may be beneficial for the highly efficient degradation capacity of IC reactor.   Figure 5 shows the simulated results of the three models. It is easy to find that the accuracy of the ISC model is higher than that of the other two models. This finding agreed with the research results of Dai et al. 11 and Ren et al. 30 , who even pointed out that it was convenient to modify the ISC model to accurately simulate the hydrodynamic of the anaerobic bioreactor. In this study, the estimated N ISC values were acquired according to Eq. (9) calculation results, and the tank volume ratio was further estimated by minimizing the squared difference χ 2 . Therefore, exploiting the ISC model would easily improve the simulation results.  Table 3. Parameter Estimation Results for TIS Models. www.nature.com/scientificreports www.nature.com/scientificreports/ Sensitivity analysis of the ADM1 parameter. ADM1 involves 19 biochemical processes and 26 component variables 15 . All of the biochemical process rates are described by Monod uptake equations or by first-order kinetic equations 39 . Despite the fact that all of the parameters in the equations theoretically affect the output of the model, the sensitivity degree of the parameters may be far different from one to another, and many parameters have little impact. It is necessary to analyse and distinguish the important parameters which have the significant impact on the model calculation. Chen et al. 34 ever pointed out that sensitivity analysis would help to identify the important parameters and reduces the complexity of parameter tuning, who only chose k m parameters (specific Monod maximum uptake rate) of propionate and acetate for parameter tuning. Barrera et al. 33 employed local relative sensitivity analysis method to calculate sensitivity functions for the dynamic simulations. As the most sensitive parameters, k m parameters of propionate, acetate, hydrogen, and Yield of hydrogen were used for calibration in the anaerobic digestion with sulphate reduction of cane-molasses vinasse. The sensitivity analysis was a useful method to identify the dominant parameters, to reduce the model complexity and to determine the main processes 40,41 .
In the research, lab synthetic wastewater was used, and the glucose was the only substrate in the influent. Glucose, as a simple monosaccharide, undergoes easy uptake by sugar degraders (X su ) and will be decomposed directly into VFAs (butyrate, propionate, and acetate, etc.) and H 2 15,40 . The VFAs are the main ingredients of the effluent COD (COD eff ). Thus, to accurately predict the COD eff , the kinetic parameters of the uptake process rates of butyrate, propionate and acetate were used in the calculation.
Using ADM1 incorporated with the ISC hydraulic model, a sensitivity analysis of the parameters was executed by AQUASIM. The absolute-relative sensitivity analysis of the k m parameter (specific Monod maximum uptake rate) and the K s parameter (Monod half saturation constant) relating to these ingredients are shown in Fig. 6. The sensitivity of k m in the process rate equation is higher than that of K s in the same process rate equation. The k m and K s values of the uptake process in the acetate and propionate uptake processes showed higher sensitivity than their counterparts in the butyrate and sugar uptake processes, especially in the k m value of the acetate uptake process. Therefore, these four parameters were used for parameter adjustments, which would be considered in the next step of the mathematical simulation. The sensitivity curves of both k m and K s parameters for acetate and propionate are shown in Fig. 7.

Mathematical model development and validation.
Referring to the parameter estimation results for the ISC hydraulic model (Table 3), we chose the medium ratio value 1:2:5 for the developed ISC model (a 2 litre tank, 4 litre tank and 10 litre tank were joined together) to simulate the hydrodynamic of IC reactor. Then, ADM1 incorporated with this ISC hydraulic model were set as the simulation mathematical model.
Anaerobic reactor system was sensitive to influent fluctuation, such as high organic loading shock, hydraulic shock, and even toxic wastes shock etc. The disturbances might adversely affect the quality of effluent what should be monitored and predicted earlier. Consequently, to validate the mathematical model in wastewater COD eff prediction, both the stable running test and overloading shock test were carried out in our experiments. On table running test, the influence was maintain the basic COD concentration (3 g/L), and the HRT was maintain 24 h; on hydraulic shock test, the influence was six times flow rate lasting 12 hours, with basic COD concentration; and on the organic loading shock, the influence was six times COD concentration (18 g/L) lasting 8 hours, with 24 h HRT.
According to the outcome of the sensitivity analysis result, we adjusted the k m_ac , k m_pro , K s_ac and K s_pro values in this mathematical model simulation. Other parameters were adopted from the recommended values of the IWA task group 15 . These adjusted values are listed in Table 4.
The simulation result of the IC reactor start-up was on stable running, when the influent COD concentration was maintained at 3 g/L and HRT at 24 h (shown in Fig. 8). After adjusting the maximum uptake rate and the half saturation constant, the simulation COD effluent agreed with the experimental effluent COD. Comparatively, the www.nature.com/scientificreports www.nature.com/scientificreports/ simulation with the original specific Monod maximum uptake rate and Monod half saturation constant, did not fit the experimental result well, especially at the initial stage before the reactor system was running stably ( < 20 d).
Short-time overloading, such as a huge pollutant influent or water flow shock, are common in industrial wastewater treatments 42,43 . Based on regular loading simulation, we further simulated the over-organic and over-hydraulic shock tests. Organic loading shock lasting 8 hours was carried out six times, the result of which is shown in Fig. 9. The simulation COD eff matched quite well with the experimental COD eff . Over-hydraulic shock lasting 12 hours was also conducted six times. As shown in Fig. 10, the simulation result agreed with the hydraulic shock experimental COD eff . Both of these well modelled results validated that the mathematical model of ADM1 incorporated with the ISC hydraulic model can predict well the IC reactor effluent COD.    www.nature.com/scientificreports www.nature.com/scientificreports/

Conclusion
A lab-scale internal circulation (IC) anaerobic reactor (26 L) was set up, and the hydrodynamic characteristics of this IC reactor were studied. Three kinds of tank-in-series models were used to simulate the hydrodynamics of this reactor, and the gradually increasing-size CSTRs model (ISC model) was more suitable to fit the hydraulic character than the other two models. Based on the calculations of RTD dimensionless variance, 3 tanks in this estimate had been used by the TIS model. We chose the medium ratio value 1:2:5 of the ISC model to simulate the hydraulic dynamic of this reactor. Based on this, the ADM1 was combined with the ISC model to predict the effluent COD. Both stable running and overloading shock tests were executed, and the simulation results agreed with the experimental data. The proposed model in this research might be valuable to monitor reactor effluent effectively and to supply a useful tool to design and operate a full-scale anaerobic reactor.