Experimental and Model-Based Analysis to Optimize Microalgal Biomass Productivity in a Pilot-Scale Tubular Photobioreactor

A dynamic coarse-grained model of microalgal growth considering light availability and temperature under discontinuous bioprocess operation was parameterized using experimental data from 15 batch cultivations of Nannochloropsis granulata in a pilot-scale tubular photobioreactor. The methodology applied consists of a consecutive two-step model parameter estimation using pooled, clustered and reorganized data to obtain initial estimates and multi-experiment fitting to obtain the final estimates, which are: maximum specific growth rate μmax = 1.56 d−1, specific photon half-saturation constant KS,ph = 1.89 molphgX-1d-1, specific photon maintenance coefficient mph = 0.346 molphgX-1d-1 and the cardinal temperatures Tmin = 2.3°C, Topt = 27.93°C and Tmax = 32.59°C. Biomass productivity prediction proved highly accurate, expressed by the mean absolute percent error MAPE = 7.2%. Model-based numerical optimization of biomass productivity for repeated discontinuous operation with respect to the process parameters cultivation cycle time, inoculation biomass concentration and temperature yielded productivity gains of up to 35%. This optimization points to best performance under continuous operation. The approach successfully applied here to small pilot-scale confirms an earlier one to lab-scale, indicating its transferability to larger scale tubular photobioreactors.


INTRODUCTION
Microalgae as a phylogenetically non-related group of organisms are adapted to a wide variety of inhabited ecosystems (Hoffmann, 2010;Guarnieri and Pienkos, 2014). The mostly photosynthetically active organisms produce a number of substances, such as carotenoids and fatty acids, that are specific to primary producers. A number of them, such as ketocarotenoids (e.g., astaxanthin) or long-chain ω-3 fatty acids (e.g., eicosapentaenoic acid, C20:5, EPA), are of commercial interest to the food and feed industry (Posten, 2012;Minhas et al., 2016).
In comparison to plants, microalgae are in general structurally less complex and show higher area productivities . Their gross biochemical composition can be influenced without genetic modifications (Gu et al., 2012;Draaisma et al., 2013). These specific characteristics cause growing interest in autotrophic biotechnological production processes based on microalgae. Their products are of particular interest to markets such as aquaculture, nutraceuticals and cosmetics (Khan et al., 2018;Malcata, 2018). This is accompanied by worldwide growing research efforts (Endres et al., 2018;Garrido-Cardenas et al., 2018;Lippi et al., 2018).
Tubular photobioreactors are an important type among closed photobioreactor systems used for the phototrophic production of microalgal biomass worldwide (Takache et al., 2009;Grewe and Griehl, 2013;Karemore et al., 2015;Olaizola and Grewe, 2019). Biomass productivity and hence space-time yield (often referred to as volumetric yield) of a photoautotrophic bioprocess are influenced by many abiotic factors, which also change during the course of the day as well as throughout the cultivation cycle (Bernard et al., 2016). Light availability can vary by orders of magnitude along with biomass concentration, layer thickness and external light intensity and is therefore difficult to control, e. g. in outdoor discontinuous bioprocesses. Temperature can also vary in a wide range of 2, . . . , 30 • C (Zittelli et al., 1999;Ras et al., 2013), in particular at larger scale, while pH is easier to control even in industrial applications.
The work presented here employs a dynamic coarsegrained model to describe the light and temperature dependent specific microalgal growth rate within a pilot-scale tubular photobioreactor in order to finally predict biomass growth and optimize biomass productivity. Within the model used, light-limited growth is described depending on the specific light availability rate q ph . A similar approach has already been successfully applied to optimize steady-state biomass productivity within a lab-scale tubular photobioreactor under turbidostatic operation (Weise et al., 2019). The present work extends this approach to dynamic microalgal growth within discontinuous processes. While under turbidostatic operation q ph can be controlled well, under discontinuous operation it can only be kept within reasonable ranges. This is investigated here for a repeated batch process, with the inoculation biomass concentration c X,0 and the cultivation cycle time t cyc as process parameters that can be influenced during operation.
The current work additionally includes temperature dependencies into the model. Therefore, the empirical Cardinal Temperature Model with Inflexion (CTMI), first introduced by Lobry et al. (1991), is employed. The CTMI was already successfully applied to describe microbial and microalgal temperature dependencies (Rosso et al., 1993;Bernard and Rémond, 2012;Barbera et al., 2019).
Therefore, the following growth kinetics model for the description of both dependencies, i. e. for light and temperature, was also used here (Equation 1) (Bernard and Rémond, 2012): ) represents the optimum specific growth rate for a specific light availability rate q ph [mol ph g X −1 d −1 ] at the strain-specific optimum temperature T opt [ • C], whilst φ(T) [-] (Equations 5 to 7) describes the influence of the temperature T [ • C] on the optimum specific growth rate.
and K S,ph [mol ph g X −1 d −1 ] are the maximum specific growth rate and the specific half-saturation constant for photon availability, respectively. The specific photon maintenance coefficient m ph [mol ph g X −1 d −1 ] represents the light energy used for purposes other than biomass synthesis (Pirt, 1986). This kinetics describes the experimentally observed saturation of the growth rate when light availability increases within the investigated range (Darvehei et al., 2018). Regarding Nannochloropsis spp., the literature reports maximum specific growth rates µ max in the range of 0.86, . . . , 1.6 d −1 (Sandnes et al., 2005;Spolaore et al., 2006). Specific photon availability rates q ph are often found in the range of 0.25, . . . , 2.2 mol ph g X −1 d −1 (under lit conditions; flat panel photobioreactors) (Zijffers et al., 2010;Kandilian et al., 2014;Janssen et al., 2018).
The biomass-specific photon availability rate q ph (Equation 4) is calculated from process and geometrical characteristics, in particular from the biomass concentration c X [g X m −3 ], the light intensity at the reactor surface I 0 [mol ph m −2 d −1 ], the illuminated reactor projection surface A [m 2 ] and the total reactor liquid volume V L [m 3 ]. The underlying relationships have been successfully applied and demonstrated in several publications (Bernard, 2011;Kliphuis et al., 2011;Blanken et al., 2016). Equation (4) applies only to light-limited growth conditions and was originally developed considering flat-panel photobioreactors (Schediwy et al., 2019). Therefore, the tubular reactor compartment is considered to be a flat-panel equivalent with an average light path length l ∅ (Equation 3). Since V L has to be equal in both assumptions, the illuminated reactor projection surface is used for q ph calculation in Equation (4).
Within Equation (4) and T opt > T min + T max 2 (7) Outside T min and T max [ • C], which represent the minimum and maximum temperature of the growth range, no growth is assumed (Equation 5). In order to obtain the inflected asymmetric shape of the CTMI, T opt values need to be closer to T max as to T min (Equation 7). Bernard and Rémond (2012) verified both the above model properties (Equations 5, 7), which are necessary for a successful model application in this context. Optimum growth temperatures are stated in the literature within a range of 20, . . . , 29 • C (Sukenik, 1991;Bartley et al., 2015;Abirami et al., 2017). However, to the best of our knowledge, no optimum growth temperature for N. granulata has been published yet. Model parameter estimation was carried out in two consecutive steps: Initial parameter estimates were generated in a first step from pooled, clustered and reorganized data of the process. These initial estimates were then used in a second step for the final parameter estimation based on the original experimental time series data.
Based on the parameterized model, a numerical optimization regarding the biomass productivity Pr was carried out with respect to the process parameters inoculation biomass concentration c X,0 , cultivation cycle time t cyc and temperature T. Also, suggestions are provided with respect to the transfer of the proposed approach to other tubular photobioreactors under repeated discontinuous operation.

