Not Just Numbers: Mathematical Modelling and Its Contribution to Anaerobic Digestion Processes

Abstract

Mathematical modelling of bioprocesses has a long and notable history, with eminent contributions from fields including microbiology, ecology, biophysics, chemistry, statistics, control theory and mathematical theory. This richness of ideas and breadth of concepts provide great motivation for inquisitive engineers and intrepid scientists to try their hand at modelling, and this collaboration of disciplines has also delivered significant milestones in the quality and application of models for both theoretical and practical interrogation of engineered biological systems. The focus of this review is the anaerobic digestion process, which, as a technology that has come in and out of fashion, remains a fundamental process for addressing the global climate emergency. Whether with conventional anaerobic digestion systems, biorefineries, or other anaerobic technologies, mathematical models are important tools that are used to design, monitor, control and optimise the process. Both highly structured, mechanistic models and data-driven approaches have been used extensively over half a decade, but recent advances in computational capacity, scientific understanding and diversity and quality of process data, presents an opportunity for the development of new modelling paradigms, augmentation of existing methods, or even incorporation of tools from other disciplines, to ensure that anaerobic digestion research can remain resilient and relevant in the face of emerging and future challenges.

1. Introduction

Mathematical modelling has an important and continued role to play in the design and operation of engineered biological systems (EBS), where the inherently complex interaction between biotic and abiotic components may result in unpredictable and undesirable process behaviour requiring operator intervention or advanced control architecture. Over the last four decades (circa 1980 onwards), much focus has been on the development and extension of classical models based on first-principles theory and engineering knowledge. Notably, the work of several research groups coalesced around a unified mechanistic model with the aim to represent, i.e. simulate, the dynamics of the activated sludge treatment process [1,2]. This reference model, the Activated Sludge Model No. 1 (ASM1) [3], initiated a highly productive period in wastewater engineering, where mathematical modelling became a distinct and important field within the bioprocess engineering community. Iterations of ASM1 to include phosphate accumulating organisms (ASM2 [4]) and subsequent consideration of their denitrifying ability (ASM2d [5]), facilitated a greater degree of granularity in the modelling of wastewater treatment processes. However, the increased model complexity due to additional requirements on parameterisation and model calibration highlighted the difficulty in developing rigorous mathematical models that are a trade-off between complexity and accuracy. Nevertheless, the framework of these models has remained in popular use, through refinement and augmentation, as the standard modelling tool used by wastewater engineering practitioners and researchers.

The focus of this review is the globally used anaerobic digestion (AD) process that conventionally transforms biodegradable organics (e.g., agricultural residues, food waste, wastewater treatment sludge) into energy-rich biomethane and nutrient-rich biosolids under anaerobic environmental conditions. While the standard process consists of four well-characterised conversion stages (hydrolysis → acidogenesis → acetogensis → methanogesis), the complexity that manifests in AD biochemistry, energetics, microbial diversity and process operation lends itself well to mathematical modelling, as will be explored here.

Given the relative ubiquity of the ASMs and a reawakening of AD research during the late 1990s, it was inevitable that engineers and scientists would seek to develop a similar modelling framework to facilitate dynamical interrogation of this net energy positive environmental technology. Building on the fundamental process kinetics and physico-chemistry of the ASMs, with some changes in reporting of units and use of kinetic rates, Anaerobic Digestion Model No. 1 (ADM1) [6] provided a valuable in silico tool to support an expanding field of academic research. ADM1 has remained the standard model in AD research and in its integration with ASM1 as the Benchmark Simulation Model No. 2 [7,8], has provided the ability for plant-wide modelling of wastewater treatment systems [9]. This process coupling extends the capacity of sub-component models for simulation, design and optimisation, to include whole plant Life Cycle Assessment [10], future-proofing the transition from conventional wastewater treatment plants to water resource recovery facilities [11,12], and the augmentation of novel processes and technologies [13,14].

ADM1 is known to have severe limitations. As with ASM1, it does not handle phosphorus and its associated biochemistry, nor does it consider sulphur and iron species or, until recently, the bioavailability of trace elements, which are important components for stable anaerobic digester operation. Extensions to the original model have subsequently addressed these omissions [15,16,17,18]. However, unlike the iterations of the ASMs, there has been no standardised progression of ADM1, but rather bespoke modifications and simplifications tailored towards process specific objectives [19,20], more efficient computation in engineering software [9,21], or the incorporation of new empirical knowledge and parameter data [22,23]. In their excellent 2015 review of AD modelling application and needs, Batstone et al. [22] state that up to that point, forming consensus towards the establishment of a standardised ADM2 is an attractive proposition, but establishment of both the science and underlying mechanics of the model remained transitory.

In this review, we revisit the origins of mathematical modelling for application to anaerobic digestion and comment on the current state-of-the-art with respect to the theory, applications, and technologies that comprise or exist alongside models for research and practice. In Section 2, we place the most widely known and used model, ADM1, in a wider historical context and discuss its relevance to current practice. Large, highly parameterised models have well-understood limitations and requirements for managing uncertainty. Although exponential increases in computational performance (cf. Moore’s Law) have significantly reduced the simulation time of complex and highly dimensional models, simplifying assumptions and lumping of parameters are often employed in mathematical models to aid study of anaerobic systems. However, in Section 3 we discuss the potential insights that reduced or minimal models, which often represent a sub-component of the larger process, can provide, particularly in combination with analytical tools common to a branch of mathematics known as dynamical systems theory. Deterministic modelling remains generally the standard in engineering practice. Although the use of data-driven models and machine learning for monitoring and control of biological systems has developed in a similar time-frame, it has largely remained under-used by the engineering community. In Section 4, the growth of digital tools and data modelling for application to anaerobic digestion is described and comment on the increasing need for making use of the increasingly diverse sources of data and meta-data that have potential to reveal process knowledge ex-situ and across process scales. In Section 5 we conclude by taking a wide-lens view of the future of mathematical modelling for anaerobic digestion, focusing on the role of models as anaerobic processes diversify and metabolic pathways are re-purposed for alternative energetic and chemical products.

Engineered anaerobic systems have a vital function for future energy security, and the use of mathematical models will continue to be important and much more intrinsic in plant design, operation and control. This review celebrates what has gone before and the increasing role that mathematical modelling plays in design, control, operation and analysis of engineered anaerobic systems (See Figure 1 for a chronological chart of published studies on AD modelling). The review provides some commentary on anticipated future trends, as greater empirical and theoretical knowledge becomes available for derivation of more representative and rigorous models, often by incorporating ideas employed more commonly in other fields. The review aims to be accessible to readers from a diverse background and, as such, we limit the use of mathematical equations to elementary examples that are qualified in the text. More thorough reviews of the mathematical derivations of the methods presented here are cited in the text.

2. A First-Principles, Engineering Approach to Modelling Anaerobic Digestion

The history of anaerobic digestion modelling from an engineering perspective, which is a mathematical description of the microbial dynamics and biochemical processes resulting in sequential biodegradation of complex organic material into biogas, begins over 50 years ago.

The pioneering work of John F. Andrews, an environmental engineer at Clemson University, South Carolina, formed the touchstone for incremental development of practical mathematical models that could be understood and used by both researchers and practitioners to augment their experimental studies. Andrews used the structural elements of Ordinary Differential Equations (ODEs), which had facilitated the dynamical study of microbiology using mathematical constructs by researchers such as Lotka, Volterra, Haldane, Monod, Rashevsky, the forefathers of a discipline now commonly referred to as mathematical biology [24]. It is interesting to note that these key players came from broad scientific disciplines; Haldane was a geneticist, Monod a biochemist, Volterra a mathematician, Lotka a biophysicist, and Rashevsky a theoretical physicist.

