Aspects on methanogenic biofilm reactor modeling

M. C. Mussati¹, P. Aguirre², M. Fuentes³ and N. Scenna⁴

INGAR-Instituto de Desarrollo y Diseno/CONICET, Avellaneda 3657, (3000) Santa Fe, Argentina.
(¹ mmussati;² paguir;³ mfuentes;⁴ nscenna)@ceride.gov.ar

Abstract — A previously developed deterministic steady state module for modeling methanogenic biofilm reactors has been revised to enlarge the model application range and to deal with system dynamics. Two models for the hydrolysis of non-active biomass representing extreme alternatives: without biomass hydrolysis (model A) and with complete and instantaneous hydrolysis of non-active suspended and attached biomass (model B) were investigated. Both models resulted to be suitable for simulating highly anaerobic loaded systems. However, only model B showed good agreement between experimental and calculated values at low organic loading rates. The values of the specific biofilm detachment rate and the specific microbial death rate of the original model were re-estimated for model B based on a set of step-type disturbances on the organic loading rates. At loading rates ranging between 2 to 4 g COD per day per liter of expanded bed applied to a lab-scale fluidized bed reactor, the parameter estimates were 0.0269 Lg^-1d^-1 and 0.061 d^-1 for the specific biofilm detachment rate and the specific death rate, respectively, with a 95% chi-square confidence level.

Keywords — Anaerobic Digestion. Biofilm Reactor. Process Dynamics. Biomass Hydrolysis.

I. INTRODUCTION

A deterministic dynamic model can be used for different purposes. Through numerical simulation, the dynamical behavior of processes can be predicted, empirical model parameters can be determined, experimental campaigns for discriminating the main mechanisms governing the process dynamics can be defined, and full-scale industrial processes can be optimally designed. Large changes in the environmental conditions experienced under transient behavior such as system start up or operation modes provoked by changes on the influent flow-rate, chemical oxygen demand (COD) and composition, pH and concentration of inhibitory species, may deteriorate the system efficiency and, even worse, cause the process failure. Thus, the development of a dynamic bioreactor model becomes evident to address these needs. Within this context, a mathematical model for the dynamic simulation of anaerobic bioreactors containing biofilm and degrading complex substrates was presented in Mussati et al. (1998). The model allowed simulating different reactor configurations such as single continuous stirred tank reactor (CSTR) and in-series CSTRs to model non-ideal flow pattern. Afterwards, the model was used to investigate the most common loading strategies used for starting up anaerobic wastewater treatment systems (Mussati et al., 2003). Recently, a more rigorous and refined steady state model of an anaerobic methanogenic biofilm reactor-module that accounts for the biological interactions of four microbial groups, ionic equilibrium in solution, gas-liquid transfer phenomena and biofilm processes was presented by Mussati et al. (2005). The model consists of a CSTR type that allocates an inert support material, whose specific surface is taken into account. The biofilm model assumes a homogeneous biofilm of uniform thickness and constant density with no mass transfer resistance. The biofilm detachment process rate is modeled as a second-order function on the biofilm thickness and a first-order function on the mass fraction of the fixed biomass concentration of each microbial group. The balance equations for non-active biomass in liquid and biofilm are included. The model predictions have been satisfactorily compared with steady state experimental data reported in literature from a one-phase methanogenic biofilm system treating an acetic acid-based synthetic effluent and from a two-phase system with combined suspended (acidogenic) and attached (methanogenic) microbial growth treating a food industry waste-water composed by two residual process streams. These satisfactory results were obtained dealing mainly with highly loaded systems.

However, at low COD loading rates the dynamic biofilm reactor model derived from the steady state reactor module failed in predicting the effluent (transient) biomass concentration and total COD for treating a liquid stream composed mostly by acetic acid in an anaerobic fluidized bed reactor. Nevertheless, a good agreement between experimental and predicted values was observed for critical process variables influencing the system performance, such as pH, VFA and biogas production rate. Thus, the focus was on the model of the biological processes where biomass is directly involved; more specifically on the kinetic model and parameter values for the acetoclastic methanogenic step. Sensitivity studies have shown that changes on the model parameter values provide no better predictions of the suspended biomass concentration and total COD at low loading rates. As conclusion, a refinement of the original steady state module model is required to capture more accurately the process dynamics when simulating low loaded anaerobic biofilm reactors.

