Synergising biomass growth kinetics and transport mechanisms to simulate light/dark cycle effects on photo‐production systems

Light attenuation is a primary challenge limiting the upscaling of photobioreactors for sustainable bio‐production. One key to this challenge, is to model and optimise the light/dark cycles so that cells within the dark region can be frequently transferred to the light region for photosynthesis. Therefore, this study proposes the first mechanistic model to integrate the light/dark cycle effects into biomass growth kinetics. This model was initially constructed through theoretical derivation based on the intracellular reaction kinetics, and was subsequently modified by embedding a new parameter, effective light coefficient, to account for the effects of culture mixing. To generate in silico process data, a new multiscale reactive transport modelling strategy was developed to couple fluid dynamics with biomass growth kinetics and light transmission. By comparing against previous experimental and computational studies, the multiscale model shows to be of high accuracy. Based on its simulation result, an original correlation was proposed to link effective light coefficient with photobioreactor gas inflow rate; this has not been done before. The impact of this study is that by using the proposed mechanistic model and correlation, we can easily control and optimise photobioreactor gas inflow rates to alleviate light attenuation and maintain a high biomass growth rate.

and/or fluid dynamics have been proposed to assess the performance of PBRs.
On the one hand, macroscopic models that combine the Monodtype empirical biokinetic equations with light transmission functions (Socher et al., 2016) have been previously used to optimise the geometry of pilot-scale PBRs (Ali et al., 2019) and analyse biomass productivity in different sizes of open raceway ponds (Amini et al., 2018;Park & Li, 2015). However, effects of light/dark cycles on biomass growth and reactor design have been barely included. Light/dark cycles are induced by fluid dynamics and light attenuation. A short light/dark cycle caused by an intense culture mixing can reduce cells' residence time in the photic zone near the irradiated PBR surface and the dark zone in the PBR interior, thus preventing photosynthetic proteins from over-reduction and overoxidation for algae and cyanobacteria species (Janssen et al., 2003;Luo & Al-Dahhan, 2004). This improves the light utilisation efficiency and microbial photosynthetic activity (Han et al., 2000;Vejrazka et al., 2012). Although there have been a few studies simulating this phenomenon (Li et al., 2014;Pruvost et al., 2008), they did not integrate light/dark cycles into the biokinetic models. Instead, external computation of cell light exposure history was performed separately and then averaged to solve the cell growth kinetic model. This simplification may introduce large errors when applied to large scale PBRs.
On the other hand, microscopic kinetic models derived based on photosynthetic units (PSUs) have also been embedded into computational fluid dynamics (CFD) along with light transport models to simulate spatio-temporal changes of intracellular PSU states and culture biomass growth Gao et al., 2018;Papáček et al., 2014). Although these models have achieved satisfactory results when compared with experimental data, their underlying coupling strategy was computationally expensive and lacks theoretical support. This is because the PSU models are primarily developed to simulate intracellular photosynthetic reactions without the inclusion of macroscale biomass concentration changes (Han, 2002;Wu & Merchuk, 2001). Thus, when embedded into the CFD model, this multiscale reactive transport model is only valid over a small time period when biomass concentration changes are negligible. As a result, the model can only be run within a small time step to estimate the spatial ratio between different PSU states under the current light intensity distribution. This ratio is then used to estimate a volume averaged growth rate to calculate biomass concentration. The updated concentration is finally used to calculate local intensity distribution at the next time step. Hence, the lack of direct coupling between intracellular kinetics and fluid dynamics causes a high computational cost which severely limits the application of this approach for PBR design.
Overall, at the meso/macroscale, most biokinetic models assume that fluid dynamics and light transmission are not correlated.
Despite the fact that fluid dynamics directly affect the frequency of light/dark cycles, there is lack of accurate models integrating fluid dynamics into cell growth kinetics. Furthermore, theoretical connections between macroscale empirical biomass growth kinetics (e.g., Monod-type models) and intracellular photosynthetic reaction kinetics (e.g., PSU models) are still missing. It is unclear how the parameters in an empirical model are linked to the metabolic reactions. Thus, it is challenging to modify the model structure when trying to incorporate more physical/biological information.
Therefore, to address these challenges, the present study proposes a new physical model linking the synergistic effects of light/ dark cycles and fluid dynamics into biomass growth kinetics. This model is proposed based on a theoretical derivation using intracellular PSU reaction mechanisms. The remainder of this paper is structured as follows: the biokinetic model derivation and CFD coupling strategy for PBR simulation are presented in Section 2. New correlations between gas inflow rate (indicative of fluid dynamics) and light related kinetic parameters along with their practical application are discussed in Section 3, followed by the final conclusion and suggestions for future work. PSUs are light saturated thereby using the absorbed energy to begin photo-production (photochemical quenching), and (iii) an inhibited state x 3 , where the PSUs are temporarily damaged under high light intensities but can be recovered to the closed state after cellular self-reparation. This state transitioning kinetics is described by Equations (1) to (4). (3) where α and β are reaction rate constants for light-dependent reactions from x 1 to x 2 and from x 2 to x 3 , respectively, and γ and σ are reaction rate constants for light-independent reactions from x 2 to x 1 and from x 3 to x 2 , respectively.
Equation (5) can be exactly presented as the commonly used Aiba model which accounts for photo-limitation, photosaturation and photo-inhibition (Aiba, 1982). In spite of the wide application of this model, k s and k i are often considered as empirical terms, and their respective biological meaning has never been clarified. This derivation, therefore, filled this gap. It is seen that these terms represent the ratios of different intracellular reaction rate Under such conditions, cells are exposed to light/dark cycles which limits biomass growth. The frequency of these light/dark cycles are largely influenced by culture mixing which is dictated by fluid dynamics. As illustrated in Figure 1b, at each position in a PBR with a specific local light intensityI l , cells can move back and forth due to culture mixing, thus being exposed to either a "lighter" region where they can absorb more light or a "darker" region where less solar energy resource is available. Thus, the actual light intensity that cells experienced during their growth could be significantly different from the local light intensity I l . To account for this effect, we defined a new coefficient, effective light coefficient ηo estimate the "effective" local light intensity (η•I l ) based on engineering principles. Therefore, Equations (1) to (3) are modified to Equations (6) to (8): Solving Equations (6) and (8) in the same manner results in a new biokinetic model taking into account light/dark cycle as Equation (9). The effective light coefficient depends on the culture mixing which is utterly determined by the mixing energy if the PBR configuration is fixed. In practice, as mixing energy is difficult to be measure directly but is closely linked to the gas inflow rate for many PBRs (if gas inflow rate is used for culture mixing), it is more convenient to link gas inflow rate with the effective light coefficient. Hence, by changing the gas inflow rate, we can control biomass growth rate during biomass cultivation. The correlation between η and gas inflow rate will be determined in Section 3.4. In principle, if culture mixing in a PBR is carried out using other approaches, then gas inflow rate can be replaced by other parameters such as impeller rotation rate: To estimate η and investigate how it is correlated with the PBR fluid dynamics, two in silico case studies are presented in the following sections.