Cultivation Conditions
The strain of Nannochloropsis granulata (Karlson et al., 1996) used in this study was previously isolated and its identity confirmed by 18s rDNA analysis carried out by SAG Göttingen. All experiments were performed using sterile brackish water medium with a salinity of 19 psu, a nitrate concentration of 12 mM and a phosphate concentration of 0.3 mM provided by Prof. Otto Pulz, IGV GmbH, Nuthetal, Germany. Both, nitrate and phosphate concentrations were measured using ready-to-use cuvette test kits (nitrate: WTW 252073; phosphate: WTW 252075, Xylem Analytics Germany, Weilheim, Germany) after centrifugation of samples at 4, 500 * g for 20 min. Cell-free supernatant was diluted with double distilled water 1:100 for nitrate and 1:10 for phosphate analysis according to the manufacturer's instructions. Absorbance values and ion concentrations were determined in a photoLab ® 6100 VIS Photometer (Xylem Analytics Germany, Weilheim, Germany). Nitrate and phosphate were added manually on demand from sterile stock solutions of NaNO 3 (300 g l −1 ) and KH 2 PO 4 (100 g l −1 ) if values fell below 50% of the initial medium concentrations (c NO − 3 = 1,021 mg l −1 , c PO −3 4 = 98 mg l −1 ). Maintenance cultures of N. granulata were kept in cell culture flasks at 14 • C and 50 µmol ph m −2 s −1 light intensity (12 h light/12 h dark).
Preculturing was carried out in 1.7 l double jacked glass bubble columns (height to diameter ratio of 4), temperated at 22 • C and aerated at 0.1 vvm with a 2% CO 2 to air ratio (v/v). Precultures where illuminated using fluorescent tubes (Philips Cool White, 36 W, eight tubes per column, horizontally arranged). The light intensities at the inner surfaces of the columns were on average 249 µmol ph m −2 s −1 , measured by a spherical micro quantum sensor (US-SQS, Walz GmbH & Co KG, Ulm, Germany). During preculture, illumination was provided to the columns continuously for 24 h.
The 30 l photobioreactor was inoculated using precultures from four individual bubble columns which were originally inoculated by the same maintenance culture.

Experimental Set-Up
The experimental work carried out within this study was conducted using a 30 l pilot-scale tubular photobioreactor (PBR30, IGV GmbH, Nuthetal, Germany) consisting of a compartmented glass tube (borosilicate glass, inner diameter of 39.6 mm) and a stainless steel system vessel (see Figure 1). Illumination was provided permanently over 24 h to the glass tube of the reactor by fluorescent tubes (OSRAM Cool White, L 18 W/840) and five LED panels placed along the length of the compartmented glass tube (Figure 1). The LED panels had a  Figure 1A), I 0 (y i z j ) -surface light intensity of the respective compartment, y i -ith position in y-direction, z j -jth position in z-direction.
Frontiers in Bioengineering and Biotechnology | www.frontiersin.org TABLE 1 | Upper part: Light intensity at the reactor surface of the compartments used within the modeling approach, ⋆ Measuring positions (compare Figure 1A); Center part: Reactor-specific and strain-specific constants; Lower part: Set parameter constraints.  Sorokin and Krauss (1958); b Zijffers et al. (2010).

