Adapting Dynamic Mathematical Models to a Pilot Anaerobic Digestion Reactor

A dynamic model has been adapted to a pilot anaerobic reactor fed diary manure. Both steady-state data from online sensors and laboratory analysis and dynamic operational data from online sensors are used in the model adaptation. The model is based on material balances, and comprises four state variables, namely biodegradable volatile solids, volatile fatty acids, acid generating microbes (acidogens), and methane generating microbes (methanogens). The model can predict the methane gas flow produced in the reactor. The model may be used for optimal reactor design and operation, state-estimation and control. Also, a dynamic model for the reactor temperature based on energy balance of the liquid in the reactor is adapted. This model may be used for optimization and control when energy and economy are taken into account.


Introduction 1.Anaerobic digestion of animal wastes
Anaerobic digestion (AD) of animal wastes can produce biogas with methane to be used as an energy source, and a liquid effluent containing valuable nutrients.Moreover, AD reduces methane emission, odours and contaminants.AD bioreactors are effective as they allow for relatively high load rates (feed rates) and small reactor volumes.AD reactors may become unstable, i.e. a persistent decrease of gas production, because of inhibitory effects on methane-forming microorganisms due to large concentrations of volatile fatty acids (VFA) and ammonia and too low pH.Instability can also occur because of washout of microbes when the feed (load) rate is too large.
Various theoretical and practical aspects of AD processes are described in e.g.Tchobanoglous et al. (2003) and Deublein and Steinhauser (2010).A presentation of AD of animal wastes (dairy, beef, poultry, and swine) is provided e.g. by Husain (1998).

Possible applications of mathematical models
Foss Biolab, Haugen et al. (2012), is a pilot biological plant at Foss dairy farm in Skien, Norway, for nutrient and energy recovery from animal waste.The aims of this paper are to adapt a dynamic mathematical model of the anaerobic digestion (AD) processes of the reactor in the plant able to predict the methane gas flow produced in the reactor, and to adapt a dynamic model able to predict the reactor temperature.
Possible applications of a mathematical model of the AD processes in the reactor are as follows.
• Analysis of the dynamic and steady-state behaviour of the AD processes primarily based on simulations.Using simulations can provide process insight which would otherwise be practically difficult to obtain.
• Optimal design and operation of a full-scale reactor, i.e. designing optimal reactor size and calculating optimal feed rate according to proper optimization criteria, Edgar et al. (2001).
• Design of state estimators, also denoted softsensors, which are algorithms calculating the state of the process continuously, Simon (2006).The most common state estimation algorithm is the Kalman Filter which exists in different versions.A state estimator is an essential part of model-based predictive controllers, see below.It may also be used as a soft-sensor on its own as a substitute for real measurements.
• Model-based tuning of industry-standard PID controllers (Proportional + Integral + Derivative), Seborg et al. (2004), for biogas flow control to keep the produced biogas flow at or close to a given setpoint.
Possible applications of a mathematical model of the reactor temperature are: • Optimal reactor design and operation taking energy and economy into account.A combination of the AD process model and the reactor temperature model will be necessary to solve this optimization problem.
• Tuning a temperature controller for the reactor.

Outline of this paper
Section 2 gives a description of the pilot AD reactor at Foss Biolab.Section 3 describes adaptation of a dynamic mathematical model of the AD process to the reactor.Section 4 describes mathematical modelling of the reactor temperature.A discussion is given in Section 5, and conclusions are given in Section 6. Nomenclature including abbreviations is given in Appendix A. Laboratory analysis methods for relevant components are described in Appendix B. For easy reference, a summary of the modified Hill's model adapted to the pilot reactor is given in Appendix C. MATLAB (by The MathWorks, Inc.) has been used as computational tool for this paper.

Overview
Figure 1 shows a Piping & Instrumentation Diagram (P&ID) of the biological process plant at Foss dairy farm.Input to the plant is dairy manure diluted with 25% water and filtered, and outputs are fertilizer and biogas consisting of 70-75% methane.The plant is monitored and controlled with a PC running LabVIEW.The main parts of the plant are as follows (numbers refer to Figure 1).1.A reservoir for raw dairy manure with approximately 25% added water.
2. A sieve to separate the manure into two fractions of similar total solid mass: > 70 % of the volume is wetter fraction, and < 30 % is dryer fraction.The dryer fraction is used for vermicomposting.
3. A high rate anaerobic digestion reactor fed filtered cow manure as substrate for production of energyrich biogas that contains mainly methane.The effective reactor volume is approximately 250 L.
4. A nitrification reactor of approximately 200 L fed AD reactor effluent to produce liquid fertilizer and pellets fertilizer from formed foam.
This paper concerns the AD reactor in the process line.
The present reactor has been operational since April 2012, while a previous similar reactor was in operation from August 2011 until April 2012.