| Introduction to the in silico PBRs
The respectively, so that total volumetric gas flowrate is same in the two scenarios. The gas inflow rates selected here correspond to a typical range of superficial gas velocities used in PBRs (Huang et al., 2016;Yang et al., 2014;Yu et al., 2009). The PBRs were unidirectionally illuminated with fluorescent lamp source (Wu & Merchuk, 2001) at an incident light intensity of 300 μmol m −2 s −1 .
This is often used for microalgae experimental cultivations (Li et al., 2014;Rio-chanona et al., 2017). Nutrient concentration (i.e., CO 2 ) was assumed to be sufficient and therefore their effects are not included in the model.

| CFD simulation
Coupling CFD with biokinetic model for process data generation remains a challenge as discussed in the introduction. To solve this issue, a new approach was proposed here to save computational cost while guaranteeing a high simulation accuracy.

| Selection of the gas-liquid multiphase model
A two-phase fluid (water and air) model assuming that the microbial cells follow the water phase without interference is applied for the multiphase PBR system using the classical Eulerian-Eulerian approach (Papacek et al., 2018). This results to the continuity and momentum equations (10) and (11), respectively: where subscript = k L G , stands for liquid or gas phase, ρ u r , , k k k corresponds to the phase volume fraction, density, and velocity vector, respectively. The pressure and gravitational vector are denoted as g p, , respectively. μ k e , is effective viscosity which includes both molecular viscosity, μ k m , and turbulent viscosity μ k T , . μ k T , is modelled using the standard ε − k turbulence model (Pfleger et al., 1999). I is the unit tensor.
Drag effect being more dominant is the only interphase momentum force which is considered herein. It was modelled with the Schiller-Naumann drag model as it is widely used in bubble column PBRs (Nauha & Alopaeus, 2013). Detailed formulations of the standard ε − k turbulent model and the Schiller-Naumann drag model can be found in Rampure et al. (2003) and Luo and Al-Dahhan (2011).