The early work by Andrews built on the empirical and theoretic research of these pioneers to derive mathematical expressions that were founded on first-principles chemistry encapsulated by the Michaelis-Menten enzymatic reaction rate Equation [25]

$$\dot{P}=\frac{V_{max}\left[S\right]}{K_M+\left[S\right]},$$

(1)

where $\dot{P}$ is the rate of product formation, $V_{max}$ is the maximum achievable rate of the chemical system, $K_M$ is a constant representing the substrate concentration, $\left[S\right]$, at a reaction rate half the $V_{max}$. Through the early work of Lotka and Volterra, describing the dynamics of biological systems as a pair of non-linear ODEs (i.e., the predator-prey equations) (see, for example [26]), and the empirical rationalisation of the Michaelis-Menten equation by Jacques Monod to derive an equation for microbial growth on a rate-limiting nutrient [27], the platform for mathematical modelling of bioengineering systems was established [28]. Andrews synthesised this structural and empirical knowledge to describe microorganisms growing on substrates that are self-inhibitory above a threshold concentration (e.g., those found in industrial wastewater treatment such as ammonia, phenols, and volatile fatty acids). Andrews coupled the Haldane function describing enzymatic inhibition with the Monod equation for microbial growth to develop Lotka-Volterra type non-linear ODEs for simulating the continuous-culture dynamics of microorganisms under rate-limiting inhibition [29]

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{X}_1}{\mathrm{d}t}& =\frac{\mathrm{d}{X}_0}{\theta }-\frac{\mathrm{d}{X}_1}{\theta }+\left[\frac{\widehat{\mu }}{1+K_s/S_1+S_1/K_i}\right]X_1,\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill \frac{\mathrm{d}{S}_1}{\mathrm{d}t}& =\frac{\mathrm{d}{S}_0}{\theta }-\frac{\mathrm{d}{S}_1}{\theta }-\frac{1}{Y}\left[\frac{\widehat{\mu }}{1+K_s/S_1+S_1/K_i}\right]X_1,\hfill \end{array}$$

(3)

where $X_i,i\in [0,1]$ are the microbial biomass concentrations in the reactor influent and effluent, respectively, $S_i,i\in [0,1]$ are the substrate concentrations in the reactor influent and effluent, $\widehat{\mu }$ is the maximum specific growth rate, $\theta$ is the mean residence time in the reactor (Note. 1/$\theta$ is commonly denoted as D, the dilution rate, in current modelling practice), $K_S$ is the half-saturation constant, $K_i$ is the inhibition constant, and Y is the yield coefficient of biomass growing on the substrate (for units, see original article [29]). What is notable about this study is that it is one of the first examples of the use of digital computing to obtain solutions to a system of equations describing microbial dynamics. A digital analogue simulation program, PACTOLUS, ran the model on an IBM360 computer, employing a circuit-like block-diagram approach more commonly associated with electronic engineering.

The importance of seminal work in the late 1960s and early 1970s by researchers such as Andrews and Perry L. McCarty, a pioneering environmental engineer whose work on AD kinetics and energetics resonates with engineers and modelers to this day, cannot be understated. They developed an empirically derived basis for coupling microbiology with engineering practice using mathematical tenets. Integral to this was the understanding that process stability, and ultimately failure (e.g., through organic overloading, product inhibition), could be predicted by appropriate kinetic models supported by experimental evidence. Andrews published the first mathematical model of the anaerobic digestion process in 1969 to address the increasing reports of digester failures [30]. Similarly, McCarty and his former PhD student, Lawrence, derived process models from observations from anaerobic continuous cultures to determine limits of operation based on kinetic theory [31,32]. Anaerobic digestion as an industrial process became prominent in Europe in response to energy shortages during World War II, and in the United States, anaerobic municipal solid waste treatment had existed as a technology from around the same period. Given the complexity of feedstocks that were being used in AD, and the relative immaturity in process engineering, the high prevalence of reactors experiencing process instability and failing would be significant. It was, therefore, serendipitous that Andrews initiated a first epoch in AD process modelling to understand stability just a few years before the 1973 OPEC oil embargo precipitated an oil crisis that directly led to a drive towards alternative energy supplies, with AD at the forefront.

The dynamic modelling approach was extended further to include a method to incorporate changes in pH through the introduction of phase (gas, liquid, biological) interactions and alkaline buffering capacity (i.e., CO${}_2$—bicarbonate system) [33]. The time-dependent, predictive characteristics of this extended model provided the ability to obtain a semi-quantitative assessment of the model uncertainty, which is the difference between observed (experimental) data and model prediction for the model variables of interest. Noting that model uncertainty comprises parametric uncertainty and structural uncertainty, and both model inputs and experimental observations are subject to measurement uncertainty, then it is generally accepted that the model predictions are not a perfect representation of reality [34]. In turn, this acknowledgment of uncertainty laid the foundation for model-based process control of bioengineering systems such as the activated sludge process and anaerobic digestion [34,35,36,37,38]. This was a momentous leap-forward in the early days of programmatic computing, especially as industry and academia were starting to coalesce around methods for incorporation of Instrumentation, Control and Automation (ICA) technology for the improvement of biological treatment processes [39,40].

Over subsequent years, other research groups began to build on the work of Andrews and further extend AD models to include additional stages of the process beyond acetoclastic methanogensis, the focus of the earlier models [30,33]. Fundamentally, these later models adopt the same mass-balanced, reaction kinetics approach as Andrews, adding more processes and components to describe an increasingly more complex but also more complete conventional anaerobic digestion process. Hill and Barth included heterotrophic acidogens and the obligate methanogens previously studied by Andrews, and a function describing methanogen inhibition by Volatile Fatty Acids (VFA) and ammonia [41]. As with previous models, the authors considered only the acetate metabolic pathway based on existing empirical knowledge at the time. Notation in this model differs minimally from that provided by Andrews in Equations (2) and (3), with the dilution term ($1/\theta$) now given as the flow rate over the reactor volume ($F/V$) and the differential terms on the right-hand side stated implicitly.

A thorough mathematical description of the gas-liquid transfer of carbonate and ammonia, and the algebraic equations describing charge balancing to calculate temperature-corrected pH, is reported for the first time in this work. This detail invariably led to greater model complexity and, at the time, reliance on high specification computers and bespoke programming languages (Continuous modelling Simulation Program, CMSP, op. cit. Hill & Barth). Although the anaerobic digestion process was understood to be microbially diverse, process or population scale modelling could be simplified due to the highly structured and often compartmentalised nature of bioreactor systems. This structuring, based on agglomeration of microbial ecological functions, i.e., step-wise mass-balanced and stoichiometrically rational biochemical reactions performed by mixed but functionally distinct microorganisms, provided a firm basis from which mathematical modelling of biologically engineered systems became established in the following two decades.