Instrumentation
Below is a list of the instrumentation used with the AD reactor depicted in Figure 1.The encircled numbers refer to the diagram in Figure 1.
PC is a laptop PC running the computer program for monitoring and control implemented in Lab-VIEW.The PC is connected to sensors and actuators via USB-based I/O-devices.
P2: A voltage controlled peristaltic pump which is operated using pulse-width modulation (PWM) with a fixed cycle time of 10 min.
There are several benefits of using PWM control compared with analog control: The calibration of the pump is needed only at one flow, namely the maximum flow rate (i.e. the flow with PWM signal in the On state), and it is therefore easier to obtain any flow (in average) within the minimum and the maximum flow ranges.PWM may also reduce blocking in the feed pipeline.
FT-1: AD reactor effluent flow sensor (home-made) based on measuring the frequency of effluent charging and discharging of a cup of fixed volume.
The sensor output is normal litres of gas at NTP (Normal Temperature and Pressure), i.e. temperature 0 o C and pressure 1.013 bar. of the biological process line of the pilot plant at Foss Biolab, Skien, Norway.This paper concerns the AD reactor.
The raw measurement signal from this sensor is quite noisy with an observed standard deviation of approximately 14 L/d (litres per day) which is approximately 2% of the upper range limit (URL) which is 720 L/d.To smooth the noise, a lowpass filter with time-constant of 0.2 d is implemented in the LabVIEW program.
To smooth the measurement signal from this sensor, a first-order lowpass filter with time-constant of 1 h (hour) is implemented in the LabVIEW program.
The measurement signal is filtered with a lowpass filter with time-constant of 1 h.
Filter: Lowpass filter with time-constant of 10 min.
Filter: Lowpass filter with time-constant of 10 min.
TC-1: Reactor temperature controller which is operated as an On-Off controller (i.e.thermostat controller).The controller is an industrial standalone temperature controller.
H-1: Electrical heater for the AD reactor which is controlled using the built-in PWM option in TC-1.
The heater comprises an electrical resistor wound around the reactor inside the thermal insulation jacket.

Available data
Data used for model adaptation are offline-data from laboratory analysis and online-data from sensors.Samples for laboratory analysis have been taken regularly from the reactor since August 2011, following sampling guidelines given by Esbensen and Paasch-Mortensen (2010).A number of different variables characterizing the reactor influent and effluent are analyzed.Among these, concentration of volatile solids (VS) and concentration of total volatile fatty acids (VFA) are used for model adaptation in the present study.
Online-data include feed flow (load rate), reactor temperature, ambient (air) and feed temperature (assumed to be the same and therefore measured with one sensor), biogas flow, and methane gas concentration.The latter two provide methane gas flow.

Adaptation of a mathematical
model to the AD reactor

Model selection criteria
For our purposes, cf.Section 1.2, a model is searched for according to the following criteria: 1.The model must be able to predict the produced biogas.However, it is sufficient that methane gas flow is predicted since the useful energy content of the gas is related to the methane content only.
2. The model should be relatively simple since it is to be used in a real-time implementation of a state estimator and a model-based controller.Relatively simple models are preferable since they may be easier to adapt and maintain.
3. The model should be able to represent the temperature dependency of the dynamics of the AD process.This is because the reactor may be operated at different temperatures, although the mesophile condition which is about 35 o C is assumed to give the optimal temperature condition for the methane-producing micro-organisms, or methanogens, Tchobanoglous et al. (2003).Optimality in terms of cost may imply a temperature being different from 35 o C.

Model candidates
A number of candidates of model for AD of dairy manure were considered in the light of the above criteria.
Below is a summarized characterization of these models.
• Andrews and Graef (1971): -Model characteristics: The model is general, and does not assume any particular type of organic substrate.Substrate: Acetic acid.-Gas predicted by the model: CH 4 gas flow; CO 2 gas flow; Biogas flow as sum of these gas flows.
• Hill and Barth (1977): -Model characteristics: The model is general for animal waste, and is validated using experimental data from reactors fed poultry waste and swine waste.Hydrolysis step is included.Model parameters are expressed as Arrhenius-based functions of temperature.
- -Gas predicted by the model: CH 4 gas flow; CO 2 gas flow; NH 3 gas flow; Biogas flow as sum of these gas flows.
• Hill (1983): -Model characteristics: The model applies to waste from poultry, beef, dairy, or swine with two parameters -the biodegradability constant and the acidity constant -being unique for each of these wastes.Substrate: Volatile solids (VS).A certain fraction of the fed VS is assumed biodegradable.Hydrolysis step is not included.Model is validated using experimental data with all four wastes mentioned above.Applicable temperature range is 20 o C -60 o C based on the temperaturedependency of the maximum reaction rates according to Hashimoto et al. (1981).
-State variables (4): Concentration of biodegradable volatile solids (BVS); Concentration of volatile fatty acids (VFA) as acetate; Concentration of acidogens; Concentration of methanogens; -Gas predicted by the model: CH 4 gas flow.
• Husain (1998): -Model characteristics: Husain (1998) presented Hill's model (1983) with more details regarding chemical reactions.Husain also changed some of the model parameter values, and expressed the death rates of the acidogens and metanogens as VFA-based Monod functions instead of relating the death rates to the maximum reaction rates as Hill did.
-State variables (4) are as in Hill's model.
-Gas predicted by the model as in Hill's model.-Gas predicted by the model: Biogas flow (the individual gas components are however not stated explicitly in Zaher et al. (2009).

