NITRATE REMOVAL IN WOODCHIP DENITRIFICATION BIOREACTOR – AN APPROACH COMBINING MATHEMATICAL MODELLING AND PI CONTROL

X Mathematical model contains nonlocal condition representing the PI controller. X Mathematical approach is applicable to optimize the efficiency of nitrate removal. X Linear regression is obtained for nitrate removal rate w.r.t temperature and pH. X Linear regression is obtained for chemical reaction rates w.r.t temperature and pH. X While nitrate concentration varies over time flow rate adjustment takes place. Abstract. A mathematical model of nitrate removal in woodchip denitrification bioreactor based on field experiment measurements was developed in this study. The approach of solving inverse problem for nonlinear system of differential convection-reaction equations was applied to optimize the efficiency of nitrate removal depending on bioreactor’s length and flow rate. The approach was realized through the developed algorithm containing a nonlocal condition with an in-corporated PI controller. This allowed to adjust flow rate for varying inflow nitrate concentrations by using PI controller. The proposed model can serve as a useful tool for bioreactor design. The main outcome of the model is a mathematical relationship intended for bioreactor length selection when nitrate concentration at the inlet and the flow rate are known. Custom software was developed to solve the system of differential equations aiming to ensure the required nitrate removal efficiency.


Introduction
Release of nitrate ( -3 NO ) from agricultural sources is a significant surface water quality problem occurring in many areas around the world. This problem is particularly severe in humid climates where subsurface (tile) drainage systems have been installed.
Research based knowledge is a prerequisite for successful bioreactors design since the highest bioreactors efficiency is achieved with the maximum knowledge. Therefore, the design of woodchip denitrification bioreactors has been studied at the laboratory on both pilot and field scales (Addy et al., 2016;Christianson et al., 2011b;Greenan et al., 2009;Hoover et al., 2016;Lepine et al., 2016;Povilaitis et al., 2018). Subsequently, various design approaches have been proposed. For instance, (Christianson et al., 2011a) suggested a bioreactor design method which is based on 10-20% peak drainage flow and hydraulic retention time of 6-7 hours. This approach has received much attention in practical applications in United States. Another design concept correlated bioreactor surface area and treatment area allowing easier estimation of bioreactor volume (Verma et al., 2010). The design approach based on nitrogen mass removal concept was proposed in . There has been little discussion in the literature about the effect of bioreactor length-to-width ratio and cross-sectional shape on bioreactor performance. According to (Christianson et al., 2011a(Christianson et al., , 2012 the highest bioreactor efficiency could be achieved when the ratio is around 10. The nitrate removal between trapezoidal and rectangular cross-sections does not show any significant differences. Overall, to date there is no consensus regarding the optimal drainage bioreactor design method and optimal bioreactor parameters (Christianson et al., 2012). Various methods result in different bioreactor sizes and efficiencies (Christianson et al., 2012).
Although the number of investigations on denitrification process in woodchip bioreactors has significantly increased, a great interest to elaborate nitrate removal modelling tools for bioreactors design still remains relevant.
Mathematical modelling can often be used for better assessment of chemical transport, optimization, estimation and design of pollutants removal operations . The feedback loop control of influent flowrate based on the concentration in the effluent has been applied in previous works (Torà et al., 2014). In the paper (Halaburka et al., 2017) woodchip bioreactors were applied to remove nitrate from drainage runoff. A multispecies reactive transport model with Michaelis-Menten kinetics was developed to explain the concentration profiles of dissolved oxygen, nitrate and dissolved organic carbon and four additional models were developed based on simplifying assumptions.
Over the past decades, researchers have developed a number of mathematical models to simulate the fate and transport of nitrates, oxygen and products of the reactions in the reactors. Majority of them where simplified to enzymatic conversion of nitrate in anaerobic media and kinetic rate (zero and first order) expressions.
The main purpose of this paper is to present a mathematical model for the processes within the woodchip denitrification bioreactor applicable not only to simulate chemical transport of nitrates and oxygen, but also to control and optimize the nitrate conversion efficiency. The primary task of the model containing the control mechanism is to maintain the variable water flow for a required (set-point) output -3 NO concentration. The model could serve as a tool for better bioreactor design.
A pilot-scale cube-shape plastic denitrification bioreactor was used for the real experiments. The experiments were carried out under in turn flowing-through (i.e. continuously flowing) and non-flowing water (i.e. batch) conditions. Due to two types of experiments two mathematical models of denitrification process (for the non-flowing and for the flowing-through conditions) were proposed and analysed.
Models belong to a class of problems, namely differential equations subject to nonlocal conditions (Štikonas, 2014). The nonlocal conditions of the differential problem define the relationship between the values of the solution at the boundary, inner points and parameters of the equations. The real experiments of the denitrification process were mathematically described as a systems of two convection-reaction equations, without regard to auxiliary conditions (turbulent flow, etc.). Systems of convection-reaction and convection-reaction-diffusion equations are usually used for the denitrification process modelling as well as for wastewater treatment and simulations of reactive settling of activated sludge (Bürger et al., 2016). A similar model with the added source of carbon was proposed in (Lee et al., 2017) to predict bacterial nitrate removal in groundwater. The suggested kinetic model combines Monod kinetics and a constant denitrification rate.
The experiment with non-flowing water conditions was analysed first. Nitrate and dissolved oxygen concentrations were found through experiments. The constructed mathematical model allowed to calculate rate constants for the analysed chemical reactions.
Afterwards the system of differential equations was supplemented with a convection term. This improvement of the existing model with non-flowing water condition along with the experimental data enabled to predict and calculate the flow rate of water for treatment.
Feedback loop control uses nonlocal boundary condition as a measuring element (process value) and convection (flowing-through rate) as a controller output. Process control has been analysed in (Ivanauskas et al., 2017) using a PID (proportional-integral-derivative) controller (Åström & Hägglund, 1995). A similar model was created using a PI (proportional-integral) controller. The controller in this model adjusts the water flowingthrough rate in order to satisfy the set-point outflow concentration. Mathematical model with PI control allows to design systems capable to manage variable inflow -3 NO concentrations.

Field experiment
A pilot-scale cube-shape plastic denitrification bioreactor (1.0 m 3 volume) was placed below the ground by the excavation of a 1.2 m trench constructed at the Drainage Laboratory of Vytautas Magnus University, Lithuania (Figure 1). The bioreactor's container was fed by water from two interconnected plastic water tanks (1.0 m 3 volume each). The container was filled with mixed woodchips made from local raw materials. Alder (Alnus glutinosa) and pine (Pinus sylvestris) tree scraps dominated in the woodchips with prevailing (at 65% of the cumulative distribution) particle diameter varying from 1.1 to 3.0 cm (bulk density of 260 kg/m 3 ). The bioreactor was filled with woodchips to a depth of 1.0 m, and a saturation level of 0.90 m was maintained. A polyethylene liner was also folded over the top of the bioreactor to seal it from the soil and a mound with a 20 cm thickness was formed and sown with grass. The woodchip porosity was determined using a standard porosity determination procedure described by (Christianson et al., 2009). The analysis revealed that woodchip porosity was 56%.
The experiment was carried out under in turn flowingthrough (average retention time = 3.10 hours calculated as the ratio between bioreactor's pore volume and flow rate) and non-flowing water conditions. The -3 NO removal efficiency (determined as the difference between the inlet and outlet -3 NO concentrations divided by the inlet concentration) tests in bioreactor started on July 20, 2017 and the results cover the period until June 10, 2018. The water from the tanks was supplied to the bioreactor by gravity. The flow rate was determined by the difference in hydraulic head (max 3.6 m) between the water levels in the tanks and in the bioreactor. The inflow and outflow rates were adjusted manually using valves. The bioreactor was fed nitrate (via the addition of 3 NaNO to the water tanks) at concentrations ranging from 28.0 to 132.0 mg·L −1 with an average value of 66.1 mg·L −1 . These concentrations are typical (83% of the cumulative frequency) of the range of -3 NO values observed in drainage water under field conditions. The outflow concentrations ranged from 16.0 to 98.0 mg·L −1 (average of 42.4 mg·L −1 ). Therefore, the -3 NO removal efficiency changed from 17.5% to 70.8% (average of 37.1%).
During the study period the water temperature at the inlet and outlet ranged from 13.9 °C to 19.4 °C. The pH values in the inflow ranged from 5.2 to 8.6, and those of the outflow from 5.0 to 8.3. The dissolved oxygen concentrations at the inlet and outlet ranged from 3.2 to 4.4 mg·L −1 and from 0.0 to 1.2 mg·L −1 respectively.
The measurements were performed at various irregular time intervals varying from 16.43 to 183.3 hours by applying the same sampling procedures. In total 41 experiments were carried out, the results of which were used in this article. Sampled at the outlet, nitrate and oxygen concentrations were measured at the beginning and at the end of the experiments. The presence of -3 NO was determined via the spectrometric method using a Photometer MD600/MaxDirect (accuracy ±0.5 mg·L −1 ) system with powder reagents. The dissolved oxygen content (accuracy ±1.5%) and the water temperatures (accuracy ±0.2 °C) were measured with a portable HI-9142 (Hanna® Instruments Ltd.) multimeter, and the pH values were registered by a HI-98136 meter (accuracy ±0.1 pH).

Physical model
Model is constructed on the base of chemical processes in the bioreactor described in typical kinetics and competition equations (Tinoco et al., 1995). Cellulose degrading microorganisms grow on woodchips and produce extracellular enzyme cellulase. This enzyme catalyzes the hydrolysis of cellulose, which produces soluble monosaccharides (usually glucose) and a variety of oligosaccharides of different lengths. Cellulase producing microorganisms grow in aerobic media. When oxygen is consumed, production of the enzyme stops. However previously produced cellulase still works. Anaerobic cellulases producing microorganisms also exists. They consist of about 5-10% of aerobic cellulases producing microorganisms (Leschine, 1995). They produce extracellular huge protein complexes that adsorb on cellulose surface. They catalyse hydrolysis of cellulose, and, probably, further destruction of soluble sugars to acetic acid, or lactic acid. Complexes are strongly inhibited by glucose. Thus, number and variety of soluble carbon source depends on oxygen concentration.
Under anaerobic conditions some microorganisms can switch from oxygen to nitrate. A different number and variety of volatile nitrogen oxide products will be produced (the third path). This process will continue until the previously produced cellulase inactivates and no more soluble sugars will be produced. When nitrate has been consumed, some microorganisms can switch to sulphates. However, sulphates can only be consumed at very high concentration of sulphates. It means that the process of carbon consumption is complicated, and the rate of carbon consumption cannot be expressed by one equation.
Due to the decrease of oxygen concentration the production of soluble sugars also will decrease. Thus, for the optimal nitrate removal process a sufficient oxygen concentration is necessary.
We can regulate the nitrate removal rate by the regulation of flow rate through the bioreactor. At a very low flow rate all oxygen will be consumed at the inlet of the bioreactor, and the efficiency of the reactor will be low. By increasing the flow rate we will involve more of the bioreactor's content into the process, thus increasing the rate of nitrate consumption. At a high flow rate in the bioreactor the aerobic conditions will dominate, and the rate of nitrate removal will be low. For the precise regulation of the nitrate removal the concentration of the nitrates and concentration of oxygen should be controlled at the outlet of the bioreactor. It will be used as an essential information for the effective regulation of the reactor using a PI controller.

Mathematical model
A mathematical model of the nitrate removal in woodchip denitrification bioreactor that guarantees the required water purity was composed. The model is based on differential equations which were used for the analysis of the water treatment experimental results and the water treatment processes. The inverse problem for the system of reaction equations was solved for the analysis of denitrification processes.
Since the process of water treatment involves complex chemical reactions, in this paper the variation of nitrate and oxygen concentrations during water treatment was analyzed. Regression analysis was used to analyze the experimental results. The dependence of the parameters characterizing water treatment on temperature and pH was determined. Water flow rates which allow to reduce the -3 NO concentration to the desired level were presented.
The main chemical assumptions are: the oxygen concentration decreases according to a first order reaction and the -3 NO concentration decreases nonlinearly during the water treatment process. Since the real experiments were conducted on non-flowing water conditions and woodchip porosity was 56%, an assumption was made that the medium can be considered to be homogeneous, and it is sufficient to analyze the changes over time.

Non-flowing water model
First the data from the non-flowing water experiments were analyzed. Every experiment was taken as a single observation in a given dataset. The variables and dimensions for each experiment are described in Table 1. The non-flowing water experiments were used to analyze and compute the rate constants for the model, as well for the statistical analysis of the experiments. The inverse problem was solved for the system of two differential equations supplementing the initial and boundary conditions: O ( ) C tdissolved oxygen concentration (mg · L −1 ); k 1 -nitrate and k 2 -oxygen removal reaction rates (h −1 ); a -oxygen domination proportion; t -time (h); T -reaction duration (h).
The additional conditions for the inverse differential problem are formulated using concentrations at the beginning (t = 0) and at the end (t = T) of the experiment: The system should be solved for the rate constants k 1 , k 2 . Nitrates are consumed relatively to the oxygen concentration decline: To start with the solution for each experiment the value of k 2 was determined from the inverse problem (1b) with the initial conditions (2a) and the boundary conditions (2b): (4) k 1 was obtained from the model fitting for differential equations (1a, 1b) with the initial conditions (2a) and the boundary conditions (2b). The initial guess of * 1 k was used for the optimization algorithm: (5)

Flowing-through water model
An additional parameter (water flow) was added to the model described in the previous subchapter. Here we study the case when the water passes at a constant rate through the bioreactor. The modelled bioreactor structure was assumed to be homogenous with evenly distributed woodchips, laminar flow and unchanged inner distribution of microorganisms due to the flow. The idea was to look for the water flow rates which would guarantee the required -3 NO removal. The mathematical model based on the system of two convectionreaction equations (6) was applied.
The flowing-through denitrification bioreactor mathematical model is defined as a system of two first order convection-reaction nonlinear differential equations: ; C t x -dissolved oxygen concentration (mg·L −1 ); k 1 -nitrate and k 2 -oxygen removal reaction rates (h −1 ); α -oxygen domination proportion; t -time (h); V -water flow rate (m·h −1 ); a -length of denitrification bioreactor (m); T -reaction duration (h). The initial conditions specify the concentration distribution inside the bioreactor at the initial time moment (t = 0): where c 1 , c 2 -initial nitrate and dissolved oxygen concentrations (mg·L −1 ). The boundary conditions define the nitrate and dissolved oxygen concen trations at the points of inlet during the reaction ( 0 t T < ≤ ): NO concentration Q at the point of outlet was determined: To find the solution of the system of differential equations with the initial and boundary conditions the finite difference (explicit forward difference at time) methods (Baronas et al., 2009;Samarskii, 2001) and the Newton-Raphson (or secant) method for optimization (Kochenderfer & Wheeler, 2019) were used. Custom software was developed to solve the equations for a given models using Python programming language with a NumPy and SciPy libraries.

Flowing-through water model with monitoring
This mathematical model allows to compute the optimal water flow rate for a given nitrate inflow concentration.
The refined model allows to control the nitrate concentration in the effluent by manipulating the influent flowrate.
The mathematical model is very similar to the previous one. The differential equations are the same as (6), but the nonlocal condition representing the control mechanism is defined as a PI controller and the water flow rate V is a function of time t: The controller adjusts the water flowing-through rate ( ) V t in this model to satisfy the set-point outflow concentration , p i K K are nonnegative coefficients for a proportional and integral terms in a PI scheme. The error function ( ) e t is the difference between the required -3 NO concentration and the calculated one The PI controller continuously evaluates the error value ( ) e t and attempts to minimize it over time by adjusting the control variable ( ) V t to a new value defined by (10). If ( ) 0 e t < , the ( ) 0 V t = was accepted. The principal scheme of the computational experiment is shown in Figure 2.

Results and discussion
The results obtained from the mathematical model and statistical analysis are presented in this paragraph. The nitrate removal average rates are discussed as well. Experiments were performed under different water temperatures, acidicy (pH) and durations.

Nitrate removal rates
The nitrate and oxygen average removal rates are very important characteristics of the denitrification process. They are defined as is 0.0722. This indicates there was no strong interdependency between the nitrate removal rate and fluctuations in temperature and pH. Previous investigation (Povilaitis et al., 2018) shows, that the nature of the woodchips does not impact on the overall rate of nitrate removal. Therefore, we propose control mechanism to ensure the required nitrate concentration at the outflow.

Flow rate selection for the required nitrate outflow concentration
It is important to maximize the efficiency of the bioreactor. This can be achieved by the manipulation of water flow rate. This subchapter presents a mathematical model for the flow rate control. The computed flow rate ensures the required purity of the outflow water providing supreme flow rate.
In the experiments the nitrate inflow concentration 3 NO ( , ) C t a was in the 28 -132 mg·L −1 range. Having the upper and lower inflow concentration bounds the corresponding flowing-through rate for the bioreactor with NO concentration was computed and set to 2.22 mg·L −1 .
The choice of a set-point parameter was based on the fact of sulfate presence which can compete and act as an alternative electron acceptor when more reducing conditions develop. Consequently, a sulfate reduction normally occurs when -3 NO -concentrations have been substantially depleted (below 2.22 mg/L) . Figure 3 shows the computed results of the mathematical model for five randomly picked experiments. For a given range of the -3 NO concentration in the drainage water there can be found a flowing-through rate which allows to reach a safe outflow -3 NO concentration (2.22 mg·L -1 ). Each experiment had its own set of parameters (Table 3).
From the results it follows that under a higher inflow -3 NO concentration the flowing-through rate should be reduced. As an example, the experiment number 12 was analysed. Under a 40 mg·L −1 inflow concentration the flow rate cannot exceed 5.42×10 −3 m 3 ·h −1 , but when the -3 NO concentration exceeds 120 mg·L −1 the flowing-through rate should be limited to 4.29×10 −3 m 3 ·h −1 .

Bioreactor length selection
For the purposes of water treatment, it is important to ensure that the water purity requirements are met while the outflow volume is maximized.
In case the desired rate of outflow is fixed, this objective can be achieved by adjusting the length of the bioreactor.
The mathematical model can be used to calculate the length of the bioreactor which ensures the required outflow rate and 2.22 mg·L −1 nitrate con-centration. For this purpose the inverse problem for a system of differential Eqs (6) was solved at a given flow rate V.
The results of the numerical experiments showed linear dependence between the length of the bioreactor needed to maintain the "safe" -3 NO concentration and the flowing-through rate.
The linear formula for a given dataset is (12) measured in m: For example, if the reaction rates are k 1 = 0.00970 and k 2 = 0.05809, then the required flow rate of 1.38·10 −2 m/s is sufficient for the 1 m length bioreactor. Whereas for the 10 m length bioreactor, the maximum flow rate should not exceed 1.38·10 −1 m/s. The width and length of a bioreactor are equal to 1 meter each.

Flow rate control for variable inflow concentration
The water treatment optimization is more difficult when the inflow water has a varying nitrate concentration. This subchapter presents the results of calculations for the case where the input water nitrate concentration varies as a sinusoidal function. The results are presented in Figure 4 and Figure 5. In this particular case the inflow -3 NO concentration varies between 55-82 mg·L −1 . The outflow concentration is set to be maintained by the PI controler, the set-point is 2.22 mg·L −1 .
The Figure 5 shows the PI control result for a given parameter set. We can observe the outflow rate and -3 NO concentration to be approximately the same function but shifted in time. It is due to lag in the control signal and its propagation in -3 NO concentration change. In our example the set-point and the actual measured concentration have a spread within 20%. Each experiment is labeled with a number. The experimental data is presented in Table 3 Table 3

Conclusions
A series of field experiments in a denitrification bioreactor were performed aiming to develop a mathematical model for nitrate removal in a woodchip denitrification bioreactor. The nitrate removal process was modelled under non-flowing and flowing-through water conditions. The proposed mathematical approach is based on the system of nonlinear differential equations describing the kinetic and convection processes. The chemical reaction rates of denitrification process were calculated by solving the inverse problem for the differential system for each experiment event ( ). Moreover, the method has also been proposed to maximize the efficiency of a bioreactor efficiency by adjusting water flow rate.
The main outcome of the model is a mathematical relationship intended for bioreactor length selection when -3 NO concentration at the inlet and the flow rate are known. The calculated length ensures the required nitrogen removal efficiency in the outflow. In addition, a feedback loop control mathematical algorithm for flow rate selection was elaborated. Therefore, the distinctive feature of the developed model is the nonlocal condition representing the control mechanism which is defined as a PI controller. The model with nonlocal boundary conditions can be applied to other biotechnology applications.  Outflow C(NO3), mg ·L −1 Flowthrough rate C(NO3 ) outflow C(NO3 ) outflow set-point