This paper is then intended to present an alternative to the previously developed deterministic model by improving the modeling of the biological stages and their associated process rates to enlarge the model application range when dealing with system dynamics. The focus is mainly on revising the biofilm system model (more specifically the detachment rate coefficient) and the specific death rate of microorganisms, and considering a hydrolysis model of non-active biomass since the (original) steady state model assumes no hydrolysis of non-active biomass. This paper is organized as follows. In section II, the set-up for the bioreactor used and the experimental protocol followed are described, and the main process variables for monitoring the system performance are depicted. In section III the process model is described emphasizing the two model alternatives investigated: without and with hydrolysis of non-active biomass, named model A (original model) and model B (revised model), respectively. In section IV, the simulation results including the parameter estimation problem are discussed. Finally, the conclusions are drawn in section V.

II. MATERIALS AND METHODS

A. Bioreactor

A lab-scale anaerobic fluidized bed reactor consisting of a 2.0 m high acrylic cylinder with an inner diameter of 0.065 m was used to perform the experiments. The separation compartment placed over the column is a 0.18 m high cylinder with a 0.145 m inner diameter, where gas accumulation and particle sedimentation take place. The effluent discharge, the feed and the recycle suction point are also placed in this compartment. The setup used for bioreactor is schematized in Fig. 1 and the reactor specification data are listed in Table 1.

Figure 1. Experimental set-up for bioreactor used

Table 1. Specification data for bioreactor

B. Inert support material

A sample of sand loaded into the bioreactor was meshed using the Tyler sieve series. The particle size distribution is listed in Table 2. The material specific surface and the surface-volume mean diameter are calculated as in McCabe et al. (1993) and Perry and Chilton (1973). These values and other inert support characteristics are included in Table 3.

Table 2. Results of sand sample meshing

Table 3. Characteristics of the inert support material

C. Analytical Methods

The amount of biogas produced by the bioreactor was measured using a water replacement method. The gas composition was analyzed by a gas chromatograph (Hewlett Packard 6890) equipped with a thermal conductivity detector and a 3 m carbosphere column. Hydrogen was used as carrier gas at 20 mL min^-1. The column was operated at 150 ^oC. The injector and detector temperatures were 100 and 230 ^oC, respectively.

The effluent COD was measured according to the HACH potassium dichromate method approved by USEPA (Cat. 21259-15, 0-1500 ppm). The concentration of the released Cr³⁺ ions was determined by spec-trophotometry (Metrolab 330). The effluent pH was measured using a digital pHmeter (Horiba D-12). Some samples were also measured according to the Standard Methods for COD (APHA, 1995) to evaluate the approximation obtained with the HACH kit. The results were satisfactory.

The acetic, propionic and butyric acid (VFA) concentrations were measured by a high performance liquid chromatograph (HPLC) Hewlett Packard Model Series 1050 equipped with a UV-VIS detector (wavelength: 215 nm). A 20 μL sample volume was injected into a Spherisorb ODS-1 (C18) Classic 5U (250 x 4.6 mm) column (Alltech, Deerfield, IL). The mobile phase consisted of 50% acetonitrile-50% water, sulfuric acid 0.01% (pH 3) at a flow rate of 0.7 ml/min.

D. Experimental Protocol

The responses of the anaerobic fluidized bed reactor to step-type disturbances on the acetic acid concentration [HAc]^T_in and the volumetric feed flow rate Q were investigated. The experiments were carried out at 36±1^oC. The feed consisted initially of 2.0 g COD L^-1 (90% acetate and 10% milk powder plus an amount of sodium bicarbonate for pH adjustment) at 3.6 L d^-1. A step-type disturbance on the inlet acetic acid concentration was the first perturbation introduced (PI) increasing the influent COD with 50% (from 2.0 to 3.0 g L^-1). The feed flow rate was kept constant at 3.6 L d^-1. The disturbance was applied when reactor operated close to steady state conditions corresponding to 2.0 gCOD L^-1 and 3.6 L d^-1. After 28 days under PI conditions, the inlet COD concentration was increased from 3.0 to 4.0 g L^-1 (P2), keeping the same feed flow rate (3.6 L d^-1) during 20 days. Finally, disturbance P3 consisted of a step in the feed flow rate, from 3.6 to 7.2 L d^-1, keeping the input concentration at 4.0 g COD L^-1. The biogas, total COD, pH and acetic acid responses to disturbances applied are depicted in Figs. 2, 3, 4 and 5, respectively.