Selection of the ultimate model
In the light of the criteria for model selection presented in Section 3.1.1Hill's model, Hill (1983), is selected as the ultimate model as it satisfies all the criteria and because it is simpler than comparable models, see below.Support for this selection is found in the evaluations of various AD models made by Stromberg (2010) and Husain (1998) where they conclude favourably about this model.Brief evaluations of the other models presented above in the light of our selection criteria are given in the following.
The model by Andrews and Graef (1971) is not selected because it is validated only at the fixed temperature 38 o C. Also, we see it as drawback that it contains only one type of microorganisms while two or more types are very common in more modern models.
The model by Hill and Barth (1977) is attractive for our purposes, but is not selected here because it is more complicated than Hill's model, Hill (1983).However, this model may be selected in future projects.
The model presented by Husain (1998) is basically the selected Hill's model.
The ADM1 model, Batstone et al. (2002), is not selected because it is very complex, and because numerical challenges may be expected in an online (real-time) implementation of state estimators and model-based controllers.Lyseng et al. (2012) adapted ADM1 quite successfully to our pilot reactor in the AQUASIM simulation tool, Reichert (1998), however with poor predictions of produced biogas as the reactor temperature was changed (simulated and measured experimental gas production were clearly different).
The model by Zaher et al. (2009) is not selected here because it is relatively complex, and the model parameters are not presented with any temperature dependency.Also, it would be necessary to modify the model originally made for batch AD processes to make it applicable to our bioreactor which has a continuous load rate.

Hill's AD model
The reasons for selecting Hill's model, Hill (1983), as the ultimate model to be adapted to the pilot bioreactor at Foss Biolab are given in Section 3.1.Hill's original model is presented in Section 3.2.1.A modified Hill's model which is the model adapted to the AD reactor at Foss Biolab is described in Section 3.2.2.Hill's model comprises eqs.( 1)-( 10) below, but a few changes are made in the symbols to make the model more readable and to have symbols which are in more compliance with symbols in other AD models.The differential equations stem from mass balances of the pertinent components.Homogeneous conditions are assumed.

The original Hill's model
The alphabetic identificators written in parentheses in the following refer to the flow diagram in Figure 2.
Defining that portion of the raw waste which can serve as substrate (A): Defining that portion of the biodegradable material which is initially in the acid form (B): Mass balance of biodegradable volatile solids (C): Mass balance of VFA (D): Mass balance of acidogens (E): Mass balance of methanogens (F): Methane gas flow rate (gas production) (G): The reaction rates are as follows: The maximum reaction rates µ m , µ mc are functions of the reactor temperature as follows, Hashimoto et al. (1981): The death rates are set to one tenth of the maximum reaction rates: Although Hill's model is selected as the ultimate model, we are motivated to implement a few changes to the model as explained in Section 3.2.2.