| Rigorous biomass transfer model
To account for both light limitation and inhibition, the three-PSUstate kinetic model is used in this case. The mass transport model for biomass concentration accounting for both time and spatial dimensions is therefore expressed as Equation (12). The first term on the left stands for biomass accumulation, the second term on the left denotes biomass convection, the first term on the right denotes biomass diffusion, and the second term on the right represents biomass growth. Attention must be paid to the fact that Equation (12) does not need to include the effective light coefficient η, as this equation captures the instantaneous biomass growth at each time and each location inside a PBR. In addition, a decay term is also added in Equation (12)  of culture shear stress which is also dependent on fluid dynamics.
As the current study only investigates biomass growth, the decay rate constant is set fixed. This assumption is found to be true as long as the shear stress does not exceed a threshold, beneath which shear stress has negligible effect on cell decay (Leupold et al., 2013): where u L and ρ L are the culture velocity and density (using water density), respectively, = × − D 5.5 10 L 10 m 2 s −1 is biomass viscosity, Equation (13)

| Parameter estimation of the dynamic model for the effective light coefficient
Once data is generated from the multiscale model, it is used to estimate η in Equation (9). The effective light coefficient represents the effect of light/dark cycles on biomass growth kinetics, and this equation can be served as a mechanistic surrogate model to enable efficient process simulation and real-time control. To include light attenuation and cell decay, Equation (9) is modified as Equation (14), which can be approximated as Equation (15) to eliminate the spatial dimension effect (light transmission) (Zhang et al., 2015). In this way, overall biomass concentration inside a PBR as well as its correlation with light/dark cycles can be calculated: where L is the PBR width (total light transmission length).
However, Equation (15) is a highly nonlinear ordinary differential equation, imposing great obstacles when performing parameter estimation to calculate the value of η under different gas inflow rate using gradient-based optimisation algorithms. Therefore, a stochastic hybrid optimisation algorithm was used to estimate the value of η. This hybrid algorithm utilises the particle swarm optimisation (PSO) and the artificial bee colony (ABC) algorithms in parallel.
The hybrid algorithm is initialised with a discrete number of inputs, which are evaluated with the objective function (negative loglikelihood) and sorted into different groups. A specified percentage (50% in this study) of the best inputs are used by the PSO, with the rest used by the ABC (Karaboga & Akay, 2009). This process is repeated over a specified number of iterations (10 iterations in this study), and the best input found is returned. The advantage of this hybrid algorithm is that it uses the explorative characteristics of ABC to explore the search space of the worse inputs (thus finding potentially good solutions) and uses the exploitative characteristics of PSO (Kennedy et al., 2001) to exploit the search space of the better inputs. The computational time to attain the solution is 2.5 min in this study.

| Computational cost analysis
All the CFD simulations were executed on a workstation.  (Amini et al., 2018) and from 0.0005 to 10 s over 120 h of cultivation .

| In silico experimental data validation
To examine the current multiscale modelling strategy, several properties were analysed to compare with the literature observations. Figure 3a-c show typical instantaneous water displacement for Scenario 1 displayed on a cross sectional plane (x − y) at z = 0.0125 m, for better visualisation. It is seen that a dynamic vertical motion from the bottom (sparger) to the centre top PBR sections with randomly oscillating plume was observed similar to that of (Huang et al., 2015;Yang et al., 2014;Yu et al., 2009). This was interpreted to be other studies (Yang et al., 2014;Yu et al., 2009;Zhao et al., 2020).
To check for the flow regime induced by the rising bubbles within the FP-PBR, a Reynolds number Re is defined using Equation (16) which was reportedly used by Pfleger et al. (1999) for 3D multiphase CFD simulation of FP-PBR: In this study, the value of Re is 10,366 for the lowest gas inflow rate and 15192 for the highest gas inflow rate, thus indicating turbulent flow regimes and validating the selection of the multiphase turbulent flow model (Equation 11). This supports the chaotic velocity flow pattern observed inside the FP-PBR as visualized in Figure 3. The same analysis was also carried out to evaluate the data accuracy in Scenario 2, thus not repeated here.

| Effect of gas inflow rate on biomass growth
Only biomass growth for the first 144 h of cultivation is presented since there was no observable change in biomass concentration after this time (stationary phase). Figure 4a shows the biomass evolutiontime profiles of all test cases in Scenario 1. These profiles are observed to attain a saturation biomass concentration ranging from 4.2 to 5.0 g L −1 which is realistic for a batch photo-production process operated under good aeration conditions (Gao et al., 2018;Richmond, 1996;Zhao et al., 2020). The maximum biomass concentration increased by 16.7% from a gas inflow rate of 0.1 to 2.0 m s −1 . This was expected as increasing gas inflow rate leads to a better culture mixing and a more frequent light/dark cycle for cells to absorb solar energy, in alignment with previous experimental observations (Grobbelaar, 1994;Janssen et al., 2003;Richmond, 1996).
This can be further proved given that the liquid phase velocity along the light transmission direction increases by 307.7% with the increasing gas inflow rate from 0.1 to 2.0 m s −1 (see Figure 4b), shortening the light/dark cycles (Huang et al., 2016). The increase of biomass concentration under different cases also proves the accuracy of the multiscale model coupling strategy proposed in this study.
Most importantly, it is expected that maximum biomass concentration will only increase mildly when increasing the gas inflow rate (Zhao et al., 2020). This is because biomass growth kinetics does not directly depend on fluid dynamics, but culture mixing can affect the light exposure history that cells experienced during their growth.
As a result, an intensified culture mixing will not increase biomass concentration as dramatic as what observed when changing other key factors such as temperature, pH, or nutrient supply. Finally, Scenario 2 also yields a similar conclusion as Scenario 1, thus not presented here.