Light intensity at the reactor compartments surface
spectral composition consisting of 95% red light (λ = 660 nm) and 5% blue light (λ = 440 nm). The illumination profile of the reactor's tubular glass part was recorded at 140 measuring positions (indicated by ⋆ in Figure 1 and Table 1) at 7 positions in x-direction (width), 5 positions in y-direction (height) and 4 positions in z-direction (depth). All measurements were carried out using the US-SQS described above. Since only minor variations in the light intensities in xdirection were observed (Figure 1B), values in this direction were averaged (see Table 1). The resulting 20 I 0 values refer to their position in y/z-direction as indicated in Figure 1A. These values were used as surface light intensities of the respective reactor compartments within the modeling approach.
The pH value was controlled by the injection of pure CO 2 at the intake side of the circulation pump using a limit controller (set-point pH 7.25; manipulated variable: solenoid-controlled valve; switch-point: pH 7.25; pH hysteresis: 0.05; on-delay: 5 s). The temperature was controlled using a cooling water circuit supplied by a water bath (F12, Julabo GmbH, Seelbach, Germany) and connected to the double jacket of the system vessel (limit comparator; hysteresis: 1 K). The circulation pump frequency was set to 35.7 Hz resulting in a suspension flow velocity of 0.72 m s −1 in the glass tubes.

Experimental and Modeling Approach
The experimental data used within this work was collected from altogether 15 cultivations which were carried out as a series of discontinuous (batch) bioprocesses in order to avert the lag phase at the start of the individual experiments (see Figure 2). Cultures were harvested partially by draining parts of the suspension and replacing it by fresh, steam sterilized medium (Figure 2). The experiments were designed to be conducted at constant illumination at the reactor surface I 0 but by varying the inoculation biomass concentrations c X,0 , the cultivation cycle times t cyc and the cultivation temperatures T (see Table 2lower part and Figures 5A,B). Temperature variations during the cultivations, however, resulted from a limited cooling capacity of the reactor system under strongly differing ambient conditions.
In order to carry out the modeling, the photobioreactor was subdivided into 20 glass tube compartments of equal size (z = 1, . . . , 20). Light intensity at the respective compartment surface was mapped as shown in Table 1 and remained constant during each individual experiment and across the experiments.
As described by Weise et al. (2019), I 0 is considered to be a sum parameter for the light field description at the respective compartment surface. According to Equation (4) only light absorption is taken into account. Due to the actual reactor features (see Figure 1C), the light field is considered to illuminate the compartment y i z j from only one side with equal intensity. In addition, because of the particular reactor arrangement, as shown in Figure 1C, each physical compartment is represented by two modeled compartments. The resulting photon fluxes (i.e., light intensity multiplied by surface area) within the modeled compartments were added and then related to the total biomass concentration and the reactor volume (see Equation 4).

Data Collection
During the experiments a number of variables were determined on-line and off-line. The on-line values of the optical density at near-infrared OD NIR , the pH and the temperature T were continuously measured and recorded by a 6-channel writer (JUMO Logoscreen 500 cf, Fulda, Germany). The sampling of the off-line values was carried out daily in the morning. The biomass concentration c X was determined as salt-free dry cell weight. Samples for the dry cell weight determination were homogenized, an aliquot of 10 ml was taken, centrifuged in pre-weighted glass tubes at 4.500 * g for 20 min; the resulting supernatant was discarded and the biomass pellet resuspended in 9 ml deionized water; this washing step was carried out twice. The pellets were dried at 105 • C for 24 h and weighed after cooling in a desiccator for 45 min. Measurements were carried out in duplicate.
The original experimental time series data are available as Supplementary Material. FIGURE 2 | View of experimental process data: Experiments No. 05-07 as part of a series of repeated discontinuous cultivations; experiments are separated by partial harvest and medium refill; OD NIR , raw sensor signal OD NIR (black); c X , correlated biomass concentration (green); T, suspension temperature (red); pH, suspension pH value (blue).

Numerical Calculations
All computations were performed using the programming language python (version 3.6.5) and the additional packages astropy (version 3.0.4), numpy (version 1.15.0), pandas (version 0.23.4), SALib (version 1.3.8), scipy (version 1.1.0), and seaborn (version 0.9.0). Table 2-lower part) have been estimated from off-line c X measurements for each individual experiments (original data not shown).

Averaged Approximated Growth Rate
The calculation of the averaged approximated growth ratē µ [d −1 ] is carried out using the approximation of the ordinary differential equation of the biomass concentration for discontinuous bioprocesses (Equation 8). At first, the approximated growth rate µ (t) [d −1 ] is calculated employing the central difference approximation (Equation 9), which is subsequently smoothed using a moving average (Equation 10), where n is the total number of data points i between t − 1 h and t + 1 h. The employed central difference approximation linearizes the growth between the two observations and is therefore only applicable for small t (here: t = 2 h).

Biomass Productivity
The experimental biomass productivity for the respective experiment j was calculated as volumetric yield using the following Equation (11) (data point i, i = 0, n).

Model Validation by MAPE
The criterion used to evaluate the model's accuracy was the Mean Absolute Percent Error (MAPE) (Equation 12) (Mayer and Butler, 1993

Parameter Estimation and Model Analysis
In order to parameterize the developed model, a parameter estimation procedure was carried out in two consecutive steps: The initial parameter estimation, following the methodology proposed by Bernard and Rémond (2012), was applied in a first step to generate initial parameter estimates which were then, in a second step, used as starting values for the final parameter estimation. Both stages of the consecutive parameter estimation were carried out using the same experimental (time series) data. The initial parameter estimation (section 2.6.3) allows to efficiently screen through a wide parameter search range (see section 2.6.2) by fitting the growth kinetics directly to the reorganized data sets (see section 2.6.1). This was carried out in particular with respect to the temperature modeling in order for the unknown cardinal temperatures to converge. Since, however, the reorganization includes approximation and averaging ofμ as well as clustering and selection, the returned results are potentially imprecise.
Therefore, the final parameter estimation (section 2.6.4) was carried out by fitting the model output to the original experimental time series data sets, employing the above mentioned initial estimates as starting values.
Based on the model, parameterized this way, the numerical biomass productivity optimization (section 2.6.6) was carried out with respect to the process parameters c X,0 , t cyc , and T that are to be applied to the process.