Differences from the original Hill's model
In our study the following changes are made to the original Hill's model presented in Section 3.2.1.The resulting model is referred to as "the modified Hill's model".• The Haldane functions in the reaction rates µ and µ c in Hill's original model, eqs.( 8) and ( 9), are replaced with the simpler Monod functions: This makes the calculations with the model in the context of parameter estimation easier.Using Monod functions is consistent with the comparable model by Simeonov et al. (1996).
• In the original Hill's model the retention time of the biomass (here: acidogens and methanogens) is equal to the hydraulic retention time (HRT): The retention time of the biomass is larger than the hydraulic retention time in up-flow sludge bed reactors such as applied here, where biomass is retained by gravity, Tchobanoglous et al. (2003).
The retention time ratio b is here introduced.The retention time of the biomass, which is denoted the solids retention time (SRT), is assumed to be b times the hydraulic retention time: where the term V /(F f eed /b) expresses that the biomass flow out of the reactor is smaller than the flow of organic matter.
In the original Hill's model it is implicitly assumed that b = 1.Eq. ( 18) makes the model coherent with the standard ADM1 model, Batstone et al. (2002), in this respect though the SRT is represented differently (as an independent parameter) in ADM1.Eq. ( 18) is in accordance with the representation of SRT in e.g.Zaher et al. (2003) and Bernard et al. (2001).
Model equations in the modified Hill's model Defining that portion of the raw waste which can serve as substrate: Defining that portion of the biodegradable material which is initially in the acid form: Mass balance of biodegradable volatile solids: Mass balance of total VFA (see comment below): Mass balance of acidogens: Mass balance of methanogens: Methane gas flow rate (gas production): where the reaction rates, with Monod kinetics, are as follows: The maximum reaction rates µ m , µ mc are functions of the reactor temperature as in the original Hill's model, eq. ( 10), repeated here for easy reference: Above it is assumed that VFA is total VFA consisting mainly of propionate, butyrate, valerate, and acetate, Batstone et al. (2002).Acetate is the main VFA component, and it is used in aceticlastic methanogenesis which is the main methane-generating process.Methane is also generated in hydrogenotrophic methanogenesis.Hydrogen is generated from various components including the VFA components propionate, butyrate and valerate.To include effects of the hydrogenotrophic methanogenesis, S vf a in our model represents total VFA and not only acetate.
Figure 3 shows an overall block diagram displaying the variables and parameters of the modified Hill's model eqs.( 19)-( 28).
X acid

Input variables
Output variable In Section 3.3.3about parameter estimation the steady-state (or static) version of the dynamic model is needed.The steady-state model is obtained by setting the time-derivatives in the above differential equations equal to zero.

Adaptation of modified Hill's model to the AD reactor
In the the following subsections, 3. Cf. comment in text Cf. comment in text • B 0 defines the ratio between BVS and VS in the feed: It is assumed that a proper value of B 0 can be found from the following specific long-term test.At time t 0 = July 12, 2011 a fresh sample of the subtrate was put into an incubator having constant temperature of 35 o C.
The biogas production was registered regularly until time t 1 = Sept. 25, 2011 which is the time where the biogas production became virtually zero.The VS concentration, here denoted S vsin , of the substrate was measured at times t 0 and t 1 , cf.Table 3.
Since the biogas production is zero at t 1 it is concluded that the biodegradable part of the substrate is completely degraded at t 1 .Thus, in the long-term test,

Parameters and variables to be estimated
The following parameters and variables in the modified Hill's model, eqs.( 19)-( 28), are to be estimated: In the model k 1 , k 2 , k 5 , K s , and b appear as parameters, while X acid and X meth are state variables.It is necessary to estimate X acid and X meth since their values are not known.

Identifiability
Before trying to estimate parameters and variables, their structural identifiability should be determined, Dochain and Vanrolleghem (2001).Structural identifiability concerns the possibility to give a unique value to the unknown parameters and variables, and this property can be assessed with a number of alternative methods, e.g. the method of Generating series which is based on calculating Lie derivatives, or by analyzing the uniqueness of the parameters of the Laplace based transfer function of a linearized version of the (nonlinear) model.However, for the estimation method used, see below, the structural identifiability is obvious since independent analytical expressions for the pertinent parameters and variables using steady-state data are obtained, except for one parameter, namely k 5 , which is found by optimization using dynamic data.Hence a rigorous structural identifiability assessment is not necessary and therefore not accomplished here.
Method of parameter estimation Yield-parameter k 5 is a crucial model parameter since it is directly related to the methane gas flow rate, cf.eq. ( 25).It is decided to estimate k 5 using nonlinear least squares (NLS) estimation based on optimization using iterated simulations of the modified Hill's model.The simulations are based on the explicit Euler method implemented with for-loops in MATLAB.The optimization problem is to minimize the difference between real (measured) and simulated F meth in the least squares sense over a specific time interval which is from t = 66 d to 95 d in Figure 4. (t = 0 corresponds to April 19, 2012.) The iterations (including simulations) are executed automatically by the optimization solver. 1rom t = 72.3d the methane gas flow is being controlled by a feedback controller manipulating the feed pump.The noise in gas flow measurement seen in Figure 4 is due to blockings and power outages.Because of the feedback control, this noise imposes noise in the feed flow, via the controller.
Before each of these simulations is started (automatically by the optimization solver) the six parameters/variables k 1 , k 2 , K s , b, X acid , and X meth are calculated from the steady-state version of the dynamic model ( 19)-( 28) using steady-state operational data at a specific operating point assuming intially that k 5 has a "guessed" value which is set equal to the value in eq. ( 7) of Hill's original model: Table 5 shows values of inputs and states in the pertinent steady-state operation point which is t = 66 d in  F meth = 227.9