Figure 2. Biogas response to disturbances P1, P2 and P3

Figure 3. COD response to disturbances P1, P2 and P3

Figure 4. pH response to disturbances P1, P2 and P3

Figure 5. Acetic acid response to disturb. P1, P2 and P3

III. THE PROCESS MODEL

Since the steady state module for modeling anaerobic biofilm reactors presented in Mussati et al. (2005) was developed to simulate the anaerobic digestion process departing from glucose, which is usually considered as a representative contaminant, four anaerobic microbial groups were included: glucose fermenting acidogens, propionic acid degrading acetogens, butyric acid degrading acetogens and acetoclastic methanogens.

As an acetate-based synthetic substrate is used in this work, the acetoclastic methanogens, which metabolize the acetic acid to methane and carbon dioxide, is the only biological step involved. Hereafter, "model A" denotes the resulting dynamic model derived from the steady state module considering only the acetoclastic methanogenic stage. Although the complete steady state reactor model can be found in Mussati et al. (2005), the main model equations are included below for clarity.

Model A accounts for non-active biomass by accumulation of both suspended and attached dead biomass. It is assumed that the non-active biomass does not undergo hydrolysis. Unlike model A, model B assumes complete and instantaneous hydrolysis of non-active suspended and attached biomass. The end product of biomass hydrolysis is assumed to be acetic acid, which is the substrate for the acetoclastic methanogens.

MODEL A: Anaerobic biofilm reactor model: methanogenic stage without hydrolysis of non-active biomass.

Mass balance equations

Suspended methanogens:

(1)

The left-hand side of eq. (1) computes the time accumulation of suspended biomass. The first term of the right-hand side is related to the inlet and outlet mass flow rate of suspended biomass; the second and third terms compute the microbial growth and death rates, respectively; finally, the last term is the biofilm detachment rate.

Acetic acid:

(2)

Carbonate system (inorganic carbon):

(3)

Ammonia-ammonium system:

(4)

Phosphate system:

(5)

"Other anions":

(6)

"Other cations"

(7)

Attached methanogens:

(8)

Non-active susp. biomass:

(9)

Non-active attached biomass:

(10)

CO₂ partial pressure:

(11)

CH₄ partial pressure:

(12)

Expressions for biokinetic constants and inhibition function (Angelidaki et al., 1993; Mussati, 2000)

Maximum specific growth rate of methanogens:

(13)

Specific growth rate of methanogens:

(14)

Half-saturation constant for methanogens:

(15)

Specific death rate of methanogens:

(16)

Growth inhibition function by pH:

(17)

Expressions involved in the gas phase calculation

CO2 gas-liquid mass transfer rate:		(18)
Henry's Law for CO2:		(19)
Total vol. flow rate:		(20)
Vol. CO2 flow rate:		(21)
Vol. CH4 flow rate:		(22)
Vol. CH 4 prod. rate:		(23)
Vol. H2O flow rate:		(24)
Vol. biogas flow rate:		(25)
Gas chamber volume:		(26)

Relationships for the ionic equilibrium in solution

Conc. of dissociated and non-dissociated species:

	(27)
	(28)

Water:

(29)

Charge balance of the system:

(30)

System pH:

(31)

Temperature dependence of physical.chemical constants

Henry's constant for CO₂ (Angelidaki et al., 1993):

(32)

Dissociation constants (Angelidaki et al., 1993):
Carbonic acid:

	(33)
	(34)

Ammonium:

(35)

Water:

(36)

Antoine's equation for water vapor pressure:

(37)

Molar volume:

(38)

MODEL B: Anaerobic biofilm reactor model: methano-genic stage with complete and instantaneous hydrolysis of non-active biomass.

Unlike model A, model B assumes complete and instantaneous hydrolysis of non-active suspended and attached biomass. The end product of biomass hydrolysis is assumed to be acetic acid, which is the substrate for the acetoclastic methanogens. The conversion of non-active biomass to acetic acid is computed on an equivalent COD basis. To address this assumption, the mass balances for the non-active suspended and fixed biomass in model A (Eq. 9 and 10, resp.) are discarded in model B, but a source term is added instead to the acetic acid mass balance equation as below. Thus, models A and B represent extreme situations concerning to the hydrolysis process of non-active biomass.