| Correlating the effective light coefficient with gas inflow rate
The estimated values of effective light coefficient, η for all test cases in Scenarios 1 and 2 are summarised in Table 1 while  light for their growth due to culture mixing, and that purely using a light attention equation to calculate a PBR's "apparent" dark region may overestimate the effect of light attenuation on biomass growth kinetics and PBR upscaling. This greater effective local light intensity also suggests a better light utilisation efficiency since culture mixing also prevents cells exposed within the light zone from being overoxidised due to the long-term absorption of excessive photons (e.g., photo-inhibition). Secondly, the effective light coefficient η increases with an increasing gas inflow rate in each scenarios: a 159.8% increase from 0.1 to 2.0 m s −1 in Scenario 1 and an 86.9% increase from 0.25 to 1.0 m s −1 in Scenario 2. This implies that η is positively correlated with an increased liquid phase velocity and a shortened light/dark cycle, both induced by the increasing gas inflow rate.
When comparing the two scenarios, it is seen that the correlation between η and gas inflow rate is also dependent on the PBR geometry. This is understandable as fluid dynamics are influenced by the reactor configuration.
In theory, it is more appropriate to correlate the effective light coefficient with Re. However, as it is difficult to directly calculate Re for a multiphase system, in practice it is more convenient to correlate the effective light coefficient with gas inflow rate which fundamentally governs the PBR fluid dynamics. Therefore, Equations (17) and (18)  From Equations (17) and (18) and Figure 6, three important observations are drawn. Firstly, the correlations between η and u follow a power-law expression in both scenarios. This is similar to most correlations used in the field of transport phenomena. For example, power laws are commonly used to correlate mass transfer coefficient with gas inflow rate for gas-liquid reactor simulation (Fujasová et al., 2007;Moucha et al., 2003), and the power index is determined by the momentum and mass transfer behaviours at a molecular level. Physically speaking, the effective light coefficient represents the significance of local light transmission (local photons mass transfer) and microbial cells mass transfer, thus the correlation between effective light coefficient and gas inflow rate should obey a power-law expression. Secondly, parameters in the two correlations are different as they are dependent on the PBR geometry and power input. This is also similar to the mass transfer coefficient in a gas-liquid reactor as the power index corresponding to gas inflow rate and power input can vary over twice if the impeller configuration is different (Fujasová et al., 2007 rate beyond which the correlation will be invalid. This is either attributed to the fact that an increasing gas inflow rate will cause a higher shear stress which in turn leads to cell decay (Leupold et al., 2013), or the fact that the culture mixing pattern has been changed. Therefore, future work should also pay attention to the valid range of this correlation.

| Controlling PBR operation
Increasing incident light intensity (offsetting light attenuation) and enhancing culture mixing (shortening light/dark cycles) are two strategies to facilitate biomass growth in a PBR (Schulze et al., 2020).
Although model-based optimal control of light intensity has been reported previously (Del Rio-Chanona et al., 2019;Koller et al., 2017), optimal control of gas inflow rate to enhance culture mixing has never been achieved as such a model was never proposed before.
Therefore, by using Equation (15) (15) and (13), respectively. In this way, cells can always maintain a high growth rate and the presence of dark region will not significantly restrict biomass growth.

| CONCLUSION
In this study, the biological meaning of k s and k i in the Aiba model is clarified. Based on this knowledge, an effective light coefficient was proposed and embedded into a macroscopic biokinetic model to account for the effect of light/dark cycles on biomass growth. To estimate its value, a novel multiscale modelling strategy was developed to greatly reduce the computational time cost, meanwhile guaranteeing high simulation accuracy. Most importantly, the physical insight, practical applicability, and current limitations of the proposed effective light coefficient were thoroughly discussed. By correlating the effective light coefficient with PBR gas inflow rate, it is now feasible to optimally control the gas inflow rate at each time to maximise cell growth and mitigate light attenuation for biomass cultivation. Considering the early stage of this study and simplifications (e.g., sufficient nutrient supply, constant culture physical properties, fixed PBR configurations) used for theoretical derivation, more thorough investigation will be conducted in the future work to generalise the applicability of current work for PBR design, upscaling and real-time control.