Sensor reading
The formulas for calculating the above mentioned six parameters/variables, namely k 1 , k 2 , K s , b, X acid , X meth , from the steady-state version of the dynamic model ( 19)-( 28) are given below.
• b is calculated from the steady-state version (i.e. the time-derivative term is set to zero) of eqs.( 24) and ( 27) to give b = • K s is calculated from the steady-state version of eqs.( 23) and ( 26) to give • X meth is calculated from the steady-state version of eq. ( 25) to give • X acid is calculated with Here, r am is (in the steady-state calculations) set to r am = 3.4 (41) as it was observed in simulations for various feed rates with the original Hill's model with parameters from Husain (1998) that the ratio X acid /X meth varied only slightly around 3.4.However, this ratio is not necessarily 3.4 after the estimation procedure is finished.
• k 1 is calculated from the steady-state version of eqs.( 19) and ( 26) to give • k 2 is calculated from the steady-state version of eqs.( 22), ( 26) and ( 27) to give The simulations used in the optimization (to estimate k 5 ) are based on the dynamic model ( 19)-( 28).In the simulated model a lag of time-constant equal to is included.This time-constant represents the lowpass filter used in the real system at Foss Biolab to smooth the noisy biogas flow measurement, cf.Section 2.2. Figure 5 illustrates the estimation method used to estimate the parameters/variables (35).In this figure, p is the parameter to be estimated which is In Figure 5, the measured output (time series) is and inputs (time series) are The output of the optimization objective function f (and this output is to be minimized) is the sum of squares of prediction errors over the estimation time interval: E pred is the time-series (vector) of prediction errors: where N is the number of time-steps in the estimation time interval.In the present application, e pred (t k ) is where F methmeas is measured reactor methane flow, and F methsim is simulated F meth .In eq. ( 48), M represents the model comprising eqs.( 19)-( 28).
The optimal (best) estimate of p is the value of p which minimizes SSE e pred :

Results
In the time interval 80-95 d in Figure 4, the simulated F meth based on the estimated model is plotted together with real F meth .The simulation runs with initial state equal to the real state at t = 66 d.In this time interval the maximum difference between the simulated F meth and the real F meth is approximately 10 L/d while the maximum gas flow is approximately 235 L/d.It seems that the adapted modified Hill's model is able to predict the produced methane gas quite well.
Figure 6 shows in the upper plots real values (from laboratory analysis) and simulated values of the state variables S bvs and S vf a together with the respective real concentrations in the feed.Lower plots: Simulated values of the state variables X acid and X meth for which we do not have laboratory analysis data.
The lower plots show simulated values of the state variables X acid and X meth for which no laboratory analysis data are available.
Table 6 shows the values of the estimated variables and parameters together with standard deviations (σ), both absolute and relative, obtained with simulations as described later in the present section.

Uncertainty of estimates in terms of standard deviation
Uncertainty in the estimates can be expressed by the variability of the estimates.The variability may be calculated in a number of ways, e.g.:  Dochain and Vanrolleghem (2001).With a complex model analytical sensitivies may be calculated using computer tools for symbolic mathematics.
A possible alternative way of calculating the parameter estimation error covariances is with the calculation of so-called sigma points using the nonlinear model directly, without any linearization, as in the Unscented Kalman Filter, Simon (2006).
• Bootstrapping with parametric simulation, Davison and Hinkley (1997), which involves running a large number of simulations where the pertinent input data used in the estimation are varied randomly according to an assumed probability distribution which (should) resemble the actual distribution of the real data used in the estimation.
(Bootstrapping with parametric simulation resembles Monte Carlo simulations.)The parameter uncertainty can then be assessed from the observed variations of the estimates in the simulations typically in terms of standard deviation of the variation.
The method of bootstrapping is selected since this method is relatively straightforward and applicable.
Table 7 lists three quantities used as input data in the parameter estimation.These quantities are varied randomly and independently in the bootstrapping simulations.These quantities represent the steady-state operating point used in the parameter estimation, as explained above in the present section.Table 7 shows the pertinent standard deviations (σ) calculated from a number of laboratory analysis data sets taken over several weeks in the summer of year 2012.This time interval includes the steady-state operating point, namely day 66, see Figure 4).In the simulations it is assumed that their probability distributions are Gaussian.The resulting standard deviations of the estimates found from the bootstrapping simulations are shown in Table 6.Note that these standard deviatons express only the variations in the estimates due to variations in the three quantities given in Table 7.There are several other factors which contribute to the (total) uncertainty of the estimates, as model structure errors and systematic measurement errors.