Acetic acid mass balance equation for model B:

(39)

λ_H is the coefficient for converting non-active biomass to acetic acid on an equivalent COD basis.

IV. RESULTS AND DISCUSSION

A. Simulation results

Figure 6 depicts the dynamic simulation results obtained from models A and B for the same set of parameter values at different organic loading rates (OLRs). It can be observed that both models predict essentially the same steady state total COD values at OLRs higher than 20 g COD L^-1d^-1. However, significant qualitative and quantitative discrepancies between both models with respect to the total COD are observed at lower OLRs.

Figure 6. Total COD predicted by models A and B at different organic loading rates

As mentioned, in Mussati el al. (2005) model A (more specifically its steady state version) was successfully used for simulating the experimental steady state condition for a one-phase methanogenic biofilm system treating an acetic acid-based synthetic effluent (case study 1) and a two-phase system with combined suspended (acidogenic) and attached (methanogenic) mi-crobial growth treating a food industry wastewater composed by two residual process streams (case study 2). The OLRs were 17.0 and 102.5 gCOD L^-1d^-1 for cases 1 and 2, respectively, which are in the application range where both models predict good agreement for the steady state condition. Both models predict almost the same biogas production rates and pH values along the whole OLR range simulated (Fig. 7).

Figure 7. Biogas prod, rate and pH predicted by models A and B at different organic loading rates

Based on these results and on a sensitivity analysis on the biokinetic parameter values for model A, it was concluded that the effluent COD measurements obtained from the lab-scale methanogenic fluidized bed reactor (Fig. 2) are neither qualitatively nor quantitatively properly described by model A at the OLR range experienced (results not shown). Moreover, model parameter estimation using model A resulted in over-predicted total COD values, with the suspended biomass COD being the largest contribution to the total COD. It should be noted that the suspended biomass also comprises the detached active and non-active biomass.

On the other hand, the qualitative transient behavior and steady state predictions obtained using model B, i.e. with complete and instantaneous hydrolysis of non-active biomass, were satisfactory when simulating the experimental results obtained from the lab-scale biore-actor. The preliminary promising results were afterwards confirmed by the model parameter estimation task, which is discussed in the next subsection.

B. Dynamic Model Parameter Estimation

Here, parameter estimation is intended to determine values for the unknown parameters to maximize the probability that the model will predict the values obtained from the experiments. Assuming independent, normally distributed measurement errors ε_lik with zero means and standard deviations σ_lik this maximum likelihood goal is calculated by:

where N represents the total number of measurements taken during all the experiments; θ the set of model parameters to be estimated; L the number of experiments performed; I_l the number of variables measured in the l^th experiment; M_li, the number of measurements of the i^th variable in the l^th experiment; σ_lik the variance of the k^th measurement of variable i in experiment l; _lik and z_lik the k^th measured and k^th predicted values, respectively, of variable i in experiment l.

Since constant and fixed variance for each variable i is assumed, the objective function Φ reduces to a weighted least squares type. Each variable i can be a differential or an algebraic variable. For each experiment l, the duration, the initial condition and the control variables are fixed.

The model was implemented and solved using gPROMS (general PROcess Modeling System; Process Systems Enterprise Ltd., 2004a, 2004b). The gEST tool of gPROMS package was used to solve the maximum likelihood optimization problem (MXLKHD solver). A control vector parameterization algorithm based via single-shooting (CVP_SS) is used for solving the dynamic optimization problems. The algorithm employs the DASOLV code for the solution of the underlying DAE problem and the computation of its sensitivities.

The parameter values for model B are estimated using the transient experimental data obtained from the fluidized bed reactor described in sections II.D and II.A, respectively.

In this case, the (untreated) substrate and suspended biomass are the components that mostly contribute to the total measured COD.

During the previous simulation study (section IV.A), it was observed that the main discrepancies between predicted and measured COD values arise in computing the suspended biomass leaving the system. High sensitivity values of the predicted suspended biomass to the biofilm detachment rate coefficient k_E and the specific microbial death rate b_M were, in turn, observed from a sensitivity analysis carried out using model B. Then, the model parameter estimation consisted of identifying k_E and b_M by minimizing the deviation between the model predictions and measurements of output COD, pH and biogas flow rate corresponding to the three step-type disturbances. The rest of the model parameter values were assumed as in Mussati et al. (2005), the main of which are listed in Table 4.