Data Set Reorganization
In order to carry out the initial parameter estimation procedure described in the following section, a pooling, reorganization, clustering, and selection of the data was performed: • Calculating the averaged approximated growth ratesμ ,i,j (according to Equation 10) for the i th data point of the j th experiment (j = 1, . . . , 15).
• Concatenating data sets for q ph,i,j ,μ ,i,j and T i,j from all 15 experiments into a new single data set.
• Reassigning the newly formed entries forμ ,i and T i that share the same q ph value to a new k th data set which corresponds to this particular q ph value (k = 1, . . . , s).
• Sorting each data set k obtained by its T i,k values (for identical T i,k values, correspondingμ ,i,k values were averaged).
• Clustering the entire q ph range into 10 equidistant clusters; assigning each data set k to its respective cluster.
• Selecting one data set per q ph cluster that covers the widest individual T k range (k = 1, . . . , 10).
The data sets (k = 1, . . . , 10) reorganized this way were then applied in the first part of the parameter estimation procedure.

Search Range
An adequate search range for the parameter estimation procedure was determined based on reasonable assumptions. In particular, the temperature range was constrained to -5.0, . . . , 50.0 • C, based on the assumption that growth outside this range is not possible for the genus Nannochloropsis.
• T min ranges from −5.0 • C (lower bound) to the lowest temperature for which experimental data is available (upper bound).
• T opt ranges from −5.0 • C (lower bound) to 50.0 • C (upper bound), as long as the model's properties are fulfilled (see Equations 5, 7).
• T max ranges from the highest temperature for which experimental data is available (lower bound) to 50.0 • C (upper bound).
• µ max ranges from the highestμ value obtained from the experiments (lower bound) to the highest µ max value reported for microalgae (upper bound; see Table 1-lower part).
• K S,ph ranges from 0 (lower bound) to the highest q ph value obtained from the experiments (upper bound).
• m ph ranges from 0 (lower bound) to the highest m ph value reported in the literature (upper bound; see Table 1-lower part).

Initial Parameter Estimation
Each data set k, reorganized as described above, contains growth rate valuesμ ,i,k at temperatures T i,k . Bernard and Rémond (2012) published a strategy that consists in the following: Identifying s + 3 parameters µ opt,1 , µ opt,2 , . . . , µ opt,s , T min , T opt , T max (here s = 10, for the 10 cluster data sets) vectorised as θ T , where the parameter µ opt,k is the optimum growth rate at the optimum temperature T opt for the kth data set. Following their strategy, the underlying key assumption is that the cardinal temperatures (T min , T opt , T max ) are common to all the data sets, whilst µ opt,k is constant within the respective data set due to constant experimental lighting conditions q ph,k . The unknown s+3 parameter values were estimated in a first step by minimizing the optimization criterion J(θ T ) (Equation 13) (Bernard and Rémond, 2012). Within a second step, the remaining model parameters µ max , K S,ph and m ph (vectorised as θ q ph ) were estimated in a similar way (Equation 14). Both, Equations (13) and (14), were minimized with respect to the respective θ resulting in the minimum optimization criterion value J( ) for the optimum parameter vector according to Equation (15).
The downhill-simplex method (Nelder-Mead method) used for the model parameter estimation procedure (scipy.optimize.fmin) was started from every possible parameter value combination within the search ranges described above. The T range was screened stepwise by T = 5 • C, resulting in 70 initializations that match the model's properties (Equations 5, 7). The ranges for µ max , K S,ph and m ph were screened using 10 equidistant values along each axis, resulting in 1,000 initializations. The method was set to perform a maximum of 200 iterations for each initialization. Computations were carried out using the standard settings of the above package regarding the convergence criteria of the method.