Sensitivity of estimates with respect to assumptions
It is informative to assess the sensitivity of the estimates with respect to assumptions for the estimation method.The relative sensitivity of parameter p to parameter a is defined as where p 0 and a 0 are nominal values.p 0 is a parameter shown in Table 6. a 0 is a parameter shown in Table 2 or Table 5. p 1 and a 1 are respective values after a change is made in a.
Table 8 shows relative sensitivities as found by perturbing parameters a i with a 10% positive additive change and observing the corresponding change in the estimated parameter p.For example, S b,B0 = −0.62 in Table 8 means that a 10% additive increase in the assumed value of B 0 causes an additive change in the estimated b of −0.62 • 10% = −6.2%.
Note that X acid and X meth in Table 8 are actually estimated values for the steady-state operating point, and serve as initial values in the simulations which are run as a part of the estimation of parameter k 5 using optimization.Thus the simulated responses in X acid and X meth varies with time, see Figure 6.
None of relative sensitivities shown in Table 8 has extreme values.Hence it is concluded that the relative sensitivities do not demand a change of the assumed values of the pertinent a-parameters.

Temperature dependency
In the modified Hill's model (as in the original model) the reaction rates, eqs.( 26) and ( 27), depend on the reactor temperature T reac , cf. ( 28).In the model adaption made in Section 3. Figure 7 shows responses in F meth due to changes in T reac and F f eed for the reactor ADR1.
During the time period shown in Figure 7 both T reac and the feed flow F f eed were changed, but only the variations caused by the temperature change are of interest here.The simulations are based on a mathematical model adapted to ADR1 using the same method for model adaptation as is used for model adaptation to ADR2 as described in Section 3.3.3.
As seen in Figure 7, T reac was increased twice: These temperatures are in the mesophilic range.The changes were implemented as step-wise changes of the temperature setpoint in the temperature controller TC1, cf. Figure 1.
Since the simulated F meth and the real F meth plotted in Figure 7 show similar responses, it can be concluded that eq. ( 28) represents the temperature dependency of the real reactor quite accurately, at least in the mesophilic temperature range.