Table 1. Important AD mathematical modelling milestones towards ADM1. Key (Process): EH = Enzymatic hydrolysis, AG = Acidogenesis, AGp = propionate Acidogenesis, AGb = butyrate Acidogenesis, AGa = acetate Acidogenesis, AcG = Acetogenesis, pAcG = propionate-degrading Acetogenesis, bAcG = butyrate-degrading Acetogenesis, fAcG = fatty acid-degrading Acetogenesis, AMG = Acetoclastic Methanogenesis, HMG = Hydrogenotrophic Methanogenesis; (Industry): ${}^A$ = Agricultural, ${}^{BW}$ Biological (Activated) Sludge, ${}^G$ General AD, ${}^{IW}$ Industrial Wastewater; (Program): CSMP = Continuous Systems modelling Program, NR = Dynamic solver/program not reported. Note. Number of parameters include conventional operating/control parameters such as dilution rate and influent concentrations. * Parameter information not available.

Model	# State Variables	# Parameters	Year	Program	Notes	Reference
[AGa →AMG] ${}^A$	9	28	1977	CSMP	VFA and NH${}_3^{+}$ inhibition on MG only	[41]
[AGa →AMG + HMG] ${}^A$	4	N/A *	1982	CSMP	Complete use of H${}_2$ assumed	[45,46]
[AGp + AGb →AcG→AMG + HMG] ${}^G$	9	18	1983	NR	pH regulation through redox potential by H${}_2$	[47]
[EH →AG→pAcG→AMG + HMG] ${}^G$	10	35	1984	MIMIC	First generalised and structured model (stoichiometric matrix)	[48]
[EH →AGa→AMG] ${}^{BW}$	7	14	1986	NR	Hydrolysis of dead, particulate biomass is rate-limiting step	[49]
[AGa →AMG] ${}^{IW}$	5	23	1986	NR	VFA inhibition on AG & MG, NH${}_3^{+}$ inhibition on MG	[50]
[EH → AG → pAcG + bAcG →AMG] ${}^A$	16	33	1993	Pascal	NH${}_3^{+}$ inhibition & pH-Temperature interactions	[51]
[EH →AG→pAcG + fAcG→AMG + HMG] ${}^{BW}$	18	48	1993	Pascal	Extension of structured approach. Focus on acetate kinetics/digester stability	[52]
[EH →AGp + AGa→AcG + pSR + aSR→AMG + HMG] ${}^{IW}$	33	55	1994	METHANE	H${}_2$ inhibition considered	[53]
[EH →AG→MG] ${}^G$	5	14	1996	TUTSIM	Reduced order for process control	[38]

indicates a selection of the key mathematical models that provided necessary iterations towards a consensus dynamic model for anaerobic digestion, manifest as the Activated Sludge Model No. 1. For brevity, the list is by no means exhaustive and does not aim to present any detailed analysis of the models, which may instead be found in comprehensive reviews elsewhere [42,43,44].

Modelling efforts were bolstered by research conducted in empirical microbiology, which contributed to providing better characterisation (i.e., parameterisation) of anaerobic growth kinetics and understanding of the metabolic and biochemistry of the principal reaction pathways [54,55,56,57]. Experimental work provided the data and knowledge to determine the structure and stoichiometric properties of the anaerobic digestion process that, in turn, allowed for an emerging formalism in parallel model development. An important body of work that motivated a structural approach to AD modelling was Gujer and Zehnder’s extensive study of anaerobic digestion conversion processes [55], which built on an earlier proposed structure describing overall substrate flux from particulate organic matter to methane and carbon dioxide end-products [58]. Their study describes anaerobic digestion as six distinct conversion processes, integrating established knowledge of the acetogenesis and methanogenesis steps with emerging bioprocess kinetics and microbial research on the mechanisms of organic matter hydrolysis and acidogenesis, thus providing a complete description of the canonical AD process as understood at that time. Defining the units of substrate flux in mass of theoretical Chemical Oxygen Demand (tCOD) as opposed to a molar basis, was a consideration replicated from activated sludge modelling, in its conventional use for measurement of wastewater characteristics and given its implicit function for calculating carbon oxidation states and, thus mass balancing of biochemical reactions [6,59]. Note. Typically, for anaerobic modelling, concentrations of organic compounds are given in units kgCOD m${}^{-3}$, and inorganic carbon and nitrogen in kmole m${}^{-3}$, whereas ASMs generally use mg L${}^{-1}$. The case for using molar units, however, was argued on the basis that the inclusion of inorganic carbon (CO${}_2$) to close the mass balance for anaerobic oxidation of VFAs invalidates the stoichiometric carbon balance of the biochemical reactions [60].

The structured approach was further formalised by Bryers [48] using simple matrix algebra to derive expressions for the biological reactions in the form of a stoichiometric matrix of ($m\times n$) molar coefficients (${\nu }_{i,j}$, for $i=1,\dots m$ reactions and $j=1,\dots n$ components [reactants, biomass, gases]) to ensure a completely mass-balanced system of process equations. In considering a unit basis of mass COD, normalised yield coefficients defining $S_j$ consumed (−) or formed (+) in the ith reaction are given by

$$Y_{i,j}=\frac{{\nu }_{i,j}{M}_j{m}_j}{{\nu }_{i,lim}{M}_{lim}{m}_{lim}},$$

(4)

where $M_j$ and $m_j$ are the molecular weight and COD mass equivalence of component $S_j$, respectively, and ${}_{lim}$ indicates the limiting reactant in the ith reaction. Now, using matrix algebra, we have the following condition under the principle of mass conservation

$$\mathbf{{\rm Y}}·S=0,$$

(5)

where Y is an $m\times n$ matrix of yield coefficients and S is a vector of length n of components.

The use of a structural matrix for deriving a system of equations representing biochemical reactions was first proposed by Petersen [61] and employed in the development of a mathematical model of activated sludge systems by the IAWPRC (now International Water Association, IWA) Task Group on Mathematical modelling for Design and Operation of Biological Wastewater Treatment [3]. The Gujer matrix, as schematically summarised in Figure 2, is partitioned into three sections comprising a matrix of stoichiometric coefficients and a process rate vector, as described above, and a composition matrix that includes conservative constituents of the state variables, e.g., biomass contains both COD and nitrogen, which ensure proper mass and charge balancing as model resolution increases, i.e., the product of each row in $\mathbf{{\rm Y}}$ and each row in the composition matrix will be zero.

The conversion rates for each state variable can be simply calculated by the column-wise summation of the product of stoichiometric terms and process rates

$$\dot{R}={\rho }^T·\mathsf{\Lambda }=\sum _{j=1}^n{\rho }_j{\nu }_{i,j},i=1,\dots ,m,$$

(6)

where $\Lambda$ is the matrix of $m\times n$ stoichiometric terms and $\rho$ is a vector of n process specific kinetic rates, which take the form

$${\rho }_j=\frac{{\mu }_{j,lim}X}{Y_{j,lim}},$$

(7)

where ${\mu }_{j,lim}$ is the substrate-limited growth rate of organism X, defined by some growth function (e.g., Monod saturating kinetics, Contois density dependent growth, etc.), cf. the final term in Equation (3). Note. Hydrolysis in these early models is typically included as a first-order reaction rate (i.e., $r_{hyd}=k_{hyd}S$). The notation $\dot{R}\equiv \frac{\mathrm{d}{R}_j}{\mathrm{d}t},\mathrm{for}j=1,\dots ,n$, is used to indicate that the resultant vector of rates are time-dependent and equivalent to a system of ODEs, with additional terms to mass balance the reactor and describe any liquid-gas transfer, which can be integrated to give transient dynamics of the process with respect to each state variable.