Final Parameter Estimation
Since the initial parameter estimation described above can provide results which are potentially imprecise due to the approximation and averaging ofμ , the reorganization and clustering of the data sets k as well as the split parameter estimation of θ T and θ q ph , a second and final parameter estimation step was carried out. Based on the ordinary differential equation of the biomass concentration for discontinuous bioprocesses (Equation 16), the original experimental time series data sets j (j = 1, . . . , 15) were simulated using the initial parameter estimates. Numerical integration was carried out employing the Euler method. The time series of c X were simulated for the i discrete time points t i,j and their corresponding temperatures The objective of this final parameter estimation was to identify the m + 6 parameters c X,0,1 , c X,0,2 ,. . . , c X,0,m , µ max , K S,ph , m ph , T min , T opt , T max (here m = 15, for the 15 experiments) vectorised as θ f , where the parameter c X,0,j is the inoculation biomass concentration c X,0 for the jth data set. The optimum parameter values were estimated by minimizing the optimization criterion J(θ f ) (Equation 17). The final optimum parameter vector f (according to Equation 18) was determined in analogy to the initial parameter estimation as described above.
with c X,meas = c X,data,i,j c X,pred = c X,model,i,j (q ph,i,j , T i,j , c X,0,j , µ max , K S,ph , m ph , T min , The Nelder-Mead method was also employed for this second step of the model parameter estimation procedure (scipy.optimize.fmin). The initially estimated parameter vectors were now used as start parameter vectors. The method was set to carry out a maximum of 200 iterations for each initialization. Again, the standard settings of the above package were used with respect to the convergence criteria.
Since the values of the parameters T min and T max to be estimated were expected to lie outside the experimental data range and to account for the relatively small number of experiments, the Delete-2-Jackknife analysis (astropy.stats.jackknife_stats) was used to indicate the 95 % confidence interval of the identified parameters of the vector f (Efron and Tibshirani, 1993;Duchesne and MacGregor, 2001). Jackknife re-sampling was carried out by bootstrapping two respective data sets (experiments) from the data at a time.

Parameter Sensitivity Analysis
A variance-based sensitivity analysis according to Sobol (2001) was performed using the SALib.analyze.sobol in order to investigate the model's parameter sensitivities more deeply. This sensitivity analysis was done with respect to the biomass productivity as the main focus of this study. Within Equation (19), M represents the model output depending on the process parameters c X,0 , t cyc and T as well as the model parameter vector θ (Zhang et al., 2017) consisting of µ max , K S,ph , m ph , T min , T opt , and T max .
The Sobol method is based on the decomposition of the variance of the model output Var(M) (Equation 20) into summands of increasing dimensionality (Zhang et al., 2015(Zhang et al., , 2017. Var i is the partial variance corresponding to the first-order index of θ i of the model output M, while Var ij is the partial variance corresponding to the second-order index of the ith and jth parameter interaction (Zhang et al., 2015(Zhang et al., , 2017. k is the number of model parameters, here k = 6. M = Pr model (c X,0 , t cyc , T, θ ) Var ij + . . . + Var 1,...,k The sensitivity indices S i and S ij (Equations 21, 22) are calculated as ratios of the partial variances to the total variance (Zhang et al., 2015). Total-order indices S Ti are calculated following Equation (23) using Var ∼i which represents the variation of all parameters except θ i (Homma and Saltelli, 1996;Sobol, 2001;Zhang et al., 2017). S i quantifies the effect of varying θ i alone, while S Ti quantifies the effect of varying θ i and includes all effects caused by its interactions with all other model parameters.

First-order index
The Saltelli sampling scheme (SALib.sample.saltelli) was used to generate model parameter samples of the final model parameter estimates f (section 2.6.4) varied by ± 3 σ (see Table 2-upper part). This scheme generates N · (2k + 2) samples (here N = 10,000; k = 5). The calculated indices were classified using thresholds. Indices contributing < 0.01 to the model output variance are considered "non-sensitive, " while indices contributing ≥ 0.1 are considered "highly sensitive" (Tang et al., 2007). S i and S Ti also potentially depend on the process parameters. The experimentally set c X,0 is associated with the present light availability, while the set T may influence the sensitivity of the three cardinal temperatures. Therefore, the sensitivity analysis was carried out for 9 combinations (scenarios) of c X,0 and T taking into account data corresponding to percentages of the respective experimental range [20% ("low"), 50% ("medium"), and 80% ("high")].

Biomass Productivity Optimization
In order to optimize the biomass productivity Pr with respect to the described photobioreactor, a numerical optimization was carried out. The objective was to maximize Pr with respect to the process parameters c X,0 , t cyc , and T. Numerical integration of Equation (16) was performed using the estimated model parameters ( f ; Table 2-upper part) by employing lsoda of the ODEPACK (scipy.integrate.odeint). After each integration, Pr was calculated for the individual iteration steps using Equation (11). The process parameters c X,0 , t cyc , and T were vectorised as θ Pr containing all possible process parameter combinations, while Pr contains the optimum process parameter combination with respect to Pr. The optimum process parameter values were estimated by minimizing the optimization criterion J(θ Pr ) (Equations 24, 25) using also the downhill-simplex method (Nelder-Mead method; scipy.optimize.fmin). With respect to numerical integration and convergence criteria, standard settings of the above packages were used.

RESULTS AND DISCUSSION
The experiments carried out and presented here consisted of a series of discontinuous (batch) cultivations. Single experiments were separated by partial harvest of the biomass suspension and replacement by new medium (see Figure 2). Inoculation biomass concentrations c X,0 and cultivation cycle times t cyc were varied between the experiments (see Figure 2 and Table 2-lower part). Temperature T as well as pH were recorded during the experiments (see Figures 2, 5A,B). The original experimental time series data and the reorganized data sets are available as Supplementary Material. Light saturation was not observed during the experiments. The linear accumulation of biomass indicates limited lighting conditions (Figure 2). The light intensity at the photobioreactor surface varied between 19 and 836 µmol ph m −2 s −1 (average: 337 µmol ph m −2 s −1 ; see Table 1-upper part and Figure 1). These incident light intensities range from low values (100 µmol ph m −2 s −1 Raso et al., 2011) to rather typical values at lab-scale (700 µmol ph m −2 s −1 Pal et al., 2011;up to 850 µmol ph m −2 s −1 Wagenen et al., 2012) for Nannochloropsis sp. With respect to the scalability toward outdoor production, higher incident light intensities of 1,500 to 2,000 µmol ph m −2 s −1 are present on sunny days. These however are lowered due to the vertical arrangement of the glass tubes at industrial scale, the light absorption by greenhouse parts, etc. Therefore, the comparability and the transferability of presented model to industrial-scale outdoor conditions is limited unless further model extensions regarding higher variability in light and temperature, nightly biomass losses, etc. are carried out.
The inoculation biomass concentrations within this study (0.27,. . . , 1.03 g X l −1 ) are within a typical range for both, lab-scale and outdoor production of Nannochloropsis sp. in tubular photobioreactors (Olofsson et al., 2012;Pfaffinger et al., 2016;Pereira et al., 2018). The experimental variability regarding c X,0 is rather high within the investigated small pilotscale experiments. However, biomass concentrations throughout the cultivation cycle time are found to be typical for large-scale and outdoor production of Nannochloropsis. c X ranged 0.27,. . . , 2.79 g X l −1 , reflecting usual inoculation and harvesting biomass concentrations, respectively. To the best of our knowledge, depending on the climatic zone, biomass concentrations of Nannochloropsis outdoor cultures within tubular photobioreactors rarely exceed 3 g X l −1 (Olofsson et al., 2012;Benvenuti et al., 2015) and do not exceed 5 g X l −1 (Zittelli et al., 1999). Therefore, the c X range observed here is comparable to larger-scale and outdoor processes.
Biomass productivity Pr was calculated for each experiment according to Equation (11) in order to evaluate the model's accuracy expressed as MAPE (see Table 2-lower part) and to use it for the biomass productivity optimization. The overall context of this work was to optimize Pr as part of a series of batch bioprocesses (i.e., repeated batch).
The devised two-stage consecutive parameter estimation procedure was carried out using the reactor-specific and strainspecific constants as well as the set parameter constraints shown in Table 1.
The initial parameter estimation converged to two valid solutions ( 1 and 2 , Table 2-upper part) depending on the start parameter combinations. Both solutions differ only strongly in the estimated value of T min .
Both initial parameter estimation vectors ( 1 and 2 ) were used as start parameter combinations for the final parameter estimation. It was performed using the original time series data sets j (in contrast to the reorganized data sets k). The altogether six shared model parameters where estimated within multiexperiment fittings.
The final parameter estimation converged to the solution f (Table 2-upper part, Figure 3) initialized from each of the two initial parameter combinations ( 1 and 2 ). The φ(T) data shown in Figure 3A was normalized with respect to the corresponding estimated µ opt (Figure 3B) according to Equation (1).
The estimated values of T min and T max lie outside the range 19, . . . , 31 • C of available experimental data. Since the range between the estimated T min and the available data is much wider than the one between the estimated T max and the available data, in particular the T min estimation is afflicted with some FIGURE 3 | Model parameter estimation results: Initial parameter estimation provides solutions 1 and 2 and final parameter estimation solution f (see Table 2 uncertainty, which is further enforced by the asymmetric shape of the temperature function φ(T) as can be seen in Figure 3A showing the CTMI results. This is underlined by the jackknife 95 % confidence intervals CI and standard deviations σ of the estimated parameters ( Table 2-upper part). While for all other parameters the σ values indicate relatively certain parameter estimations, σ for T min is at the same order of magnitude as the parameter's value.
The finally estimated cardinal temperatures for N. granulata (T min = 2.3 • C, T opt = 27.93 • C, T max = 32.59 • C) lie within the specified parameter constraint ranges ( Table 1) and match laboratory experience. The CTMI parameters estimated by Bernard and Rémond (2012) for the related species N. oceanica (data by Sandnes et al., 2005) showing a good degree of similarity (T min = −0.2 • C, T opt = 26.7 • C, T max = 33.3 • C).
The results obtained here are supported by an optimum growth temperature range of 20, . . . , 29 • C given in the literature (Sukenik, 1991;Wagenen et al., 2012;Bartley et al., 2015;Abirami et al., 2017). Further, Wagenen et al. (2012) reported only very poor growth below 13.6 • C and above 32.3 • C for N. salina. Also, Sukenik (1991) found that Nannochloropsis spp. failed to grow below 10 • C and above 38 • C. In view of the fact that the parameter estimation carried out here had only rough quantitative specifications with regard to the cardinal temperatures, the values calculated (Table 2-upper part) are in strong compliance with the literature. Sandnes et al. (2005) empirically described a light-dependent optimum cardinal temperature T opt for N. oceania (Sandnes et al., 2005). Although this combined effect of light and temperature is physiologically plausible, the functional relationship of these two variables was not observed for N. granulata within the own data.
The µ opt values for 1 and 2 estimated increase monotonically with higher q ph (Figure 3B). These estimates are subject to fluctuations but show no inhibition within the observed q ph range. Therefore, a growth kinetics considering light limitation (Equation 2) without light inhibition was applied. The q ph range observed during the experiments (0.59, . . . , 1.9 mol ph g X −1 d −1 ) is comparable to q ph values given by the literature (0.25, . . . , 2.2 mol ph g X −1 d −1 ) (Zijffers et al., 2010;Kandilian et al., 2014;Janssen et al., 2018). In accordance with Equation (4) q ph is calculated as a mean light availability rate of the whole culture suspension, considering the total culture volume V L . However, local light intensities differ due to light gradients within the culture suspension, which are not spatially resolved within the model. Despite this aggregation, the devised model predicts microalgal growth and biomass productivities precisely. The estimated photon maintenance coefficient m ph = 0.346 mol ph g X −1 d −1 is comparable to values reported for C. sorokiniana: 0.16 mol ph g X −1 d −1 and D. tertiolecta: 0.32 mol ph g X −1 d −1 (Zijffers et al., 2010). Pirt (1986) reported increasing m ph values for photobioreactor set-ups with lower illuminated/non-illuminated culture volume ratios, which has to be considered during the transfer to different photobioreactor set-ups. The growth yield Y X,ph [g X mol ph −1 ], calculated following, Pirt (1986) as Y X,ph = µ opt / (q ph − m ph ), varied 0.45, . . . , 0.73 g X mol ph −1 during the experiments (based on Figure 3B), which matches the Y X,ph range of 0.2, . . . , 2.1 g X mol ph −1 reported in the literature (Zijffers et al., 2010;Dillschneider et al., 2013;Schediwy et al., 2019). The estimated parameter µ max = 1.56 d −1 is comparable also to µ max ranges between 0.86, . . . , 1.6 d −1 as found in the literature (Sandnes et al., 2005;Spolaore et al., 2006;Weise et al., 2019).
The similar parameters indicate that the model structure in particular with respect to the light-limited growth modeling shows a good transferability for tubular photobioreactors that differ in scale, glass tube diameter and the mode of bioprocess operation, although the experimental and modeling approaches differed in many respects. In contrast to Weise et al. (2019), who considered constant temperature conditions (T = 22 • C) under steady-state turbidostatic bioprocess operation, experiments presented here were conducted under variable temperature conditions and have been scaled up by ≈ 7x regarding the reactor liquid volume V L , by 10x regarding the number of modeled compartments and by 2x regarding the glass tube diameter.
In order to transfer the devised model to other reactors and different scales, it is necessary to consider homogeneous onesided illuminated glass tube compartments within the model, although the glass tubes may be illuminated homogeneously from both sides in the actual arrangement. Therefore in this context also, the physical compartments are to be decomposed mathematically into several theoretical ones. If these assumptions do not apply to the specific photobioreactor structure, the model cannot be transferred without further modification. The presented model considers basic geometric properties (e.g., inner tube radius, liquid reactor volume, illuminated reactor surface, etc.) of tubular photobioreactors. Using biomass productivity or reactor size as scale-up objectives, it is possible to estimate the dimensions (e. g. illuminated reactor surface) of a photobioreactor scale-up.
The presented model is implemented using the ordinary differential equation (ODE) for batch bioprocesses. In order to adapt the model for continuous cultivation, the ODE for batch bioprocesses can be extended by a dilution term. The growth kinetics, which are described by algebraic equations are not affected by this extension.
With respect to the model's sensitivity analysis, Figure 4 illustrates the first-order and total-order sensitivity indices of the model parameters across different scenarios. The comparatively small differences in the values of S i and S Ti indicate only limited interactions between the parameters, which was to be expected due to the model's structure (see Equations 1, 2, 7). The output of the model is highly sensitive (S i and S Ti ≥ 0.1) with respect to the parameters µ max and K S,ph . Both parameters combined contribute predominantly to the variance in the model output, since they describe the utilization of light as the sole energy source for phototrophic growth. The set c X,0 effects the sensitivities of K S,ph and m ph . Since higher c X,0 reduce the q ph , the photon half-saturation constant K S,ph becomes more sensitive. In addition, m ph becomes sensitive (S i and S Ti ≥ 0.01) only at higher c X,0 (Figure 4).
T opt is sensitive (S i and S Ti ≥ 0.01) in all scenarios, while T min and T max are sensitive only at lower and higher temperatures. Also, T min and T max become non-sensitive at temperatures close to T opt . At temperatures away from T opt , the combined sensitivity of the three cardinal temperatures contributes up to 10 % (i.e., 0.1) to the model's output variance. The set c X,0 does not influence the sensitivity of the cardinal temperatures.
All experiments were carried out in the range 20, . . . , 25 • C, except Experiment No. 10 at 31 • C. Although the temperature inside the photobioreactor was actively controlled, diurnal rhythms in the ambient temperature resulted in temperature variations within individual experiments. These are displayed in Figure 5A using kernel density approximations (violins). The width of the individual violin represents the number of available data at the respective temperature value, while small pin-like tails of the violins indicate short-term outliers of the measured variable. In analogy, Figure 5B shows the distribution of the measured pH values of the respective experiments. Variations in temperature and pH value did not occur rapidly, as can be seen in Figure 2. Despite the variations in the observed temperature and pH ranges, using the specific approach applied here, the available data could be successfully processed in order to estimate model parameters and optimize process parameters at pilotscale. Figure 5C provides violin plots that represent the absolute modeling error with respect to the biomass concentration time series using the estimated parameter vector f . It can be seen that the modeling error rarely exceeds ± 0.15 g X l −1 .
Furthermore, one good as well as one badly predicted experiment are highlighted in green (Experiment No. 05) and red (Experiment No. 13), as shown in Figure 5C. These two experiments are presented as time series in Figures 5D,E. Although the biomass concentration has been well predicted over much of the time (Figures 5D,E), some regions are underestimated or overestimated at the beginning and the end of the respective experiments ( Figure 5E as example).
The prediction accuracy MAPE (Equation 12) of the target value biomass productivity Pr (Equation 11) improved over the course of the estimation procedures. Table 2-lower part provides an overview of these experimental and modeling results as well as the parameters a and b. The initial parameter estimation of 1 and 2 resulted in MAPE( 1 ) = 7.5% and MAPE( 2 ) = 7.2%, whilst the prediction accuracy after the final parameter estimation was MAPE( f ) = 7.2%. Therefore, the accomplished overall model prediction accuracy is very satisfactory, although the Absolute Percent Error of individual experiments (APE j ) covers a wider range 0.1, . . . , 26.1%. Typical MAPE values for industrial and business data and their interpretation are: < 10% highly accurate forecasting, 10, . . . , 20% good forecasting, 20, . . . , 50% reasonable forecasting, > 50% inaccurate forecasting (Lewis, 1982;Moreno et al., 2013).
Based on the results of the parameter estimation procedure (Table 2-upper part), the parameterized model has been used to design an optimized bioprocess regime regarding the inoculation biomass concentration c X,0 , the cultivation cycle time t cyc and the temperature T with respect to optimum biomass productivity Pr opt for a series of equal batch cultivations (repeated batches). Since, the specific light availability rate q ph decreases monotonously along with biomass accumulation under constantly illuminated discontinuous bioprocess operation, q ph is only influenced by alterable c X,0 and t cyc .
Following the above objective, a numerical optimization of Pr was conducted using Equations (24) and (25) as described in section 2.6.6. The biomass productivity optimization showed two major aspects: First, the optimum growth temperature equals the model parameter T opt = 27.93 • C, which matches the intuitive expectation with regard to Equations (1) and (5). This applies independently of the lighting conditions according to the devised model. Second, the optimization results in a corner point optimum that points toward an insignificantly small cultivation cycle time t cyc → 0.
The developed approach therefore predicts the theoretical optimum biomass productivity Pr opt → 0.50 g X l −1 d −1 within the investigated bioreactor set-up for negligibly small cultivation cycle times with c X,0,opt → 1.07 g X l −1 , and hence the transition from discontinuous (e.g., batch, fed-batch or their repeated versions) to continuous (e.g., chemostat, turbidostat) bioprocess operation. The optimum inoculation biomass concentration c X,0,opt calculated this way corresponds to the steady-state biomass concentration of a settled chemostat bioprocess or to the optimum biomass concentration set-point of a turbidostat bioprocess, respectively (Weise et al., 2019). (Statement is no longer applicable due to the modifiaction of Equation 2). The above findings are also supported by Chen et al. (2012) who estimated c X,0,opt = 0.98 g X l −1 (Pr opt = 0.75 g X l −1 d −1 ) for a small-scale laboratory set-up, as well as Hulatt et al. (2017) finding Pr opt = 0.51 g X l −1 d −1 for a flat-panel photobioreactor, both under discontinuous operation. In addition, literature values for Pr at pilot-scale and industrial-scale show a range < 0.1, . . . , 0.71 g X l −1 d −1 (de Vree et al., 2015;Pereira et al., 2018), confirming the feasibility of the Pr opt values obtained here.
A simulation of the original experiments was carried out using the actually applied cultivation cycles times t cyc of the experiments (see Table 2-lower part) as well as the optimized process parameters c X,0,opt and T opt (according to section 2.6.6).
TABLE 3 | Numerical optimization results: Model predicted Pr opt for different bioprocess operations and pre-set t cyc at T opt with the required c X,0,opt (see Figure 6A).

Bioprocess operation
t cyc c X,0,opt Pr opt Discontinuous ( The model predicts a noticeable increase of Pr by 39% from 0.33 g X l −1 d −1 (see Table 2-lower part) to 0.46 g X l −1 d −1 . The convergence toward the predicted optimum biomass productivity is asymptotic. This results in only minor additional gains in biomass productivity when moving from short-cycled discontinuous to continuous bioprocesses and requires a simultaneous adaptation of the inoculation biomass concentration c X,0 . Within this context, Table 3 shows examples of the predicted Pr opt under discontinuous (repeated batch) bioprocess operation for fixed cultivation cycle times t cyc (1 d, 3.5 d, and 7 d) starting at the required inoculation biomass concentration c X,0 and under optimum temperature T opt . As a consequence, solely reducing t cyc from 7 d to 3.5 d would outweigh the increase in operational efforts by a higher Pr (0.37 → 0.46 g X l −1 d −1 ; with c X,0,opt 0.17 → 0.43 g X l −1 ) under discontinuous bioprocess operation. Figure 6A illustrates these relations in more detail by presenting the model's prediction with respect to biomass productivity Pr depending on the cultivation cycle time t cyc and the inoculation biomass concentration c X,0 for the optimum cultivation temperature T opt . According to the model's prediction, only moderately high biomass productivities are gained for cultivation cycle times greater ≈ 5 days, as well as for inoculation biomass concentrations above ≈ 1 g X l −1 (Figure 6A). Notably, the optimum productivity is reached at the shortest cultivation cycle time, which illustrates the said corner point optimum regarding the cultivation cycle time t cyc .
In accordance with these findings, Figures 6B-D additionally provides graphical representations of the biomass productivity loss when varying the process parameters c X,0 , t cyc and T. Figures 6B-D illustrate scenarios for three different cultivation cycle times: 1 d, 3.5 d (1/2 week), and 7 d (1 week) with respect to discontinuous bioprocess operation (repeated batch). Considerations regarding lag phases are neglected here, as can be seen in Figures 2, 5D,E. All three scenarios show a single optimum of biomass productivity at the respective optimum temperature T opt and the optimum inoculation biomass concentration c X,0,opt .
The predicted Pr opt increases by about 35% from 0.37 g X l −1 d −1 at t cyc = 7 d (c X,0,opt = 0.17 g X l −1 ) to 0.50 g X l −1 d −1 at t cyc = 1 d (c X,0,opt = 0.84 g X l −1 ). In general, the optimum biomass productivity Pr opt increases when the cultivation cycle time is shortened. In addition, when this time is shortened, a higher c X,0,opt is required to obtain the optimum productivity.

CONCLUSION
The work presented here provides a transferable methodology to model microalgal growth covering light availability and temperature based on experimental data from cultivation runs in a small pilot-scale tubular photobioreactor under discontinuous operation in order to subject it to biomass productivity analysis and optimization.
The established model with its estimated parameters accurately predicts light and temperature dependent growth of Nannochloropsis granulata. The parameter ranges are supported by the literature. In general, the model parameter T opt is much closer to T max than to T min , thus the CTMI displays a strong asymmetry. Temperatures above T opt therefore lead to a steep decline in the growth rate and also the biomass productivity. Since an accurate temperature control is hardly to provide under large-scale or outdoor conditions, these processes should be operated below the targeted T opt .
Model-based numerical biomass productivity optimization for repeated discontinuous operation points toward best performance under continuous operation. The optimization for repeated discontinuous operation yields reduction of cultivation cycle time and increase of inoculation biomass at optimum temperature. The calculated optimum inoculation biomass concentrations c X,0,opt and the corresponding optimum biomass productivities Pr opt are confirmed by different publications. Furthermore, biomass productivities for laboratory-scale, pilotscale and industrial-scale reported in the literature support the feasibility of the Pr opt values obtained here. Applying these optimized process parameters would deliver a noticeable increase in biomass productivity.
The successful application of this approach here to small pilot-scale under discontinuous operation, following a previous investigation into lab-scale under continuous operation, indicates its potential transferability also to larger scale tubular photobioreactors covering both light and temperature dependent microalgal growth and biomass productivity.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
TW contributed to the conception and design of the study, the collection and assembly of data, the analysis and interpretation of the data, the drafting, critical revision and final approval of the article. CG contributed to the collection and assembly of data, the drafting, critical revision and final approval of the article and to obtaining funding. MP contributed to the conception and design of the study, the interpretation of the data, the drafting, critical revision and final approval of the article and to obtaining funding.