Reactor temperature model
The reactor temperature T reac is actually a (dynamic) state variable although it is assumed to be a parameter in the modified Hill's model, cf.eq. ( 28).In some situations it is useful to have a dynamic model describing the dynamic behaviour of T reac , e.g. in optimization of the reactor design and operation where energy and economical cost are taken into account, and in modelbased tuning of a temperature controller for the reactor.A dynamic model describing T reac is now derived and adapted the model to the real bioreactor (ADR2).
T reac depends on a number of variables and parameters, e.g. the manipulated supplied power to the electrical heater for the reactor, the temperature of the reactor feed, the feed flow rate, etc.It is reasonable to assume that T reac can be modelled with an energy balance for the liquid of the reactor.Due to the mixing which takes place in the reactor it is assumed that there are homogeneous conditions in the reactor.The energy balance can be written as The liquid in the reactor is assumed having the same thermal characteristics as water.All model parameters in eq. ( 55) except the thermal conductivity G are assumed known, cf.Table 9.
Table 9: Assumed known parameters of reactor temperature model.c = 4200 J/(kg K) ρ = 1000 kg/m 3 V = 250 L variables and parameters in eq. ( 55) are assumed to have SI units, although some are presented in different units in the present section.
Measured time series of all variables in eq. ( 55), i.e.T reac , P heat , F f eed and T room , are available.F f eed is constantly 55 L/d.T room varies between 12.7 o C and 15.8 o C with mean value 13.7 o C. T f eed is assumed to be equal to T room .
G is estimated from experimental data.Three alternative methods are applied: 1. Least squares (LS) method with static model: By assuming that all variables, including T reac , in eq. ( 55) have constant values, dT reac /dt can be set to zero.Solving the resulting static model for G gives from which G can be estimated with the LS method.
2. LS method with dynamic model: The model ( 55) is linear in the unknown parameter G: Here, dT reac /dt is calculated using simple numerical differentiation 2 : G is estimated from eq. ( 57) with the LS method.
3. Nonlinear least squares (NLS) method with dynamic model and additional lag: It is assumed that the original model (55) describes the dynamic properties of T reac .It is also assumed that there may be additional dynamic phenomena due to energy capacitance in the heating element (the coil) and in the reactor wall.Furthermore, a measurement filter in terms of a discrete-time algorithm resembling a time-constant system with time-constant T f = 600 s is actually in operation 2 We used the diff function in MATLAB.
in the AD reactor.This filter also adds dynamics to the temperature behaviour.To obtain a total model able to represent these additional dynamics, the original model ( 55) is augmented with the following general "lag model" incorporating timeconstant dynamics: where θ temp [d] is a time-constant.One motivation for this augmentation is the expectation that the estimation of G in eq. ( 55) will be improved if also θ temp in eq. ( 59) is estimated.Thus, both G and θ temp are estimated.The measured temperature is now represented with T reac lag , while T reac is actually unknown which makes it difficult to apply the ordinary LS method for estimation.Instead, NLS estimation is used based on optimization using iterated simulations.The procedure is the same as was used to estimate parameter k 5 in the modified Hill's model, cf.Section 3.3.3.Therefore, Figure 5 applies, but with the following changes.The measured output is and inputs (time series) are where T room = T f eed .The prediction error, e pred (t k ), is where T reac lag sim is simulated T reac lag .
The model used in the iterated simulations executed by the optimization function consists of the differential equations ( 55) and (59).
The parameter vector to be estimated is The following guessed (initial) parameter values for the estimation are used: The guessed value of G guessed stems from some calculations made with the previous AD reactor (ADR1).10.
To check the quality of the estimates, both SSE (sum of squared error) which is defined by eq. ( 48), and RMSE which is the Root Mean Squared Error index, Varmuza and Filzmoser (2009), are calculated: SSE is the same function as is used to solve the NLS estimation problem, but it does not have the same unit as the prediction error, e. RMSE resembles the standard deviation, and it has the same unit as the preduction error, e.Therefore, RMSE may be more useful than SSE as a performance index for model validation, while SSE is more useful in solving estimation (optimization) problems since it is a square function, Varmuza and Filzmoser (2009).
Figure 8 shows simulated T reac for three cases of parameter estimation, cf.Table 10, together with real (measured) T reac , P heat , and T room = T f eed .(The oscillations in P heat are due to the On/Off control signal from the temperature controller.)The respective time intervals of the timeseries of data used for estimation in the three cases are shown in the upper plot in Figure 8.
Notes to the plots in Figure 8: • Around t = 161.4d and 161.9 d the changes in the measured T reac appears earlier than the changes in P heat assumed to cause the changes in T reac .The sampling time is T s = 15 min = 0.0104 d, which is in the same range as of the periods of the changes of P heater .
• The simulated T reac is clearly larger than the measured T reac between t = 160.6 and 161.0 d.The difference may be due to disturbances for which the controller compensates by increasing the (average) P heater .These disturbances are not due to the observed reduction of T room since the actual values of T room are used in the simulation.We can not identify these disturbances, but they may be unknown variations in T f eed .The results shown in Table 10 and Figure 8 show no large difference between the three estimation methods used here, and all of the simulated responses of T reac resembles quite well the measured T reac .Still, the best result is obtained with the third method (NLS with augmented dynamic model).Hence, the ultimate estimated value of G is selected as G = 1.96 • 10 5 (J/d)K (66) The thermal time-constant, θ thermal , can now be calculated from the model (55).Using the pertinent numerical values for the parameters, In some applications of the thermal model derived in the present section, for example in optimization of reactor design and operation where thermal energy is taken into account, it will be sufficient to use only the main part of the model.The main part is given by eq. ( 55).However, in temperature controller tuning more appropriate controller parameter values can be expected by taking into account also the lag-model (59).

Discussion
The modified Hill's model which has been adapted to the pilot bioreactor fed dairy manure is a relatively simple model compared with alternative models since the model does not contain neither ammonia, alkalinity, nor pH as variables.These variables are more important in reactors fed manure from or poultry because their values may have higher impact on the stability of such reactors.
The modified Hill's model is assumed to be sufficiently accurate as a basis for optimal reactor design and operation, state-estimation and control for a reactor fed dairy manure where the main output is the produced methane gas flow.In applications requiring a prediction of hydrogen or carbondioxide gas production alternative models must be used.
The parameters B 0 (biodegradability constant) and A f (acidity constant) are estimated from data from one experiment only.Ideally, more experimental data should have been used.
In the original Hill's model, Haldane functions are used in the reaction rates.Instead, Monod functions are used, mainly to simplify model adaptation.These simplifications have support in some literature references.
In the orignal Hill's model, the solids residence time and the hydraulic residence time were assumed to be equal.In our study, this assumption caused problems with the model adaptation (results are not shown here), while assuming different residence times, related with a proportionality factor, worked well.
The dynamic model based on energy balance describing the temperature in the liquid phase of the bioreactor assumes homogeneous conditions in the reactor.A model with acceptable accuracy was adapted under this idealized assumption.However, the model adaptation was improved by including an additional timeconstant lag in the model.This lag can be regarded as a representation of inhomogeneous conditions, or spatial variations, in the reactor.