This highly structured matrix formulation for constructing biochemical reaction models was adopted as standard within the engineering community, leading to a general framework for modelling the activated sludge process, specifically related to nitrogen removal (ASM1) [3]). It was particularly open to integration in simulation software, especially programs with superior abilities for rapidly solving matrix linear algebra problems, such as Matlab and FORTRAN. Naturally, this led to several bespoke simulation software to provide engineers with the tools to rapidly develop and interrogate the models through process simulation, often in a modular way (e.g., Biowin, STOAT, GPS-X, SIMBA, AQUASIM) [1]. Although the Anaerobic Digestion Model No. 1 was published over ten years after the first ASM, work on developing a general anaerobic digestion model employing a similar structural and notational form emerged earlier. Siegrist et al. used a Gujer matrix of 14 processes and 18 state variables to model the digestion of mesophilic sewage sludge [52] and Christ et al. used this approach with a better description of the hydrolysis constants for the digestion of municipal organic residues [62]. Other groups incorporated greater degrees of resolution in existing deterministic models of the process, for example by separation of the complex organic substrates (carbohydrates and proteins) for hydrolysis [51,63,64].

The Anaerobic Digestion Model No. 1, developed by the IWA Task Group for Mathematical modelling of Anaerobic Digestion Processes [6], consolidates much of the preceding theory, empirical knowledge and mathematical structure referenced in this section. The original model represents a comprehensive five-stage anaerobic digestion process (including an initial disintegration step), with 19 sub-processes, comprising biochemical reactions (defined through the use of the Gujer matrix) incorporating degradable substrates and soluble/insoluble inerts, and physico-chemical reactions (non-biologically mediated processes such as ion association/disassociation, phase transfer). Biochemical processes are (1) an initial disintegration of composite material (e.g., lysed biomass) to particulates, which undergo (2) enzymatic hydrolysis (forming carbohydrates, proteins and lipids), (3) acidogenesis of simple sugars and amino acids to organic acids, CO${}_2$ and H${}_2$, (4) acetogenesis of long-chain fatty acids and organic acids (butyrate, valerate, propionate) to acetate and hydrogen, and finally (5) acetoclastic and hydrogenotrophic methanogenesis. The model may be solved as a system of 32 ODEs, where some lumped variables are decoupled (e.g., inorganic carbon $S_{\mathrm{ic}}\equiv [S_{\mathrm{CO}2},S_{\mathrm{HCO}3-}]$). Biochemical reactions are assumed to follow the 1^st-order rate law for extra-cellular processes and Monod-type substrate uptake for intra-cellular kinetics. Alternatively, the model can be defined as a system of Differential Algebraic Equations (DAE), with the biochemical reactions integrated explicitly as a system of 26 ODEs (7 groups of bacteria and archaea, 19 catalysing processes), three liquid phase physico-chemical processes, and eight implicit algebraic variables (per reactor). This approach is a means to overcome the inherent ‘stiffness’ problem due to the fast hydrogen and pH dynamics (now explicitly calculated algebraically), which results in some numerical solvers encountering difficulty in solving the system of ODEs.

As noted earlier, the model is far from perfect but its flexibility and structure means it is suitable as a tool to be used at the design stage of AD plant development, for qualitative scenario testing, or for benchmarking control schemes [65,66,67]. As with the ASMs, modifications to the model have been implemented based on process requirements, such as operating mode [68,69], or to include new biochemical processes [70,71]. A more comprehensive transition to a consensus ADM No. 2 that addresses the main issues with the original model has yet to emerge, perhaps due to the lack of clear consensus or motivation, or that its inherent structure provides ample opportunity for modification of components for bespoke applications. Many of the suggested approaches relate to the lack of identifiability and heavy parameterisation of ADM1, which makes its use for process monitoring and model-based control difficult [72]. It is interesting to note that models that have sought to overcome the structural verbosity are greatly simplified, relying on fewer state variables and lumping parameters where practical. These models, such as the widely studied Bernard model (sometimes referred to as AMOCO or AM2) [73] and its subsequent mathematical treatments [43,74,75,76], relate back to the pioneering early work of Andrews and Graef, allow the process to be identifiable and enable process control schemes to be designed for improved or optimal process performance. Given the substantial quantity of literature on the topic of model identifiability, i.e., model selection and parameterisation, we refer the readers to an excellent review of the subject related to anaerobic digestion [77].

In the next section, we take a closer look at how reduced modelling approaches and their mathematical analysis are being used to provide further insights into the dynamic and steady-state behaviour of anaerobic digestion.

3. Reduced-Order Models Unlock the Power of Mathematics

It is common practice in modelling to start with a simple description of the system under investigation before building in complexity to achieve a more complete representation of reality. However, it should be noted that greater complexity does not necessarily equate to a higher degree of accuracy [78]. Terminology such as simplicity and complexity are often abstract and subjective, but in the frame of engineering models such as those discussed here, models with a greater number of parameters, state variables (dimension), and processes, or those that are highly non-linear, can be defined as complex. In mathematical sciences, this complexity is often referred to as model order and is the principle driver for the richness of dynamical behaviours captured in microbial community models [79].

The early anaerobic digestion models were low-order and aggregated, developed with the empirical knowledge of the biochemical reactions at that time, and constrained by the available computing resources and numerical solvers for integrating the system of differential equations. The dramatic advances in computational power and efficiency have provided engineers with an ever-increasing ability to simulate bioprocesses at greater levels of complexity and speed. However, there are still significant trade-offs between model resolution, accuracy and performance. For example, consider a mathematical model describing biofilm growth in a biological reactor. There are several approaches that can be taken to represent the evolution of this process in time, where the simplest would be to assume a homogeneous, quasi steady-state biofilm and to develop a model as a system of ODEs, for which the solution is trivial (e.g., [80]). The underlying assumption of a completely mixed system is not correct in practice, as biofilms are subject to diffusive mass transfer limitation, requiring consideration of spatial heterogeneity. Describing biofilm processes at higher dimensions, using Partial Differential Equations (PDEs) to solve the diffusion-limited mass transfer of biomass, gas and liquid (i.e., using the second-order Navier–Stokes equations or via the simpler boundary layer approximation), is non-trivial and significantly increases the computational requirements to find solutions [81]. A trade-off between model complexity and empirical uncertainty, which is the extent and accuracy of data to support model parameterisation, is one of the fundamental challenges that limits anaerobic digestion modelling [82].

Model order reduction methods such as decomposition, interpolation, low-order approximation and simplifying physics, are generally applied in engineering when high order models are computationally costly, or simulation become impractical. In EBS, mathematical biologists are often interested in the qualitative behaviour of the system rather than quantitative predictions, and the simpler models often allow for a tractable and complete analysis using techniques re-purposed from dynamical systems theory. It is important to ensure that the simplifying process retains conservation of the properties of the full-order model and has a low approximation error, otherwise analysis of the simple model becomes a study of a largely theoretical system with no further insight into the actual system. Anaerobic digestion is an ideal process for reduced-order modelling as it is highly compartmentalised allowing for easy decoupling of system interactions into smaller sub-processes with fewer state variables and parameters. Physico-chemical processes that may operate at significantly different timescales to the biochemical reactions may also be excluded.