Table 4. Model parameter values

Figure 8 shows the predicted COD values when estimating k_E and b_M using model B and considering all available experimental data (P1, P2 and P3).

Figure 8. Parameter estimation from disturbances P1, P2, P3

A good agreement between calculated and experimental values is observed for disturbances P1 and P2. The discrepancies related to disturbance P3 can be explained by the short residence time of the liquid, where the validity of the hypothesis of complete and instantaneous hydrolysis of non-active biomass may be questionable. Note that disturbance P3 consisted of increasing the feed flow rate by 100%. However, the steady state value predicted for P3 is acceptable. The parameter estimates resulted to be 0.023 Lg^-1d^-1 and 0.071 d^-1 for k_E and b_M, respectively.

Afterwards, the parameter estimation consisted of discarding the first experimental points of P2 due to an undesirable failure acted on the system, and disturbance P3, except for its first measurements (Fig. 9). In this case, the parameter estimates were 0.0269 Lg^-1d^-1 and 0.061 d^-1 for k_E and b_M, respectively. The 95% chi-square test for goodness of fit (0.05 level of significance) was used to accept or reject the proposed model at long retention time. As the sum of the weighted residual terms (38.0) is less than the 95% chi-square value (58.12) a good fit is obtained. The correlation matrix shows a high correlation of the parameters. The off-diagonal element value of the correlation matrix is -0.7 (Table 5). The main computational statistics reported by gPROMS are listed in Table 6.

using model B

Figure 9. Parameter estimation from disturbances P1, P2 and first days of P3 using model B

Table 5. Parameter estimates for Model B using P1, P2 and first days of P3

Table 6. Computational statistics

Finally, the evolution of the main process variables predicted by model B is shown in Fig. 10, 11 and 12.

Figure 10. Predicted and measured total COD and acetate COD

Figure 11. Predicted and measured biogas prod, rate

Figure 12. Predicted and measured pH

The slightly over-predicted biogas production rates and the under-predicted pH levels are consistent with the assumption of complete and instantaneous hydrolysis of non-active biomass to acetic acid. Indeed, as biomass is hydrolyzed to acetic acid, a higher biogas production rate is predicted from the "excess" of acetic acid, which, in turn, decreases the system pH.

Based on this, with a more refined model for the hydrolysis of non-active biomass, i.e. with a finite hydrolysis rate, the model predictions would better approximate the experimental results. This point will be further addressed.

V. CONCLUSIONS

Difficulties in predicting the process dynamics for a methanogenic biofilm reactor at low organic loading rates using a model derived from a previously developed steady state module, which was successfully applied at loading rates higher than 20 g COD per day per liter, were overcame by revising the hydrolysis model of the non-active biomass and re-estimating the model parameters.

Good steady state and transient model predictions were obtained for disturbances on the inlet substrate concentration at low organic loading rates (2 to 4 g COD per day per liter of expanded bed) hypothesizing complete and instantaneous hydrolysis of non-active biomass to substrate (acetic acid).

The biofilm detachment rate coefficient k_E for the fluidized bed reactor model and the specific death rate b_M were satisfactorily estimated with a 95% chi-square confidence level. The correlation matrix shows a high correlation between them. For disturbances on the feed flow rate at higher loading rates (8 g COD per day per liter) acceptable steady state predictions were obtained. In that case, the dynamic discrepancies between predicted values and measurements may be explained by the questionable validity of the hypothesis of complete and instantaneous hydrolysis of non-active biomass. This is supported by the good results obtained when only first measurements of P3 were considered.

The model parameter estimation from dynamic (transient) data is helpful since they allow obtaining more information about the system behavior. This is so because they correctly predict not only the values of the steady state corresponding to all the disturbances applied but also the dynamics. This fact warned about the need for a modification of the biomass hydrolysis model. However, a critical aspect in using dynamic models for parameter estimation in continuous anaerobic bioreactors for wastewater treatment should be pointed out: the long period required for experimentation and, consequently, the increased risk of undesirable failure of the experimental setup. This may lead to poor parameter estimates if the failure is not appropriately modeled and included into the estimation process.

A further refinement of the model of the non-active biomass hydrolysis, i.e. with a finite hydrolysis rate, is required for handling cases in between the two extremes, i.e. with and without hydrolysis.