Conclusions
A dynamic model has been adapted to a pilot anaerobic reactor fed dairy manure using steady-state and dynamic operational data.The model is a modification of a model originally developed by Hill (1983).The model is based on material balances, and comprises four state variables, namely biodegradable volatile solids, volatile fatty acids, acidogens, and methanogens.Simulations compared with measured methane gas flow indicate that the model is able to predict the methane gas flow produced in the reactor.
The steady-state data used for the model adaptation are feed flow (loading rate), reactor temperature, methane gas flow, and laboratory analysis values of influent and effluent VS and VFA concentrations at one specific steady-state operating point.The dynamic data used are feed flow, reactor temperature and methane gas flow over a time-interval of 15 days.
Also, a dynamic model for the reactor temperature based on an energy balance of the liquid is adapted to the pilot reactor.The model is able to predict the reactor temperature.A combination of this model and the model of the the anaerobic processes can be useful in optimization of reactor design and operation when energy production and economical costs are taken into account.Furthermore, this model can be used for temperature controller tuning.

Figure 1 :
Figure 1: Piping & Instrumentation Diagram (P&ID)of the biological process line of the pilot plant at Foss Biolab, Skien, Norway.This paper concerns the AD reactor.

Figure 2 Figure 2 :
Figure2shows the steps of the AD process on which Hill's model is based.The alphabetic identificators used in Figure2refer to the model equations presented in the present section.
and k 5 , replace the original parameters (yields) 1/Y , (1 − Y )/Y , 1/Y c , and k meth (1 − Y c )/Y c .Parameters k 1 , k 2 and k 5 are estimated from experimental data, except k 3 which is calculated from the parameter values in the original Hill's model.

Figure 5 :
Figure 5: The estimation of parameter vector p = k 5 is based on optimization (minimization) using iterated simulations.

Figure 6 :
Figure 6: Upper plots: Real values (from laboratory analysis) and simulated values of the state variables S bvs and S vf a together with the respective real concentrations in the feed.Lower plots: Simulated values of the state variables X acid and X meth for which we do not have laboratory analysis data.
3.3 T reac was kept constant at 35 o C. Simulations with modified Hill's model, though not shown here, indicate that a change in T reac gives a dynamic response in F meth .We do not have experimental results with varying T reac for reactor ADR2 which is used in the present study.However, experimental results exist for reactor ADR1 which was in use at Foss Biolab from August 17, 2011 until April 19, 2012.It is fair to assume that the temperature dependency as expressed in Hill's model holds equally well for ADR2 as for ADR1 since the physical appearances of the two reactors are similar and the operation and feed (manure from the same dairy livestock) are similar.

Figure 7 :
Figure 7: Reactor ADR1: Responses in F meth (middle) due to changes in T reac (lower) and F f eed (upper).

Figure 8 :
Figure 8: Upper plot: Real (measured) and simulated T reac , and time intervals for estimation (marked 'o').Middle plot: Real P heat .Lower plot: T room = T f eed .

Table 2 :
Parameters B 0 and A f are found from a laboratory test as described below.Their values are shown in Table2.Parameters assumed to have known values

Table 3 :
Long-term test to find parameter B 0 Time S vsin Biogas prod.t 0 = July 12, 2011 29.25 Non-zero t 1 = Sept. 25, 2011 21.92 Zero It is assumed that a proper value of A f can be found from the approximate steady-state operating point around June 10, 2012 which is used for estimation of the model parameters (except parameter k 5 which is estimated from dynamic responses) using the steady-state version of the dynamic model.The pertinent values needed to calculate A f are shown in Table 4. From eqs.

Table 5 :
Values of inputs and states in the steady-state operation point (t = 66 d; t = 0 = 19.April 2012) used for model adaption.Units are listed in Appendix A.2.

Table 6 :
Values of estimated variables and parameters, and absolute and relative standard deviations (σ).Units are listed in Appendix A.2.

Table 7 :
Standard deviations of quantities which are varied (randomly) in bootstrapping with parametric simulation used to assess variations of estimated parameters.σ Svs in [g/l] σ S bvs [g/l] σ S vf a [g/l]

•
At time t = 60.5 d: From approximately 24 o C to approximately 30 o C.

Table 12 :
Values of inputs and states in one example of a steady-state operation point.F f eed = 45 L/d T reac = 35 o C S vsin = 30.2g/L S bvs = 5.2155 g/L S vf a = 1.0094 g/L X acid = 1.3128 g/L X meth = 0.3635 g/L F meth = 196.1560L CH 4 /d