The work of Andrews and colleagues (see Section 2) laid the early foundations for the application of dynamical analysis to anaerobic systems. They presented not only a template for mathematical modelling and dynamic simulation of the process, but were the first to discuss the opportunities for automatic process control to address system instability or failure caused by acidification or functional inhibition [33,36]. The analysis at this time involved exhaustive and time-consuming numerical simulation (by today’s standards) of a model over hundreds of operating points and process conditions. Although ADM1 consolidated state-of-the-art modelling and provided a reliable, if imperfect, representation of the process, others drew on these earlier ideas to address specific challenges such as the effect of reactor start-up conditions on process stability [83] and reactor control under process disturbances [84] or shock loading events [85,86].

A major benefit of with representative and qualitatively reliable reduced models, which are tractable to mathematical interrogation, is their availability for deeper study of the non-linear dynamics and the implications on system performance. In effect, they represent in silico tools that can complement time-consuming experimental studies, allowing for rapid investigation of local and global system behaviour to aid process or experimental design, derive new knowledge and understanding, and even suggest new operating protocols. It must be reiterated that models cannot completely replace practical experiments, nor prove comprehensively the underlying process mechanisms in all but the simplest of systems [87].

Anaerobic digestion models have been extensively studied by mathematical biologists due to their generally reducible structure, the rich dynamics emerging from relatively few biotic interactions, and the modularity of components facilitating an array of hypothesis driven studies. Mathematical analysis provides an ability to study the asymptotic behaviour of dynamical systems, such as AD, i.e., we are able to observe whether the process tends to some steady-state (stable equilibrium) or undergo unforced cycling (functional periodicity), for a given set of parameters and initial conditions of the process [88].

Several mathematical analyses of existing AD models appeared in the 2000s. The Bernard model (see Section 2), describes a two-step (acidogenesis-methanogenesis) process and is tractable to dynamical systems analysis of its steady-state behaviour and stability. One of the principal benefits of simpler models are their availability for model-based process control, where a perfect representation of the typically non-linear system is not required [89]. An initial bifurcation analysis of the model at equilibrium allows for identification of the key process variables that can significantly impact process behaviour and, thus, determine the main candidates for control variables and their critical values. Rincón et al. employed the model to design an adaptive controller to avoid biomass washout in an anaerobic digester through identification of the normal form bifurcation (i.e., local behaviour around the system equilibrium) with dilution rate as the bifurcating parameter [90]. An analytical determination of the model equilibria and their stability using general growth functions (i.e., only the functional response and not the mathematical structure of biomass growth is assumed), showed the existence of bistability dependent on initial conditions, where the process may converge to the nominal (desired) operating point or undergo acidification and complete washout of the biomass [74]. Recently, the Bernard Model was again used, and a rigorous analytical approach employed for open-loop control of Biochemical Oxygen Demand (BOD) within a target design interval [91]. The model has also been used to investigate the correct determination of overloading tolerance, defined as the distance between a stable interior equilibrium and the limit point of the system, for correct process monitoring [75,92]. To address issues of process representation with this reduced-order model, especially at low microbial population sizes, a stochastic process has been suggested as a way to account for its structural and parametric uncertainty [93], demonstrating a promising, if not fully realised, potential.

Researchers have studied the dynamical behaviour of EBS more generally [94], analysing the role of ecological interactions (e.g., competition, mutualism), rather than the specific process, on the existence, uniqueness and stability of the system [95,96,97]. These methods have also been applied to specific processes, such as the ammonia removal (SHARON) process [98], Anammox systems [99], anaerobic chlorophenol mineralisation [100,101], and conventional anaerobic digestion [102,103].

Although global stability analysis is challenging for models of increasing order, an initial numerical analysis of the full ADM1 (44 state variables, 29 ODEs, 15 DAEs) was achieved using a prediction-correction continuation algorithm for solving the differential equations and the Newton-Raphson method for simultaneous solution of the algebraic equations [104]. They showed, through performing a one-parameter bifurcation that a limit point exists when either the dilution rate increases, or substrate inlet concentration increases above or decreases below critical threshold values. A positive feedback loop between acetic acid accumulation and acetate degrader washout due to inhibition was identified as the critical mechanism in all three cases. Interestingly, a reduced model that allowed for the derivation of analytical solutions matched the qualitative behaviour of the full ADM1, as long as key features (i.e., ammonia inhibition) are preserved in the simpler model). Further study of the Börnhoft model have shown the role of inhibition on multiplicity of equilibria [105], expanded knowledge of the richness of AD dynamics, and counter-intuitive result suggesting biogas production may not always be optimal at the co-existence equilibrium [106]. A recent study not only provides new theoretical results for the model but examines the effects of stochasticity to account for both environmental and parameter uncertainty (as with [93] and the Bernard Model). Interestingly, the author’s find that the process is more prone to failure at start-up and becomes resilient to stochastic effects as the microbial communities establish. This is perhaps an intuitive result given real-world examples of start-up failure in biological reactors, for example due to acclimation issues (e.g., [107]).

A relatively simple example of bifurcation analysis, for illustrative purposes, is provided in Appendix A, using the model first developed by Andrews (See Section 2) [29].

Despite a significant body of work demonstrating the power of mathematical analysis for enhancing AD process understanding and enabling practical controller design, there still remains a substantial barrier to the wider use of these techniques for engineering purposes. In part this is due to a dissonance between the level of mathematical acuity required to conduct the analysis and the substantial process knowledge in the domain of engineers and scientists. Although there do exist groups who have core competency in both fields, more formal engagement through communication, collaboration and bottom-up empirical studies to support theoretical observations, will be necessary if it is to gain traction in the wider EBS community.

4. Empiricism, Data and the Digital Future

It may be argued that the mechanics of anaerobic digestion is well-understood, and the maturity of phenomenological models facilitate their use for process design, operation and control. In data poor environments, where measurements are infrequent, erroneous or hard to obtain, mechanistic models are an obvious choice. Nevertheless, improvement in sensor technology and fast and cheap tools to characterise the microbial communities in digesters through multi-omic analyses, has gifted scientists and engineers with a wealth of data and knowledge of their ecology [108,109,110]. When used effectively, this can provide a greater resolution of understanding, helping to derive models that link microbial community composition and function effectively [111,112]. For example, metagenomic data as microbial biomarkers have recently been used in combination with a data-driven model for prediction of ammonia and phenol inhibition [113].

Unlike data-rich industries that employ extensive process monitoring, there is a sense that data-driven models for systems that are open to high levels of uncertainty and transient or irregular variation in their operation is under-exploited. Although process monitoring and automatic control have become standard, for example in largest wastewater treatment systems, there is still a tendency to rely heavily on existing knowledge and empirically driven decision-making, where the plant operators are the experts, while much of the data gathered at the plant is unused. This has come to be known as ‘dark data’, which is data collected by an industry that is either archived or discarded in lieu of capacity or motivation for its use at that time [114]. It can be argued that for AD plants, dark data is largely irrelevant given the quantity of data typically acquired, but as industries are increasingly integrated, relying on automation and knowledge extracted from data, so the opportunities offered by digital technologies, which can manage large and diverse amounts of information, become apparent [115].