The dynamic model can be used for process design and optimization of the start up and operation of anaerobic methanogenic biofilm systems for a wider application range than the previously developed models.

NOMENCLATURE

A, B, C, D: Coefficients in correlations.
A^-: Anions different to OH^-, Ac^-, Pr^-, But^- ions and ionic species of carbonate and phosphate systems
C⁺: Cations different to hydrogen and ammonium ions.
CO_2d: Dissolved carbon dioxide.
H: Henry's Law constant
K_T: Gas-liquid mass transfer coefficient (d^-1)
k_E: Biofilm detachment rate coefficient (L g^-1d^-1)
k_E^*: Biofilm detachment rate coefficient independent of the support material (g d^-1 L^-1 cm^-2)
k: Specific substrate utilization rate (g sub. gVSS^-1 d^-1)
K₁ K₂, K₃, K_j: First, second, third and jth acid dissociation constant, respectively.
K_W: Water dissociation constant
K_i: Inhibition constant (mol L^-1)
K_S: Monod's saturation constant (mol L^-1)
NH_3D: Dissolved ammonia
P_biogas: Volumetric biogas flow rate (L d^-1)
P^*: Volumetric production rate (mol L^-1 d^-1)
p, P_T, P_V: Partial pressure, total gas pressure, vapor pressure, resp. (atm).
pK_h, pK₁: Upper and lower pH-coefficients, resp., in the pH inhibition function
Q: Volumetric flow rate (L d"¹)
sv^o: Standard molar volume (22.4 L mol^-1 at 25 C)
t: Time (d)
T: Temperature (C)
V: Liquid volume or reactor volume (L)
X^F_T: Total (active plus non-active) biofilm conc.(g L^-1)
X^T_T: Total (active/non-active and suspended/fixed) biomass cone, (g L^-1)
Y: Yield coef., g biomass/mol of chemical specie
τ: Retention time (d)

Symbols

[ ] : Concentration (mol L^-1)

Sub-indexes / Supra-indexes

na: Non-active biomass
F: Biofilm
M: Methanogenic stage; methanogens
m: Number of protons released.
n: Total number of protons
S: Suspended
T: Total

ACKNOWLEDGEMENTS
The financial support from CON1CET, ANPCyT and UNL of Argentina is acknowledged.

REFERENCES
1. Angelidaki, I., L. Ellegaard and B.K. Ahring, "A mathematical model for dynamic simulation of anaerobic digestion of complex substrates: focusing on ammonia inhibition", Biotechnol. Bioeng., 42, 159-166 (1993).         [ Links ]
2. APHA. Standard Methods for the Examination of Water and Wastewater. 19th Edition. American Public Health Association (1995).         [ Links ]
3. McCabe, W.L., J.C. Smith, P. Harriott, Unit Operations of Chemical Engineering. New York: McGraw-Hill Book Co (5th ed) (1993).         [ Links ]
4. Mussati, M., P. Aguirre and N. Scenna, "Modeling of real biological reactors for the treatment of complex substrates. Dynamic simulation." Computers Chem. Engng, 22, S723-S726 (1998).         [ Links ]
5. Mussati, M.C. "Modeling of anaerobic biofilm reactors. Application to wastewater treatment systems". Ph.D. dissertation (in Spanish), Universidad Nacional del Litoral, Argentina (2000).         [ Links ]
6. Mussati, M., M. Fuentes, P. Aguirre and N. Scenna, "Optimizing load policy in anaerobic biofilm reactors for wastewater treatment." Latin American Applied Research, 33, 301-305 (2003).         [ Links ]
7. Mussati, M.C., M. Fuentes, P. Aguirre and N. Scenna, "A steady-state module for modeling anaerobic biofilm reactors." Latin American Applied Research, 35, 255-263 (2005).         [ Links ]
8. Perry, R.H., C.H. Chilton, Chemical engineers' handbook. New York: McGraw-Hill (5th ed). (1973).         [ Links ]
9. Process Systems Enterprise Ltd., gPROMS Advanced User Guide, Release 2.3, London (2004a).         [ Links ]
10. Process Systems Enterprise Ltd., gPROMS Introductory User Guide, Release 2.3.1, London (2004b).         [ Links ]

Received: May 11, 2005.
Accepted: February 7, 2006.
Recommended by Subject Editor J. M. Pinto.