Machine learning and intelligent systems (i.e., Artificial intelligence, AI) have become, seemingly, pervasive in industry and society in general. However, there are risks associated with their wholesale adoption, and not because they threaten to break Asimov’s three laws of robotics, as popular media may like to imagine. The main issues for concern, especially in motivating uptake by the water and AD industries, are (i) the potential for application redundancy by using ‘smart’ tools that provide little or no benefit where conventional engineering practice is more than adequate (human vs. machine expertise), (ii) hastily appropriated technologies without undertaking due diligence on model relevance/adequacy or consideration of the long-term system behaviour (garbage in, garbage out/inherent non-linearity and uncertain dynamics), (iii) advanced models and algorithms may require significant cost and expertise to implement and maintain. Nevertheless, it is apparent that developments in data acquisition and computing capacity facilitates the cautious use of ‘smart models’ for improvement of EBS. Indeed, many of the tools used now are based on fundamental theory developed over a century of statistical and mathematical research. For example, popular dimension reduction and regression techniques such as Principal Component Analysis (PCA) and Partial Least Squares (PLS) have been re-purposed for chemometric process monitoring [116,117] and control [118] of anaerobic digesters. These multivariate methods are used to transform high-dimensional data by matrix factorisation into a lower-dimensional space that retains the maximal information (variance) present in the original dataset (X for PCA) or maximal covariance between predictor and response datasets (X and Y for PLS). The relationship between the original variables and samples projected into this reduced-dimensional space represents a model of the system or process described by the original data

$$\begin{array}{cc}\hfill \mathbf{X}& ={\mathbf{TP}}^T+\mathbf{E}:[\mathrm{PCA}\&\mathrm{PLS}],\hfill \end{array}$$

(8a)

$$\begin{array}{cc}\hfill \mathbf{Y}& ={\mathbf{UQ}}^T+\mathbf{F}:\left[\mathrm{PLS}\right],\hfill \end{array}$$

(8b)

where T and U are the projections (or scores matrices) of X and Y, respectively, P and Q are the weights (or loadings matrices). E and F are the residual matrices. If the purpose of these models is for process monitoring and fault or disturbance detection, it is vital that the data used to construct them, and the number of principal components (factors in PLS) to retain (${\mathbf{T}}_l,{\mathbf{U}}_l$, where $l<m$, where m is the original number of variables in the X-block), represent its normal operating behaviour. New data acquired from the process is projected into the model and statistical measures (e.g., Hotelling’s T${}^2$ statistic) are used to identify samples that exceed some pre-determined threshold. The presumption is that these points represents non-normal process behaviour where variable contributions to the loadings may be examined to determine the source(s) of this exception.

PCA has been used as a data pre-processing step for extracting the correlation structure between variables. This variable selection method has been used as a preliminary step in the use of machine learning models to predict biogas yield from municipal AD [119]. A PCA model was used to detect disturbances in a sewage sludge fed AD process using PLS regression to chemical composition of the digestate from Near Infrared spectral data [120].

However, conventional multivariate models are limited by their stationarity [121], i.e., the model is a representation of the process over a fixed domain (e.g., time). To be useful for monitoring of dynamic processes, such as wastewater treatment and AD, adaptive/recursive methods have been developed for real-time monitoring and multivariate statistical process control, whereby the data model is updated with new information as it is acquired [122,123,124].

Unlike PCA, regression models, such as PLS, can be used to relate empirically derived data (measured explanatory variables) such as spectral data to one or more dependent (predicted) variables such as chemical concentrations, through their covariance structure. These regression models become useful in process operation when key performance indicators or control variables are difficult to measure. A model may be developed based on available and reliable measurements that are correlated with the target variable, acting as a soft-sensor for process monitoring and control. Like PCA methods, issues such as stationarity, multi-collinearity, and bias persist with ordinary PLS (OPLS). Robust method adaptations and rigorous data handling protocols are generally available (and required) to handle issues that may be encountered when using multivariate tools in practice [125,126]. Multivariate soft-sensors have been used quite widely by the wastewater community but fewer examples relating to the anaerobic digestion process have appeared [117,127,128]. This may be reflective of the general preference for using mechanistic models in AD or, anecdotally, the suspicion by practitioners that data-driven models are intangible and unrepresentative of the process. This may be exacerbated by the varied quality and quantity of data historically recorded at the plants themselves.

Nevertheless, ‘smart’ AD systems as a sub-class of a wider smartening of industry through the so-called 4th Industrial Revolution is an interesting and potentially valuable direction that practitioners may explore. This is especially pertinent in relation to the larger ‘Smart Water’ initiatives and repurposing of AD technology in bio-based manufacturing and biorefineries.

A renewed interest in Artificial Neural Networks (ANN) and Fuzzy logic, methods that have a long association with the bioprocessing and manufacturing industries, as well as AD itself [129,130], may propagate research and debate on the role of advanced data-driven models for application to AD moving forward. However, at present most applications of machine learning models have focused on conventional purposes, such as methane prediction [131,132,133].

In the final section we discuss how modelling is evolving by incorporation of interdisciplinary ideas and methods and show that no single approach is a catch all for engineered biological systems and anaerobic digestion, specifically.

5. Multi-Disciplinarity and the Future of Anaerobic Digestion Modelling

Mathematical modelling will play an increasingly important role in the design, monitoring and optimisation of EBS such as anaerobic digestion. As new challenges and opportunities emerge, technology develops and understanding increases, models will adapt and integrate the new knowledge acquired in the pursuit of smarter, cleaner and transformative processes. A recent review has highlighted possible trajectories for mathematical modelling of AD systems, both vertically (i.e., scale of models and connections across these scales) and horizontally (i.e., integration of processes, new anaerobic technologies) [22].

5.1. Hybrid Models

In the previous sections we have seen how disparate modelling approaches have been employed in the design, monitoring and control of AD. The selection of the appropriate model to use is dependent on its purpose, the quality, quantity and type of available process knowledge and data, the cost of model development and implementation versus the benefit of employing a model, among other considerations. It is a truism to state that no model is perfect, and often no single model is uniquely suited to the task in hand. For example, when there is lack of confidence in process measurements or data for characterising a process is unobtainable but a good mechanistic understanding has accrued, then a first-principles model seems logical. On the other hand, for many complex biological systems, only surface understanding of the system is known, and no consensus exists on the precise structural properties and interaction mechanisms, then a data-driven approach may be suitable. There is still debate, intensified by the current digitalisation epoch, as to whether machine learning is ready to entirely replace mechanistic modelling. As previously alluded to, the rapid adoption of AI and machine learning tools across many industries brings to mind the gold rush era. Some may succeed and prosper, while others may inevitably discover only fool’s gold. This is exemplified in the medical field, where early adoption of machine learning by clinicians is superseded by their acquisition in tackling more complex problems found in cell biology research. To quote a recent polemical opinion on this subject [134]:

Fundamental biology should not choose between small-scale mechanistic understanding and large-scale prediction. It should embrace the complementary strengths of mechanistic modelling and machine learning approaches …

Hybrid modelling is not a new concept, in fact many applications of what would be considered ‘pure modelling’ are a combination of structured and unstructured methods, e.g., consider parameter estimation techniques, such as the Extended Kalman Filter for fitting mechanistic models. However, explicit use of distinct modelling approaches has rarely been considered for application to bioprocesses. For anaerobic digestion, it appears viable that combining mechanistic and data-driven methods, together with extant process knowledge and empirical data would provide a means to overcome some of the shortcomings of either modelling approach alone [132].

An early example of hybrid modelling for AD employed a neural network to model process kinetics in combination with the Bernard Model described in Section 3 to ensure the ANN is constrained by the law of mass conservation (see Equation (5)) [135]. A similar approach using a kinetic model for lignocellulose fermentation and a PLS soft-sensor showed reasonable prediction potential over long time horizons [136]. The study is notable in its claim that it represents a simple demonstration of the digital twin concept. Although we do not discuss digital twins in this review, the concept that physical systems from individual process units to whole cities may be replicated digitally is becoming fairly ubiquitous, in tandem with concepts such as Industry 4.0 mentioned previously. Although it could be argued the use of buzzwords may attract attention to otherwise conventional research, all models may be classed as digital representations of real-world phenomena and some of these will feature in the establishment of digital twins in the future.

A neural network was used to predict AD performance and a Particle Swarm optimisation (PSO) algorithm is used to maximise methane yield, essentially employing a (meta)heuristic approach to evaluate a fitness function until a stopping criterion is reached (based on temperature and Volatile Suspended Solids (VSS), in this case) [137]. Nature-based methods, such as PSO, are not covered in this review, but more information on their emerging relevance to bioprocess monitoring, can be found in [138]. PCA was used as a feature extraction method to greatly reduce the dimensionality of ADM1 with relatively low loss of information, producing an identifiable model for supervisory control and optimisation of a winery effluent treatment plant [139]. A novel hybrid approach, augmenting ADM1 with information at the cellular level including intra-cellular fluxes and microbial activity, has recently been described. Flux Balance Analysis (FBA), a method to predict metabolic function from genomic information, was used to replace the Monod growth function in the ADM1 with an FBA growth prediction based on intra-cellular flux of metabolites [140]. The method allows for the dynamic prediction of optimal microbial composition, where conventional FBA provides steady-state analysis.

5.2. Thermodynamic Modelling

There has been considerable debate within the scientific community on the aptitude of purely kinetic models to describe the behaviour of systems that are subject to energetic constraints [28], especially in systems where microbial growth is directly linked to free energy availability. In anaerobic digestion, omission of energetic considerations in the mathematical description of the process dynamics will lead to qualitative and quantitative errors in model predictions, significantly undermining their usefulness. Indeed, the interplay between microbial kinetics and energetic constraints becomes important when the system approaches thermodynamic equilibrium, i.e., moving from exergonic sufficiency (net release of free energy) towards an endergonic state (net absorption of free energy) where biochemical reactions are completely inhibited (See Figure 3 for an example of the impact of hydrogen on thermodynamic inhibition of anaerobic propionate degradation). The rich microbial ecology that emerges in low energy environments, where thermodynamic inhibition is likely to constrain metabolic conversion, is both of interest and importance in the preservation of viable communities in AD [141,142] and the study of thermodynamic principles in biological systems has been an active area of research for decades [143]. Indeed, McCarty’s formative description of microbial growth energetics as the relationship of free energy of catabolism and the energy necessary for cell synthesis, the cell growth yield [144,145], is perhaps as important today as Andrews’ early conceptualisation of kinetic models for AD has been in the development of model formalisms for process control and optimisation over the last 20 years.

Various strategies have been proposed for the incorporation of thermodynamics into existing kinetic models or for the development of a universal theory. The Monod growth function is purely an empirical construct. Despite showing that the equation can relate with thermodynamic principles through the substrate affinity constant, i.e., $K_S={\mathrm{e}}^{-\Delta G^0/RT}$ [146], it is done so by making the assumption that only the limiting substrate and not microbial concentrations impact the growth dynamics. The Gibbs free energy change as the thermodynamic driving force for biochemical reactions forward was first coupled with a Michaelis-Menten type irreversible kinetic model describing anaerobic propionate degradation by Hoh and Cord-Ruwisch [147]. They surmised that reactions occurring close to thermodynamic equilibrium are reversible, driven by bidirectional reaction rates for the enzymatic conversion of substrates ($\left[S\right]$) to products ($\left[P\right]$), and the reverse [25,148], given by the Briggs-Haldane mechanism

$$\begin{array}{l}\mathrm{E}+\mathrm{S}\underset{{\mathrm{k}}_{-1}}{\overset{{\mathrm{k}}_1}{\rightleftharpoons }}{\mathrm{E}}_{\chi }\underset{{\mathrm{k}}_{-2}}{\overset{{\mathrm{k}}_2}{\rightleftharpoons }}\mathrm{E}+\mathrm{P}\end{array}$$

where E is the enzyme concentration, E$\chi$ is the enzyme-substrate-product complex, and ${\mathrm{k}}_1$, ${\mathrm{k}}_2$, ${\mathrm{k}}_{-1}$ and ${\mathrm{k}}_{-2}$ are the forward (${\mathrm{k}}_{+}$) and reverse (${\mathrm{k}}_{-}$) rate constants, which are combined to give the overall equilibrium constant ${\mathrm{K}}_{\mathrm{eq}}={\mathrm{k}}_1{\mathrm{k}}_2/{\mathrm{k}}_{-1}{\mathrm{k}}_{-2}$. Under substrate saturation conditions, the Haldane kinetic scheme (shown in Equations (2) and (3), within square brackets) dominates any thermodynamic constraints. However, as the substrate is consumed, a trade-off between kinetics and thermodynamics becomes apparent [141], and a rate equation for the reversible reaction can be derived

$$\begin{array}{c}\hfill \mu =\frac{\widehat{\mu }S\left(1-\Gamma /{\mathrm{K}}_{\mathrm{eq}}\right)}{K_S+S\left(1+k_r\Gamma /{\mathrm{K}}_{\mathrm{eq}}\right)},\end{array}$$

(9)

where $k_r={\mathrm{k}}_2/{\mathrm{k}}_{-1}$, $\Gamma$ is the mass-action ratio ($\left[P\right]/\left[S\right]$) at dynamic equilibrium). Given that the overall free energy change of reaction ${\Delta }_{r}G=RTln\Gamma /{\mathrm{K}}_{\mathrm{eq}}$ [149], then the general model proposed becomes [147]

$$\begin{array}{c}\hfill \mu =\frac{\widehat{\mu }S\left(1-{\mathrm{e}}^{{\Delta }_{r}G/RT}\right)}{K_S+S\left(1+k_r{\mathrm{e}}^{{\Delta }_{r}G/RT}\right)},\end{array}$$

(10)

In modelling the growth rate using the theory of reversible enzymatics, the energy required for driving the reaction and invested in maintaining other cellular functions not related to growth is overlooked, i.e., energy used for microbial growth is not the maximum free energy available in the system [150,151,152]. There has been significant developments in the understanding of thermodynamic principles for general application to microbial life, notably refinements in the derivation of metabolic energy requirements [153,154] and attempts to develop general or universal models to address shortcomings in established models, e.g., limits on metabolic parameter identification [155], and even the use of Maxwell-Boltzmann statistical physics to derive a fundamental link between microbial growth and energy states at the molecular level of substrates (i.e., the Microbial Transition State) [156,157].

However, despite a considerable amount of research on the application of thermodynamics to microbial systems, there are relatively few examples of direct application to the anaerobic digestion process itself. Oh and Martin used principles of thermodynamic equilibrium with chemical activity and fugacity together with a stoichiometric model to study energetic requirements of methanogensis [158], and recently Delattre et al. have proposed a kinetic-thermodynamic model to characterise a synthetic anaerobic tri-culture community performing sulphate reduction and methanogenesis [159]. Despite this, it is certain that as more consistent approaches are developed and theory more widely understood, appropriate thermodynamic models will be developed and tested for application to AD systems with the aim to improve performance and understand energetic limitations of different metabolic pathways. A recent paper has demonstrated that potentially, analytical methods from mathematics (see Section 3) can be used, despite the increased complexity, to elucidate the role of thermodynamics in biological growth models [160].

The thermodynamics of anaerobic digestion is inherently linked to the function of microbial communities and its inclusion in mathematical models has been shown to better predict inhibition of microbial growth under non-saturated substrate concentrations. Furthermore, there is an emerging trend for end-product diversification in AD systems, such as sustainable biorefineries producing high-value bio-based products [161,162,163]. The control and optimisation of these ‘manufacturing’ facilities will require models that represent metabolic pathway diversity and selection, especially the reversibility of novel catabolic pathways, through a general and robust framework describing the thermodynamic-kinetic trade-offs [142,152]. We anticipate significant integration of thermodynamic approaches with conventional engineering practice and microbial ecology in the coming years.

5.3. Computational Fluid Dynamics

Although ADM1 may be regarded as the most complete phenomenological description of the AD process from a bioprocess perspective, where the conditions of microbial growth are rigorously defined, there is an underlying assumption, as in many chemostat-type models, of spatial homogeneity, i.e., a completely mixed environment, or zero spatial dimensionality. This makes calculation of mass retention time and transfer rates relatively straightforward. However, in practice these ideal conditions are not observed and, at full-scale, the heterogeneous distribution of metabolites and biomass can lead to sub-optimal performance, dead-zones and unwanted reactor behaviour (e.g., foaming, physico-chemical gradients) [164,165].

Modelling the mass transport and hydrodynamic behaviour of fluids and particles in anaerobic digestion has been an emerging research field in recent years. Computational Fluid Dynamics (CFD) is an established modelling tool that uses the mathematics of fluid mechanics (e.g., Navier–Stokes, Euler equations) and numerical analysis to solve the governing equations on a 1- to 3-dimensional discretised grid [166]. As computational capabilities have increased (see Section 1), so the possibilities for using CFD to influence design and performance of AD reactors, especially related to mixing efficiency and identification of flow issues, have become a viable tool for multi-disciplinary research teams [167] and modelling consultancies (see Figure 4 for an example of CFD employed at full-scale). This is reflected by the increasing number of studies and software platforms that use CFD for AD research, including the performance of gas-mixing [168], physical mixing [169,170], process design and optimisation [171].

The use of CFD in anaerobic digestion processes, whether for reactor mixing and design or process optimisation, is likely to become a valuable tool for process engineers. At present it is used as a complementary approach with models that describe biochemical reactions through dynamic integration. However, recent studies have looked at using CFD models with biological kinetics and advection-diffusion transport of soluble metabolites and biomass, in a similar way to the methods for modelling spatially realised biofilm growth [172,173].

5.4. Individual-Based Modelling

We have seen that modelling of anaerobic digestion, and bioprocesses more broadly, have evolved and diversified as knowledge has accrued at different scales and computational capacity has dramatically increased. Although mechanistic models are likely to persist as the method of choice for engineers and practitioners, other approaches will emerge that offer those working in AD with tools to enhance their understanding and management of these systems. One such method that is highly granular and explicitly resolved is individual-based modelling. Individual-based Models (IbMs) and associated tools, such as agent-based models and cellular automata (CA), are spatio-dynamical and discrete in time, space and states. Although CA are typically used to predict geometric distributions based on local interactions, IbMs are highly resolved at the cell or individual level [79,174]. Unlike aggregated and continuous population models, IbMs proscribe behavioural physical and biological traits to individuals in a highly discretised manner, where the complex, global behaviour of the population emerges from the simpler, local properties and interactions among individuals.

IbMs have historically been used in macro-ecology but found significant purchase as a tool for modelling biofilms, building on knowledge gained from spatio-temporal, multi-species biomass-based models, such as those briefly mentioned in Section 3 [175,176]. In addition to the microbial growth properties, physical interactions between individuals, motility, environmental gradients, fluid dynamics, thermodynamics, and even stochasticity may also be included as hybrid components to a more comprehensive process, sub-process or reactor scale model [177,178,179]. Such a multi-scale model is complex and requires a large degree of multi-disciplinarity to ensure successful implementation. However, the use of IbMs is becoming more ubiquitous and recent efforts to develop large-scale simulators of EBS using hybrid IbM and statistical emulators have established a framework for studying microbial community behaviour at scale [71,180,181], opening the door for application to anaerobic digestion processes [179,182] (See Figure 5 for an example application of an IbM applied to anaerobic processes).

6. Conclusions

With five decades of anaerobic digestion mathematical modelling there may be an assumption of consensus among the scientific and engineering communities who have directly or indirectly aided the development and validation of models from theory to practice. Although it could be argued that there is a level of maturity embodied by the most well-known and used tool, the Anaerobic Digestion Model No. 1, we have seen that such highly structured, complex models do not fit all purposes. Indeed, many applications of empirically derived mechanistic models rely on simplifying assumptions to obtain a structure that is identifiable, without loss of specificity or accuracy, e.g., for process control. Here we argue that the rigorous application of mathematical analysis tools, typically in the domain of theorists, can be appropriately applied to reduced-order models of AD to investigate the qualitative behaviour of the system, identifying unexpected behaviour, guiding experimental studies, and suggesting suitable control parameter limits for improved process operation.

Many of the bottlenecks that have limited not only modelling but more general AD research, such as reliable and in situ instrumentation, limited understanding of the complex microbial ecology of digesters, and restrictive computational capacity, have generally been overcome in the last two decades. This has manifested, more widely, in the large availability of data in the bioprocess industry, opening up opportunities for AD modelling across scales, from the cellular and metabolic level, through to plant-wide integration and supervisory control. Although data-driven models and digital tools based on machine learning algorithms and AI are emerging in parallel with the Industry 4.0 ethos of “Smart Systems”, it is cautioned by the well-known aphorism “garbage in, garbage out”. There are a significant number of publications on the use of AI techniques such as Artificial Neural Networks, but these may be subject to the problem of over-fitting and loss of generality. Investment in time and energy to develop complexity is driven by an assumption that the more information the better the model. This is not strictly true, many soft-sensors employed today are based on simple multivariate relationships between a few measured variables, where a standard PLS regression model may be adequate.

However, to address more challenging issues faced by industry, it is clear that some standard modelling approaches may not be fit for purpose. The future of AD modelling will track the needs of the industry but also incorporate ideas from other fields to ensure that models are robust, reliable, flexible and adaptable to local and global shifts in the process domain. We are already seeing this happen in the water sector, where environmental challenges such as climate change, increased urbanisation and economic constraints, are driving industries to look at the potential for models to ensure that infrastructure is resilient in the face of novel demands. AD technology is a key component of the energy-food-water nexus and drivers already exist for collaboration between modelers, engineers, scientists and other stakeholders to exploit the potential of available models but also to integrate emerging and innovative methods into their thinking around AD, such as hybrid modelling techniques, cellular level and metabolic knowledge, thermodynamic constraints on microbial growth, fluid dynamic process visualisation, and fine-grain individual scale models that can incorporate all of these elements.

This review does not suggest a preference for any model or modelling approach. There is no ideal model for AD and often several models may achieve the same goal. Simply, we present a broad overview of the historical context of AD modelling and look to the future as the wealth of knowledge and information grows. The use of models for AD is strengthening and we envisage that this will become even more important in the context of the climate emergency, the need for smart management of resources and in the provision of resilient infrastructure that generate and process these resources.