# Adaptive temperature model for microalgae cultivation systems

Ali Gharib, Walid Djema, Francesca Casagli, Olivier Bernard

► To cite this version:

Ali Gharib, Walid Djema, Francesca Casagli, Olivier Bernard. Adaptive temperature model for microalgae cultivation systems. Journal of Process Control, 2024, 141, pp.103280. 10.1016/j.jprocont.2024.103280. hal-04868166

HAL Id: hal-04868166

https://inria.hal.science/hal-04868166v1

Submitted on 6 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Adaptive Temperature Model for Microalgae Cultivation Systems

A. Gharibª, W. Djemaª, F. Casagliª, O. Bernardª

ªINRIA Sophia Antipolis Méditerranée, BIOCORE Project-Team, Université Nice Côte d'Azur, 2004, Route des Lucioles BP 93, Valbonne, 06902, Sophia-Antipolis Cedex, France

## Abstract

Microalgae cultivation for energy production is a promising avenue for converting solar light into sustainable biofuel. Solar processes are however subjected to the permanent fluctuations of light and medium temperature. Accurate temperature prediction of the culture medium turns out to be critical for optimising growth conditions. In this study, we introduce a reduced-model approach derived from existing models, turning the complex heat transfer modelling problem into an identification problem. The resulting generic model, called the Simplified Auto Tuning Heat Exchange (SATHE) model, has a clear and simple structure, offering a balance between accuracy and computational complexity. The SATHE model is versatile and contains the necessary terms to catch a large variety of heat transfer problems, while the parameters can be identified from experimental data. We first prove the parameter identifiability and then propose an identification strategy, based on the gradient computation, to identify the model's underlying parameters. We further validate the SATHE model performance in two distinct reactors across various seasons. Finally, we discuss the potential of online applications with a continuous self-tuning strategy to keep optimal predictive performances. This work lays the foundation for enhanced control strategies in large-scale cultivation systems.

**Keywords:** Modeling, Temperature, Identifiability, Nonlinear, Dynamical systems, Microalgae.

## 1. Introduction

Microalgae are microscopic, autotrophic organisms found in both marine and freshwater environments. These organisms play a significant role in our ecosystem and are emerging as promising resources in various industrial sectors. From serving as a sustainable feedstock for biofuel production to their use in wastewater treatment, pharmaceuticals, and health foods, microalgae offer remarkable potential in numerous applications [1, 2].

The growth and productivity of microalgae depend on multiple factors, including physical conditions like light, nutrients, and notably, temperature. Medium temperature greatly affects their metabolism, photosynthesis, and overall growth [3]. An accurate heat transfer model for accurately predicting the temperature along various climatic conditions is therefore crucial for managing and optimising microalgae cultivation.

In the field of microalgae cultivation, existing models like [4] offer accurate predictions of temperature in outdoor raceway ponds, although their precision can wane in colder seasons. Casagli et al. [5] refined this model, achieving higher prediction accuracy. The model of [6] is also interesting, but requires the measurements of the soil temperature below the raceway to predict the temperature in the pond. It can therefore not be used for medium temperature prediction. Pruvost et al. [7] presented a thermal model discredited into layers validated for a period of 4 days but with a high level of complexity. The literature also contains simplified models that are designed for particular for closed [8] or open reactors [9, 10], validated for periods shorter than a few weeks.

However, a drawback of these models is their rigidity, mak-

ing them difficult to adapt to new reactor configurations and designs, and to prove their validity along various seasons. Customising models to adapt them to a new reactor in another context is a long and tedious task.

In recent work [11], we introduce the SATHE model, an adaptable reduced-model that can be easily tailored to different configurations, reducing the expertise required, and transforming the modelling problem into a parameter identification problem. The framework in Fig. 1 outlines our approach, building on the previous work [11]. We use a four-day dataset of measured data (weather data, fluxes, and reactor temperature) to fine-tune the SATHE model. Then, we employ this refined model to predict temperature and water depth in raceway ponds based on weather forecast data. The novelty of this model and the associated identification approach is

* the strategy for adapting the parameters to a new reactor and climatology using a small amount of data

* the regular parameter re-adaptation, ensuring validity periods of more than one month even for different climatologies

* the reduction in complexity, which makes it more convenient for model-based control purposes

This paper builds on the work presented in [11] and provides an in-depth description of the methodology used to derive a reduced-model from the heat transfer model developed by Casagli et al. [5] for a raceway reactor located in Narbonne (France), and made of concrete laying on the ground. Furthermore, we investigate the identifiability of the reduced-model to
<figure><img src="image_2.png" /><figcaption>Figure 1: Schematic framework for predicting temperature, depth, and microalgae growth in a bioreactor</figcaption></figure>

ensure accurate parameter estimation.

Meteorological conditions influence reactor temperature even more in smaller systems and shallow raceways due to their lower thermal inertia [5, 4, 12]. Consequently, we demonstrate the efficacy of the adaptive reduced-model in predicting the performance of a second reactor, for a different configuration located in Milan (Italy) with a raised raceway made of polypropylene.

Even advanced heat transfer models, such as those proposed by [4, 5], exhibit limitations in accurately predicting temperatures during winter. In this study, we thoroughly examine the predictive capabilities of the developed reduced-model across all seasons for two different reactors.

Sophisticated control strategies, including optimal control and model predictive control, are more efficient when using accurate models and precise predictions for computing future control inputs [13]. Accurate predictions are paramount, directly influencing the efficacy of control strategies. In our study, we evaluate the predictive performance of the reduced-model, emphasising its pivotal role in online applications. This assessment is further strengthened by a comparative analysis with an alternative model.

Through this demonstration, we validate the adaptability and performance of the reduced-model across different seasons, reactor sizes, conductivity with the ground and configurations. This innovative approach can enhance temperature control in microalgae cultivation systems, improving growth rates and increasing productivity.

The paper is organised as follows. In Section 2, we present the mathematical model for microalgae growth and an established temperature prediction inside a reactor. In Section 3 we detail how the reduced adaptive model is derived from established models. Section 4 is dedicated to the identifiability and identification of the parameters in the reduced-model. Sections 5 and 6 study the analysis of the model's performance in two distinct reactors across different seasons. Finally, in Section 7, we discuss the potential application of the model in online modern control approaches.

## 2. Modelling Microalgae Growth and Raceway Pond Temperature Prediction

### 2.1. The general model describing microalgae dynamics

Microalgae growth is influenced by chemical, biological, and notably, physical constraints. Here, we focus on the influence of physical factors, given their important role in growth. Therefore, we consider a model that describes the growth dynamics based on physical parameters assuming that all the other conditions, such as nutrients, are not limiting growth. The average growth rate $\bar{\mu}(I_0, T, x_m)$ is determined by integrating the specific growth rate $\mu(I(z), T, x_m)$, following [12], over the cultivation depth $L$

$$ \bar{\mu}(I_0, T, x_m) = \frac{1}{L} \int_{0}^{L} \mu(I(z, x_m), T) dz, \quad (1) $$

where $x_m$ represents the microalgae biomass, $\mu$ is the specific growth rate dependent on light intensity $I$, and temperature $T$. The light intensity $I(z, x_m)$ follows the Beer-Lambert Law [14], which describes the light attenuation along the depth of the cultivation system. It is defined as

$$ I(z) = I_0 \exp(-\xi x_m z), \quad (2) $$

where $\xi$ represents light attenuation rate per biomass unit, and $z \in [0, L]$. The overall model, proposed in [15], is given by

$$ \left\{ \begin{aligned} \dot{x}_m &= \bar{\mu}(I, T, x_m)x_m - \frac{q^{in}}{V}x_m, \\ \dot{T} &= \mathcal{T}(T, T_a, I, q^{in}, \dots), \\ \dot{V} &= \mathcal{V}(T, T_a, q^{in}, q^{out}, \dots), \end{aligned} \right. \quad (3) $$

where the term $q^{in}$ denotes the inflow rate, $q^{out}$ the outflow rate, $T_a$ the ambient air temperature, and $V$ represents the volume of the cultivation system.

Predicting the temperature evolution $T$ and volume $V$ is crucial for optimising microalgae growth in cultivation systems. In the upcoming section, we analyse the heat transfer model proposed in [16], which describes the evolution of temperature and depth in a raceway pond and we extract the underlying model structure.
## 2.2. The full-model of temperature prediction inside raceway ponds

In this section, we introduce the equations describing the heat transfers in the system, previously documented in [5]. Subsequently, this model will be referred to as the full-model throughout the paper. The model comprises three states that describe the heat transfer within the raceway, the changes in raceway depth, and the temperature dynamics in the soil. The full-model is described by

$$ \rho_w c_{p_w} \frac{d(VT_p)}{dt} = Q_{ra,p}(T_p) + Q_{ra,s}(I_0) + Q_{ra,a}(T_a) + Q_{h,rain}(T_a, T_p) + Q_{h,evap}(\eta, T_a, T_p, v_0) + Q_{h,in}(T_{in}, T_p) + Q_{h,out} + Q_{conv}(T_a, T_p, v_0) - Q_{conv,cond_w,co}(T_p, T_{co}) \quad (4) $$

$$ S \frac{dh_L}{dt} = -S \frac{q_w}{\rho_w} + S q^r + q^{in} - q^{out} \quad (5) $$

$$ \rho_{co} V_{co} c_{p,co} \frac{dT_{co}}{dt} = Q_{conv,cond_w,co}(T_p, T_{co}) + Q_{cond} \quad (6) $$

where $T_p$ (K) is the reactor temperature, $\rho_w$ (kg m⁻³) is the water density, $c_{p_w}$ (J kg⁻¹ K⁻¹) is the specific heat capacity of the water, $S$ (m²) is the reactor surface, $h_L$ (m) the pond depth and $T_{co}$ (K) is the temperature of the concrete (co), which is the reactor support in this case. This model uses the data from observed weather data, which are the light intensity $I_0$ (W m⁻²), the relative humidity $\eta$ (-), the wind speed $v_0$ (m s⁻²), the rain flow $q^r$ (m³ m⁻² s⁻¹) [16]. A list of the heat terms can be found in Tab. 1.

<table><caption>Table 1: Description of Heat Transfer Terms</caption><thead><tr><th>Term</th><th>Description</th><th>Unit</th></tr></thead><tbody><tr><td>Q<sub>ra,p</sub></td><td>Radiation from reactor surface</td><td>W</td></tr><tr><td>Q<sub>ra,s</sub></td><td>Total solar radiation</td><td>W</td></tr><tr><td>Q<sub>ra,a</sub></td><td>Radiation from air to water</td><td>W</td></tr><tr><td>Q<sub>h,evap</sub></td><td>Evaporation flux</td><td>W</td></tr><tr><td>Q<sub>conv</sub></td><td>Convective flux at water surface</td><td>W</td></tr><tr><td>Q<sub>h,in</sub></td><td>Heat flux from influent water</td><td>W</td></tr><tr><td>Q<sub>h,out</sub></td><td>Heat flux to effluent water</td><td>W</td></tr><tr><td>Q<sub>h,rain</sub></td><td>Heat flux from rain</td><td>W</td></tr><tr><td>Q<sub>conv,cond_w,co</sub></td><td>Conductive/convective heat exchange between water and concrete</td><td>W</td></tr><tr><td>Q<sub>cond</sub></td><td>Conductive heat flux with ground at reactor bottom</td><td>W</td></tr><tr><td>q<sub>w</sub></td><td>Evaporation rate</td><td>kg s<sup>-1</sup> m<sup>-2</sup></td></tr><tr><td>P(ζ)</td><td>Saturated vapour pressure</td><td>Pa</td></tr><tr><td>h<sub>conv</sub></td><td>Convection coefficient</td><td>W m<sup>-2</sup> K<sup>-1</sup></td></tr></tbody></table>

In this work, we focus on the terms $Q_{ra,p}$, $Q_{ra,s}$, $Q_{ra,a}$, $Q_{h,in}$, $Q_{h,evap}$, and $Q_{conv}$, as these terms will be reused in Section 3. The radiation from the reactor surface $Q_{ra,p}$, solar radiation $Q_{ra,s}$, air radiation $Q_{ra,a}$ and the heat flux associated with the influent water $Q_{h,in}$ are calculated by

$$ Q_{ra,p} = - \varepsilon_w \sigma T_p^4 S \quad (7) $$

$$ Q_{ra,s} = (1 - f_a) I_0 S \quad (8) $$

$$ Q_{ra,a} = \varepsilon_w \varepsilon_a \sigma T_a^4 S \quad (9) $$

$$ Q_{h,rain} = \rho_w c_{p,w} q^{rain}(T_a - T_p) \quad (10) $$

$$ Q_{h,in} = \rho_w c_{p,w} q^{in}(T_{in} - T_p) \quad (11) $$

where $\varepsilon_w$ (-) is the emissivity of water, $\sigma$ (W m⁻² K⁻⁴) the Stefan-Boltzmann constant, $S$ (m²) the reactor surface, $f_a$ (-) photosynthetic efficiency, and $\varepsilon_a$ (-) the emissivity of the air. The evaporation term $Q_{h,evap}$ in (13) and the convective flux $Q_{conv}$ in (15) both rely on the flow regime. The flow is considered as laminar when the Reynolds number $Re_L < 3 \times 10^5$, or turbulent when $Re_L > 5 \times 10^5$. The Reynolds number is characterised as

$$ \begin{aligned} Re_L &= \frac{L v_z}{\nu_a} \\ v_z &= v_0 (z_r / z_0)^\alpha \end{aligned} \quad (12) $$

where $v_z$ (m s⁻¹) is the wind velocity at the height $z_r$ (m) above the ground with a known wind velocity $v_0$ (m s⁻¹) at the height $z_0$ (m), $L$ (m) is the characteristic length of the raceway, $\nu_a$ the air kinematic viscosity (m² s⁻¹), and $\alpha$ is the power law exponent and is a function of obstructions surrounding the pond. The evaporation flux $Q_{h,evap}$ is calculated by

$$ \begin{aligned} Q_{h,evap} &= -q_w L_w S \\ q_w &= K(P(T_p) - \eta \times P(T_a)) \frac{M_w}{R} \\ P(\zeta) &= \lambda_8 \exp(\lambda_9(\zeta + \lambda_{10})^{0.5}) / \zeta \end{aligned} \quad (13) $$

with the water latent heat $L_w$ (J kg⁻¹), $M_w$ the water molecular weight (kg mol⁻¹), and the ideal gas constant R (Pa m³ mol⁻¹ K⁻¹). The mass transfer coefficient K (m s⁻¹) is depending on

$$ K = S h_L D w_a L \quad (14) $$

$$ \begin{aligned} Sh_L &= \begin{cases} S h_L^a = 0.628(Re_L)^{0.5}(Sch_L)^{1/3}, & \text{for } Re_L < 3 \times 10^5 \\ S h_L^b = 0.035(Re_L)^{0.8}(Sch_L)^{1/3}, & \text{for } Re_L > 5 \times 10^5 \\ S h_L^a(2.5 - 5 \times 10^{-6}Re_L) & \text{for } 3 \times 10^5 \le Re_L \\ +S h_L^b(5 \times 10^{-6}Re_L - 1.5) & \le 5 \times 10^5 \end{cases} \\ Sch_L &= v_a / D w_a \end{aligned} $$

where $v_z$ (m s⁻¹) is the wind velocity at the height $z$ (m) with a known wind velocity $v_0$ (m s⁻¹) at the height $z_0$ (m). $Dw_a$ (m² s⁻¹) is the mass diffusion coefficient of water vapour in air. Sherwood number $Sh_L$ (-) and the Schmidt number $Sch_L$ (-) are both dimensionless. The following term represents the convective flux at the water surface

$$ Q_{conv} = h_{conv}(T_a - T_p)S \quad (15) $$

$$ h_{conv} = N u_L \lambda_a / L $$

$$ N u_L = \begin{cases} N u_L^a = 0.628(Re_L)^{0.5}(Pr)^{1/3}, & \text{for } Re_L < 3 \times 10^5 \\ N u_L^b = 0.035(Re_L)^{0.8}(Pr)^{1/3}, & \text{for } Re_L > 5 \times 10^5 \\ N u_L^a(2.5 - 5 \times 10^{-6}Re_L) & \text{for } 3 \times 10^5 \le Re_L \\ +N u_L^b(5 \times 10^{-6}Re_L - 1.5) & \le 5 \times 10^5 \end{cases} $$
$$Pr = v_a / \alpha_a$$

where $\lambda_a$ (W m⁻¹ K⁻¹) is the air thermal conductivity, $Pr$ (–) the Prandtl number, and $\alpha_a$ (m² s⁻¹) is the air thermal diffusivity. $Nu_L$ (–) is the Nusselt number, which is also dependent on the Reynolds number $Re_L$. For further details, we refer to [4, 16].

The full-model is complex, involving numerous parameters and demanding an in-depth understanding of the specific reactor configuration it applies to. In the next section, we will detail our approach to extract the model structure so as to derive a simplified reduced-model (the SATHE model) from this comprehensive full-model.

### 3. Reducing Complexity: Deriving the Reduced-Model for Temperature Prediction

In our pursuit of a more streamlined, generic and efficient model for temperature prediction, we start with a process of simplification. This involves identifying the key terms that have the most influence on the heat balance within the system. By focusing on these critical components, we aim to create an accurate reduced-model while minimising complexity. According to Bechet et al. [4] and to Casagli et al. [16], the most relevant terms in the heat balance are $Q_{ra,p}$, $Q_{ra,s}$, $Q_{ra,a}$, $Q_{h evap}$, and $Q_{conv}$. Focusing only on these terms, and simplifying their parametrization, results in a much simpler structure involving fewer parameters. Then, the other terms have structures which are the sum of slowly varying terms (for example due to the evolution of the ground temperature for the conduction flux), or that are linear with respect to $T_p$ (like the conducto-convective flux between water and concrete) which will be compensated by the linear term in the generic model.

In Table 2, we compare the representation of the various terms in both models. The subsequent section outlines the simplification process applied to the terms of the full-model. The resulting reduced-model is shown in (24), which includes the adaptive parameters $\theta_i$ and constant values $\lambda_i$, which result from the simplification and multiplication of constants of the most relevant terms.

<table><caption>Table 2: SRP (Specific to raceway pond), CS (compensated in the SATHE model), AAR (adaptive to any reactor)</caption><thead><tr><th>Term</th><th>full-model</th><th>SATHE model</th></tr></thead><tbody><tr><td>Q<sub>ra,p</sub></td><td>SRP</td><td>AAR</td></tr><tr><td>Q<sub>ra,s</sub></td><td>SRP</td><td>AAR</td></tr><tr><td>Q<sub>ra,a</sub></td><td>SRP</td><td>AAR</td></tr><tr><td>Q<sub>h evap</sub></td><td>SRP</td><td>AAR</td></tr><tr><td>Q<sub>conv</sub></td><td>SRP</td><td>AAR</td></tr><tr><td>Q<sub>h,in</sub></td><td>SRP</td><td>AAR</td></tr><tr><td>Q<sub>h,rain</sub></td><td>SRP</td><td>SRP</td></tr><tr><td>Q<sub>conv,cond,w,co</sub></td><td>SRP</td><td>CS</td></tr><tr><td>Q<sub>cond</sub></td><td>SRP</td><td>CS</td></tr><tr><td>q<sub>w</sub></td><td>SRP</td><td>AAR</td></tr></tbody></table>

We reparametrised the terms $Q_{ra,p}$ (7), $Q_{ra,s}$ (8), and $Q_{ra,a}$ (9) in a similar manner. These terms share the same structure consisting of (physical) constants multiplied by either the state $T_p$ or a weather input value ($T_a$ or $I_0$). The multiplication of these (physical) constants yields to the values $\lambda_i$.

$$ 
\begin{aligned}
B(T_p) &= \lambda_2 T_p^4 \\
D(T_a) &= \lambda_4 T_a^4 \\
E(I_0) &= \lambda_5 I_0,
\end{aligned} 
\quad (16) $$

Since these constants vary across cultivation systems, we multiply these terms with parameters $\theta_i$, which will be identified (see Sec. 4). To illustrate, we provide an example of how the constants $\varepsilon_w$, and $\sigma$ in the pond radiation term $Q_{ra,p}$ in (7) yield to $\lambda_2$.

$$ 
\begin{aligned}
Q_{ra,p} &= -\varepsilon_w \cdot \sigma \cdot S \cdot T_p^4 \\
&= S \lambda_2 T_p^4
\end{aligned} 
\quad (17) $$

Since surface $S$ is a key parameter in biomass production, we decided to keep this parameter explicit in our parametrization where it appears in the denominator for the temperature dynamics (see Appendix B). A comprehensive list of the $\lambda_i$ values for the SATHE model is available in Table 3.

In the full-model, (15) represents the convective heat transfer $Q_{conv}$ between the air and the water in the raceway pond. This term depends on the Nusselt number $Nu_L$, which varies for different cases, as it is dependent on the Reynolds number and so on the flow regime. However, it turns out that the change in the Nusselt number has a neglectable impact on temperature prediction. To simplify the model parametrisation, we propose a nominal parameter value with a constant Nusselt number based on the Reynolds number (12). Here, we consider the wind velocity to be constant. This simplification is based on an average wind velocity of $\bar{v}_0 = 5.04 \text{ m s}^{-1}$, calculated from a dataset spanning 443 days. This assumption results in Reynolds numbers consistently exceeding $4 \times 10^5$, indicating a state of turbulent flow for the specific parameters. Consequently, we use as a nominal parameter, the Nusselt number given by equation (15) for turbulent flow which leads to:

$$ \tilde{Q}_{conv} = S \lambda_3 (T_a - T_p). \quad (18) $$

We split $\tilde{Q}_{conv}$ into two terms which will both include a compensation for the other linear terms, linear with $T_p$ (such as conductive, or conducto-convective heat fluxes), or with $T_a$ (such as conducto-convective heat fluxes). Therefore, the convection term (see (15)) appears to be included in the difference between the two terms:

$$ A(T_p) = \lambda_1 T_p, \quad C(T_a) = \lambda_3 T_a \quad (19) $$

The accurate prediction of $Q_{h evap}$ is crucial as it has a high influence on both states $h$ and $T_p$. Therefore, the $Q_{h evap}$ term in (4), (13) and (14) is just slightly simplified in the reduced-model, due to its high complexity. Like the convective heat transfer coefficient $h_{conv}$, the parameter $K$ in (14) is also dependent on the flow regime. Here again, we choose a nominal
parameter value assuming a consistently turbulent flow and retain the wind velocity term, which originates from the Reynolds number. This simplification leads to the expression of the nominal parameter:

$$ \tilde{K} = 0.0021 v_0^{0.8}. \quad (20) $$

When we combine this term with the remaining constants $\frac{M_w}{R}$, and $L_w$ the result is

$$ \begin{aligned} \tilde{Q}_{h, evap} &= S \lambda_6 v_0^{0.8} (P(T_p) - \eta \times P(T_a)) \\ \tilde{Q}_{h, evap} &= S F(T_p, T_a, v_0, \eta) \end{aligned} \quad (21) $$

The basic structure of the term $P(\zeta)$ is a nonlinear non polynomial function of the temperatures. Combining (21) and (13), the function $F(T_p, T_a, v_0, \eta)$ is then given by:

$$ F(T_p, T_a, v_0, \eta) = \lambda_6 v_0^{0.8} (P(T_p) - \eta \times P(T_a)) \quad (22) $$

This function is then used to compute the evaporation term (13)-(14).

Conversely, numerical analyses reveal that the impacts of the fluxes $Q_{conv,cond,w,co}$ and $Q_{cond}$ on the predicted temperature are minimal. These terms will be indeed hidden behind other terms (mainly constants and linear with respect to $T_a$ and $T_p$) that will be automatically identified during the parameter estimation step. It results that equation (6) of the full-model does not appear explicitly but is embedded in the other model terms. The influent heat term $Q_{h,in}$ (11) is a simple term that is valid for any reactor configuration and depends on inflow rate $q^{in}$ and the reactor temperature $T_p$. Here, we make the assumption that the inflow temperature $T_{in}$ matches the ambient temperature $T_a$. Additionally, when computing the time derivative of $VT_p$ in (4), the evaporation term $F(T_p, T_a, v_0, \eta)$ related to dynamics of $\dot{V}$ appears (see (B.4) in Appendix B). This yields the following equation

$$ \begin{aligned} Q_{flow}(T_p, q^{in}) &= (q^{in} + q^{rain})(T_a - T_p) - \\ & \qquad T_p \theta_6 S q_{flow, ev} \\ q_{flow, ev} &= \lambda_7 F(T_p, T_a, v_0, \eta) \end{aligned} \quad (23) $$

where $q^{rain} = q^r \cdot S$ (m$^3$ s$^{-1}$) is the total rain volume arriving at the raceway pond. The inflow rate $q^{in}$ and the outflow are $q^{out}$ drive the process dynamics and are often the main control variables of the process. Similarly, the liquid depth can also be an efficient control variable to regulate temperature. Therefore, the term $Q_{flow}$ and two inputs $q^{in}$ and $q^{out}$ are maintained in the versatile reduced-model since they will be used as control variables of the system. The resulting model is then given by

$$ \dot{T}_p = (\theta_1 A(T_p) + \theta_2 B(T_p^4) + \theta_3 C(T_a) + \theta_4 D(T_a^4) \\ + \theta_5 E(I_0) + \theta_6 F(T_p, T_a, v_0, \eta))/h + Q_{flow}(T_p, q^{in})/(hS) \quad (24a) $$

$$ \dot{h} = \theta_6 \lambda_7 F(T_p, T_a, v_0, \eta) + q^{in}/S - q^{out}/S + q^r \quad (24b) $$

where the quantities $\lambda_i$ involved in the computation of the terms $A(.)$ to $F(.)$, together with the surface of the reactor $S$ are the known nominal values.

As evidenced in (24), the terms $A$ to $F$ are each rescaled by a respective parameter which are adapted to each case study. By adjusting these parameters, the model can better align with the actual system behaviour, compensating for the simplifications made and enhancing the accuracy of predictions. Model (24) is called the Simplified Auto Tuning Heat Exchange (SATHE) model.

The vector of parameters involved in the SATHE model is therefore defined as follows:

$$ \theta = [\theta_1, \theta_2, \theta_3, \theta_4, \theta_5, \theta_6] \quad (25) $$

<table><caption>Table 3: List of all values for the constants $\lambda_i$ involved in the SATHE model.</caption><thead><tr><th>λ<sub>i</sub></th><th>value</th><th>λ<sub>i</sub></th><th>value</th></tr></thead><tbody><tr><td>λ<sub>1</sub></td><td>-2.009e-06 (m s<sup>-1</sup>)</td><td>λ<sub>7</sub></td><td>0.002 (K<sup>-1</sup>)</td></tr><tr><td>λ<sub>2</sub></td><td>-1.318e-14 (m s<sup>-1</sup> K<sup>-3</sup>)</td><td>λ<sub>8</sub></td><td>1.035 (-)</td></tr><tr><td>λ<sub>3</sub></td><td>2.009e-06 (m s<sup>-1</sup>)</td><td>λ<sub>9</sub></td><td>0.976 (-)</td></tr><tr><td>λ<sub>4</sub></td><td>1.055e-14 (m s<sup>-1</sup> K<sup>-3</sup>)</td><td>λ<sub>10</sub></td><td>-230.543 (-)</td></tr><tr><td>λ<sub>5</sub></td><td>2.337e-07 (m<sup>3</sup> K J<sup>-1</sup>)</td><td></td><td></td></tr><tr><td>λ<sub>6</sub></td><td>-2.620e-06 (K<sup>2</sup> Pa<sup>-1</sup> m<sup>0.2</sup> s<sup>-0.2</sup>)</td><td></td><td></td></tr></tbody></table>

The number of parameters to be determined reduces significantly to 6 in the case of the generic adaptive SATHE model, in contrast to the 27 parameters required by the full-model [16] while considering the same surface area $S$. Additionally, the SATHE model has 2 states as opposed to the 3 states in the full-model [11].

The SATHE model employs a dynamic parameter adjustment mechanism, allowing it to flexibly respond to changing environmental conditions and system configurations. In the subsequent section, we present a rigorous identifiability proof, demonstrating that a single parameter value can provide a given temperature profile. Additionally, we define the parameter identification process, illuminating how our adaptive model refines parameters for optimal performance.

# 4. Model identifiability and parameter identification

## 4.1. Identifiability

Identifying the parameters of a mathematical model from experimental data is a fundamental issue in many scientific fields. However, not all models can have their parameters uniquely identified from data. In that case, it means that different parameter values can lead to the same predictions, making it impossible to infer the physical values of the parameters. This is known as the identifiability problem, and it is particularly relevant in the context of complex systems, where models may have many correlated parameters and non-linear dynamics [17, 18]. The parameter vector $\theta$ is said to be identifiable from the output $y$ if the following property is satisfied:

$$ \begin{aligned} \dot{x}_1 &= f(x_1, u, \theta_1) \text{ and } y_1 = h(x_1, u, \theta_1) \\ \dot{x}_2 &= f(x_2, u, \theta_2) \text{ and } y_2 = h(x_2, u, \theta_2) \\ \text{Then } y_1 = y_2 &\implies \theta_1 = \theta_2, \end{aligned} \quad (26) $$
where $x_1$ and $x_2$ are two state vectors associated to two different parameters $\theta$, $u$ the inputs and $y$ the output.

Structural identifiability is a fundamental concept in system modelling, assessing whether the parameters of a mathematical model can be uniquely determined from input-output data, considering an idealised, noise-free scenario. It investigates into the underlying mathematical structure, examining whether each parameter's influence on the system's behaviour is distinguishable. This theoretical analysis is crucial, as it establishes whether the model possesses non-redundant parameter estimates, essential for making accurate predictions based on real-world data [17, 18].

In our specific scenario, model (24) is linear with respect to its parameters. Previous research [19] has presented specific results on proving identifiability for this class of systems when the full state (here $T_p$ and $h$) is assumed to be measured. These results can be extended for time-varying systems that are linear-in-parameters with combinations of the functions of the state and input variables.

The primary objective in our identifiability problem is to determine if there exists a unique solution to the equations

$$ \begin{pmatrix} \dot{T}_p(t_1) \\ \dot{T}_p(t_2) \\ \vdots \\ \dot{T}_p(t_N) \end{pmatrix} = \mathcal{A}(t_i)\theta + \begin{pmatrix} \mathcal{B}(t_1) \\ \mathcal{B}(t_2) \\ \vdots \\ \mathcal{B}(t_N) \end{pmatrix} \quad (27) $$

where

$$ \mathcal{A} = \frac{1}{h} \begin{pmatrix} A(t_1) & B(t_1) & C(t_1) & D(t_1) & E(t_1) & F(t_1) \\ A(t_2) & B(t_2) & C(t_2) & D(t_2) & E(t_2) & F(t_2) \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ A(t_N) & B(t_N) & C(t_N) & D(t_N) & E(t_N) & F(t_N) \end{pmatrix} \quad (28) $$

and $\mathcal{B} = Q_{flow}$, given by $N$ measurements relating temperature $T_p$ to its derivative $\dot{T}_p$.

**Property 1.** There does not exist any linear relationship between $A(t)$, $B(t)$, $C(t)$, $D(t)$, $E(t)$ and $F(t)$ so that, for $N > 6$, matrix $\mathcal{A}$ is of full rank.

Property 1 means that (27) can be inverted providing a unique parameter $\theta$.

**Proof:** Matrix $\mathcal{A}$ is composed of functions involving terms including $T_p$, $T_p^4$, $1/T_p$, and independent inputs. Let us assume that there is a linear relationship between these terms:

$$ \gamma_1 A(t) + \gamma_2 B(t) + \gamma_3 C(t) + \gamma_4 D(t) + \gamma_5 E(t) + F(t) = 0 \quad (29) $$

Differentiating this expression provides the following equation:

$$ \dot{T}_p = \frac{T_a(-\gamma_3 C_0 - 4\gamma_4 D_0 T_p^3 + F_0 \eta P'(T_a)) + I_0 \gamma_5 E_0}{\gamma_1 A_0 + 4\gamma_2 T_p^3 B_0 + F_0 P'(T_p)} \quad (30) $$

Now, equalising (30) and (24a) provides a constraint that the independent input $T_a$ and $I_0$ should satisfy, which is not possible. This demonstrates that $\mathcal{A}$ is of full rank. ■

Even if a model is structurally identifiable, it may not be practically identifiable, meaning that the available data may not be sufficient to practically estimate the model parameters. Therefore, it is important to prove practical identifiability to ensure that the experimental data collected can be used to accurately estimate the model parameters [17, 18].

Practical identifiability analysis involves investigating the sensitivity of the model output to changes in the parameters and assessing the correlation between the parameters. The Fisher Information Matrix (FIM) is commonly used to quantify the amount of information contained in the data about the parameters and assess the precision of parameter estimation. The FIM is defined as

$$ FIM = \sum_{i=1}^{N} \left[ \frac{\partial y}{\partial \theta}(t_i) \right]^T K \left[ \frac{\partial y}{\partial \theta}(t_i) \right] \quad (31) $$

where, $\frac{\partial y}{\partial \theta} = (\sigma_{T_p}, \sigma_h)^T$ represents the sensitivity functions, and $K$ is the inverse of the measurement error covariance matrix [17]. By showing that the FIM has full rank, we prove the practical identifiability of the model.

Closer information and the full identifiability proofs can be found in Appendix A.1 and in [11], where Gharib et al. have already investigated the practical model identifiability.

## 4.2. Parameter Identification

In this section, we present the parameter identification problem written as an optimisation problem [11]. In practice, it is rare to measure pond depth. To make the approach easily applicable to many case studies, we base the parameter identification only on the temperature measurements which are always implemented. The problem is stated as

$$ \min_{\theta} \quad \bar{E}(\theta) = \sum_{i=1}^{N_H} (\delta(t_i, \theta))^2, \quad (32a) $$

$$ \text{s.t.} \quad \mathbf{\Theta} = \theta : \theta_k \ge 0, \quad \text{for } k \in 1, \dots, 6 \quad (32b) $$

$$ \delta(t_i, \theta) = T_p(t_i, \theta) - T_{p,m}(t_i), \quad \text{for } i \in 1, \dots, N_H, \quad (32c) $$

where prediction errors $\delta$ are calculated as the difference between the predicted temperature $T_p(t_i, \theta)$ and the measured temperature $T_{p,m}$ at each time instance $t_i$. The objective is to determine the parameter values, denoted as $\theta_k$, that minimise temperature prediction errors, subject to the parameter constraints in (32b). The initial parameter-set is $\theta_{init} = [1, 1, 1, 1, 1]$.

To solve this optimisation problem, we use the Levenberg-Marquardt algorithm (lsqnonlin function in Matlab). Analytical gradients are provided in (A.1). Leveraging these analytical gradients helps achieve efficient parameter optimisation and improve prediction accuracy.

The parameter uncertainty is derived first, from the sensitivity system (see (A.1)) from which we obtain the FIM, (see (31)). Then we compute the inverse of the FIM:

$$ V = FIM^{-1}. \quad (33) $$
<figure><img src="image_3.png" /><figcaption>Figure 2: Temperature and depth prediction for the reactor in Narbonne: using 4 days for parameter optimisation in summer. Experimental data (blue line), SATHE model calibration results (black dashed line), and SATHE model validation results (red line).</figcaption></figure>

<figure><img src="image_4.png" /></figure>

We use the diagonal entries $V_{ii}$ to calculate the parameter's confidence intervals

$$ \theta_i \pm t_{\alpha;N-p} \sqrt{V_{ii}} \quad (34) $$

where $t_{\alpha;N-p}$ are values from Student distributions [17]. In this study, we assume a temperature measurement error with standard deviation of 0.1°K. With a confidence level of $\alpha = 90\%$ we arrive to $t_{\alpha;N-p} = 1.645$. For propagating the parameter uncertainty on state prediction we use cooperativity properties [20, 21]. Assuming a confidence interval for each parameter $\theta_i \in [\theta_i; \bar{\theta}_i]$ derived from the previous step, we compute the resulting confidence interval $[\underline{T}_p; \bar{T}_p]$ for the state $T_p$ as follows from (27):

$$ \underline{T}_p = (\bar{\theta}_1 A(\underline{T}_p) + \bar{\theta}_2 B(\underline{T}_p^4) + \theta_3 C(T_a) + \theta_4 D(T_a^4) + \theta_5 E(I_0) + \bar{\theta}_6 F(\underline{T}_p, T_a, v_0, \eta))/\bar{h} + Q_{flow}(\underline{T}_p, q^{in})/(\bar{h}S) \quad (35a) $$

$$ \bar{T}_p = (\bar{\theta}_1 A(\bar{T}_p) + \bar{\theta}_2 B(\bar{T}_p^4) + \bar{\theta}_3 C(T_a) + \bar{\theta}_4 D(T_a^4) + \bar{\theta}_5 E(I_0) + \bar{\theta}_6 F(\bar{T}_p, T_a, v_0, \eta))/\underline{h} + Q_{flow}(\bar{T}_p, q^{in})/(\underline{h}S) \quad (35b) $$

Notably, this way of propagating the uncertainties on the temperature is a conservative way, since we conjugate for worst-case scenarios.

Expanding on the adaptive model and its identification capabilities, we shift our focus to validating the model predictions across two distinct reactors and various seasons.

## 5. Raceway in Narbonne

This section presents the results of applying the versatile model to predict the temperature dynamics in an outdoor reactor located in Narbonne, in the south of France. For quantifying the error along the temperature prediction of the SATHE model $T_{p,RM}$ and the full-model $T_{p,FM}$, we calculate the Root Mean Squared Error (RMSE)

$$ \text{RMSE}(\theta) = \sqrt{\frac{\sum_{i=0}^{N_H-1} (\delta(t_i, \theta))^2}{N_H}} \quad (36) $$

and the Mean Absolute Error (MAE)

$$ \bar{\delta}(\theta) = \frac{1}{N_H} \sum_{i=1}^{N_H} |\delta(t_i, \theta)| \quad (37) $$

over 1 day, 4 days, 14 days and till the end of the respective plots. All results are listed in Tab. 4. The parameters of the SATHE model are obtained by applying the parameter identification procedure, as detailed in Section 4.2. The full-model, which was used as a basis for the SATHE model, is specifically configured for this reactor [5].

To fit the SATHE model to measured data, a dataset consisting of 200 points (equivalent to around 4 days) is used. The obtained results (Fig. 2) show that the SATHE model provides accurate predictions on short-and-mid time-horizons. It is worth noting that even though the depth is not considered in the optimisation objective (32), the SATHE model provides an accurate prediction of this variable, as illustrated in Fig. 2, where the maximum of the absolute error is 0.0192 m for a prediction of 26 days. However, the SATHE model is less accurate over longer periods due to the changing seasons (see Fig. 3).

<figure><img src="image_5.png" /><figcaption>Figure 3: Long-term temperature prediction results for the reactor in Narbonne: using 4 days for parameter optimisation in summer. Experimental data (blue line), SATHE model calibration results (black dashed line), and the SATHE model validation results (red line).</figcaption></figure>

For instance, during the cold period (approximately between days 200 and 270), the prediction error of the SATHE model
is larger since it is far from the calibration range performed
during summer [11]. A recalibration step may be needed for a
switch in seasons.

5.1. *Model performance during Winter*

In this section, we present the outcomes of parameter opti-
misation for the SATHE model during the early winter phase
(21st of December). The results visualised in Fig. 4, show-
case the model's performance. By utilising a dataset compris-
ing 200 data points, we refine the parameters of the SATHE
model. Over 4 days, the SATHE model shows accurate temper-
ature predictions, outperforming even the precision of the full-
model. Notably, beyond day 225, the discrepancy in prediction
error between the two models becomes negligible.

<figure><img src="image_6.png" /><figcaption>Figure 4: Simulation results with 4 days dataset (days 216-220) for calibration. Comparison between temperature dynamics in the winter season of the raceway located in Narbonne: experimental data (blue line), SATHE model (red line), and full-model (grey line).</figcaption></figure>

Fig. 5 and Fig. 6 demonstrate that extending the calibra-
tion dataset to 10 or 20 days does not significantly improve
mid/long-term temperature predictions within the raceway
pond. However, for short-term predictions, smaller datasets
yield to higher uncertainty for the prediction.

## 5.2. *Model performance in-between seasons*

In the sections above we showed that the model predictions
are reliable during a season when the temperature conditions do
not change a lot over weeks. Fig. 7 presents the results when
the model parameters are calibrated during a change in seasons.
The temperature prediction in the short term is here again pre-
cise and the error is even lower than the predictions from the
full-model.

However, in the long term, it is visual, that the SATHE model is
constantly underestimating the temperature. To ensure higher
prediction accuracy we propose to perform frequent parame-
ter recalibration. The prediction results after a recalibration are
shown in Fig. 8 and Tab. 4. Notably, the prediction uncertainty
here is high, is due to our conservative uncertainty propagation
approach (see Section. 4.2). With the results from the first re-
actor, we analyse the performance of the SATHE model in the
context of a different reactor configuration.

<figure><img src="image_7.png" /><figcaption>Figure 5: Simulation results with 10 days dataset (days 210-220) for calibration. Comparison between temperature dynamics in the winter season of the raceway located in Narbonne: experimental data (blue line), SATHE model (red line), and full-model (grey line).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6: Simulation results with 20 days dataset (days 200-220) for calibration. Comparison between temperature dynamics in the winter season of the raceway located in Narbonne: experimental data (blue line), SATHE model (red line) and full-model (grey line).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7: Temperature prediction results for the reactor in Narbonne in-between season change (spring to summer): using 4 days for parameter optimisation. Experimental data (blue line), SATHE model (red line), and full-model (grey line).</figcaption></figure>
<figure><img src="image_10.png" /></figure>

Figure 8: Recalibration of the reduced model Temperature prediction results for the reactor in Narbonne in-between season change (spring to summer): using 4 days for parameter optimisation. Experimental data (blue line), SATHE model (red line)

## 6. Raceway in Milan

To further validate the performance of our model, we considered another case study and reactor configuration. The reactor is an open raceway pond located in Milan, Italy. On the contrary to the previous scenario, this reactor is raised from the ground.

<figure><img src="image_11.png" /></figure>

Figure 9: Temperature prediction results for the reactor in Milan in summer/autumn: using 4 days for parameter optimisation. Experimental data (blue line), SATHE model (red line) and full-model (grey line).

Furthermore, it is of smaller dimensions (surface area of $S_{mil} = 3.8 \text{ m}^2$). For a more detailed description of the reactor configuration, we refer to [16, 22]. To counterbalance the water loss caused by evaporation, the liquid level within the reactor was regulated by water addition. Therefore, the SATHE model was used only for predicting temperature.

Fig. 9 illustrates the temperature predictions of both the SATHE model and the full-model for the dataset from day 124 till the end of the dataset (autumn till winter) of the raceway pond in Milan. Particularly, the SATHE model delivers temperature predictions for that period. However, the temperature prediction error in days 160-170 is remarkably high. However, a readaptation of the parameters by choosing a more suitable calibration window would result in more accurate predictions in the short term, as it is shown in Fig. 10.

So far we have shown that the SATHE model is delivering comparable results concerning the temperature prediction inside the

<figure><img src="image_12.png" /></figure>

Figure 10: Temperature prediction results for the reactor in Milan in winter: using 4 days (day 160-164) for parameter optimisation. Experimental data (blue line), SATHE model calibration results (black dashed line), and the SATHE model validation results (red line).

<figure><img src="image_13.png" /></figure>

Figure 11: Temperature prediction for the reactor in Milan in-between season changes (spring to summer): using 4 days for parameter optimisation. Experimental data (blue line), SATHE model (red line) and full-model (grey line).

raceway pond. Calibrating the SATHE model to the beginning of the dataset, which is in between winter and summer, Fig. 11 shows that the SATHE model delivers more accurate predictions on short-term, but then has an error oscillating around 5°C after day 7.

## 7. Self tuning model predictions

Models like the one of Casagli [16] or Béchet [4], are difficult to handle for model based control. In particular, they are less accurate for certain seasons. Assuming that on-line temperature measurements are available, the model can use the last days of measurements to automatically re-tune its parameters so as to enhance the temperature predictions. With this strategy, heat fluxes such as conductivity with the ground, or conduction-convection through the wind, variations of influent temperature, ... can automatically be reintegrated in updated parameter values. It further leverages data from the past four days to adapt its
parameters, ensuring the accuracy of temperature forecasts for future periods.

<figure><img src="image_14.png" /><figcaption>Figure 12: Close-up for recursive online prediction: using 4 days for parameter optimisation; applying the SATHE model optimised parameters for predicting the temperature of the next day. Experimental data (blue line), SATHE model (red line) and full-model (grey line).</figcaption></figure>

In Fig. 12, the model online prediction capabilities are shown. The red line represents predictions made using the SATHE model, and we compare these to predictions generated by the full-model. The first four days are used to calibrate (only) the SATHE model's parameters. Then starting on day 4 we simulate the temperature prediction for one whole day. We assume to know the initial temperature and depth at the beginning of every day for both models. The parameters of the SATHE model are recalibrated at the end of every day using the measured temperature data of the previous 4 days. As an initial guess for the optimisation, we use either $\theta_{init} = [1, 1, 1, 1, 1, 1]$ or the optimised parameters from the previous day. This is specified for the simulation and results we present in the following. The full-model's parameters do not undergo any further optimisation process, since they are already optimised for that specific reactors. Fig. 13 provides a representation of this process for the data set from the reactor in Narbonne.

<figure><img src="image_15.png" /><figcaption>Figure 13: Recursive prediction results for the reactor in Narbonne: using 4 days for parameter optimisation; applying optimised parameters for predicting the temperature of the next day. Experimental data (blue line), SATHE model (red line) and full-model (grey line).</figcaption></figure>

The MAE of the full-model from day 4 to the end of the prediction period is $\bar{\delta}_{FM,nar} = 0.880$ °C. In comparison, the SATHE model exhibits a lower MAE of $\bar{\delta}_{RM,nar,1} = 0.613$ °C, when employing an initial guess of $\theta_0$ for each recalibration.

Further refinement in performance is observed when utilising optimised parameters as the initial guess for the subsequent re-calibration step, resulting in $\bar{\delta}_{RM,nar,2} = 0.599$ °C. The observed decrease in error from $\bar{\delta}_{RM,nar,1}$ to $\bar{\delta}_{RM,nar,2}$ indicates incremental optimisation convergence day-by-day. Evaluating the maximum absolute error is somewhat arbitrary due to the presence of numerous outliers within the measured data.

In the simulation for the reactor in Milan (see Fig. 14), we observe notable differences in temperature prediction errors between the full- and SATHE model. The full-model yields an MAE of $\bar{\delta}_{FM,mil} = 0.982$ °C from day 4 to the end of the prediction period. When employing the SATHE model with an initial guess of $\theta_0$ for recalibration, we achieved a significantly lower MAE of $\bar{\delta}_{RM,mil,1} = 0.729$ °C. Furthermore, when applying optimised parameters as the initial guess for recalibration, the SATHE model shows a lower MAE of $\bar{\delta}_{RM,mil,2} = 0.717$ °C.

<figure><img src="image_16.png" /><figcaption>Figure 14: Recursive prediction results for the reactor in Milan: using 4 days for parameter optimisation; applying optimised parameters for predicting the temperature of the next day. Experimental data (blue line), SATHE model (red line) and full-model (grey line).</figcaption></figure>

Solving an optimal control problem (OCP) using for example Pontryagin's Maximum Principle (PMP), is possible only for rather simple models for which the computations involving the adjoint state can be carried out [23]. The SATHE model, with its simplified structure, eases the application of the PMP compared to the full-model and ensures that the model is adapted to the considered time period.

## 8. Conclusion

This paper introduces the conversion of a comprehensive full-model into the self-adaptive reduced SATHE model tailored for efficient heat transfer descriptions and precise temperature predictions in microalgae cultivation systems. Through the utilisation of a minimal set of parameters, we have developed a versatile model, proficient in capturing temperature dynamics with high accuracy over short to mid-term intervals and also for different reactor configurations.
<table><caption>Table 4: Temperature error quantification results for all figures: RMSE [°C] and MAE [°C] for the SATHE model (SATHE) and the full-model (Full) for the prediction over *d* = 1 day, *d* = 4 days, *d* = 14 days and till the end (*d* = all) of the respective plot. Specifying the calibration season (cali. season) and calibration span (cali. span) and the location of the reactor.</caption><thead><tr><th rowspan="2">Figure</th><th rowspan="2">Location</th><th rowspan="2">Cali. season</th><th rowspan="2">Cali. span</th><th rowspan="2">Model</th><th colspan="4">RMSE over a prediction span of d days</th><th colspan="4">MAE over a prediction span of d days</th></tr><tr><th>d = 1</th><th>d = 4</th><th>d = 14</th><th>d = all</th><th>d = 1</th><th>d = 4</th><th>d = 14</th><th>d = all</th></tr></thead><tbody><tr><td rowspan="2">Fig. 2</td><td>Narbonne</td><td>summer</td><td>4</td><td>SATHE</td><td>0.633</td><td>0.648</td><td>0.808</td><td>0.880</td><td>0.501</td><td>0.509</td><td>0.665</td><td>0.704</td></tr><tr><td>Narbonne</td><td>summer</td><td>4</td><td>Full</td><td>1.044</td><td>1.586</td><td>1.473</td><td>1.538</td><td>1.456</td><td>1.516</td><td>1.174</td><td>1.228</td></tr><tr><td rowspan="2">Fig. 3</td><td>Narbonne</td><td>summer</td><td>4</td><td>SATHE</td><td>0.633</td><td>0.648</td><td>0.808</td><td>4.671</td><td>0.501</td><td>0.509</td><td>0.665</td><td>3.779</td></tr><tr><td>Narbonne</td><td>summer</td><td>4</td><td>Full</td><td>1.044</td><td>1.586</td><td>1.473</td><td>1.823</td><td>1.456</td><td>1.516</td><td>1.174</td><td>1.449</td></tr><tr><td rowspan="2">Fig. 4</td><td>Narbonne</td><td>winter</td><td>4</td><td>SATHE</td><td>0.210</td><td>0.711</td><td>1.182</td><td>1.179</td><td>0.515</td><td>0.490</td><td>1.500</td><td>1.580</td></tr><tr><td>Narbonne</td><td>winter</td><td>4</td><td>Full</td><td>1.019</td><td>1.503</td><td>1.623</td><td>1.570</td><td>0.948</td><td>1.333</td><td>1.324</td><td>1.250</td></tr><tr><td rowspan="2">Fig. 5</td><td>Narbonne</td><td>winter</td><td>10</td><td>SATHE</td><td>0.390</td><td>1.208</td><td>1.990</td><td>2.147</td><td>0.300</td><td>0.957</td><td>1.541</td><td>1.683</td></tr><tr><td>Narbonne</td><td>winter</td><td>20</td><td>SATHE</td><td>0.552</td><td>0.601</td><td>1.916</td><td>1.979</td><td>0.515</td><td>0.490</td><td>1.500</td><td>1.580</td></tr><tr><td rowspan="2">Fig. 6</td><td>Narbonne</td><td>spring</td><td>4</td><td>SATHE</td><td>0.456</td><td>0.924</td><td>1.008</td><td>1.580</td><td>0.362</td><td>0.800</td><td>0.894</td><td>1.327</td></tr><tr><td>Narbonne</td><td>spring</td><td>4</td><td>Full</td><td>0.955</td><td>1.440</td><td>1.583</td><td>1.398</td><td>0.828</td><td>1.251</td><td>1.375</td><td>1.167</td></tr><tr><td rowspan="2">Fig. 8</td><td>Narbonne</td><td>spring</td><td>4</td><td>SATHE</td><td>0.799</td><td>0.975</td><td>0.858</td><td>0.837</td><td>0.602</td><td>0.846</td><td>0.707</td><td>0.699</td></tr><tr><td>Narbonne</td><td>spring</td><td>4</td><td>Full</td><td>1.061</td><td>1.193</td><td>1.280</td><td>1.278</td><td>0.827</td><td>0.946</td><td>1.071</td><td>1.068</td></tr><tr><td rowspan="2">Fig. 9</td><td>Milan</td><td>autumn</td><td>4</td><td>SATHE</td><td>0.741</td><td>1.254</td><td>1.207</td><td>1.978</td><td>0.630</td><td>0.814</td><td>0.852</td><td>1.271</td></tr><tr><td>Milan</td><td>autumn</td><td>4</td><td>Full</td><td>0.526</td><td>1.102</td><td>0.988</td><td>2.248</td><td>0.459</td><td>0.884</td><td>0.791</td><td>1.377</td></tr><tr><td rowspan="2">Fig. 10</td><td>Milan</td><td>winter</td><td>4</td><td>SATHE</td><td>1.263</td><td>1.737</td><td></td><td></td><td>0.960</td><td>1.421</td><td></td><td></td></tr><tr><td>Milan</td><td>winter</td><td>4</td><td>Full</td><td>1.862</td><td>4.644</td><td></td><td></td><td>1.770</td><td>4.112</td><td></td><td></td></tr><tr><td rowspan="2">Fig. 11</td><td>Milan</td><td>spring</td><td>4</td><td>SATHE</td><td>0.675</td><td>1.784</td><td>4.235</td><td>6.243</td><td>0.462</td><td>1.148</td><td>3.646</td><td>5.797</td></tr><tr><td>Milan</td><td>spring</td><td>4</td><td>Full</td><td>2.513</td><td>3.012</td><td>2.021</td><td>1.433</td><td>2.205</td><td>2.703</td><td>1.575</td><td>1.090</td></tr><tr><td rowspan="2">Fig. 13</td><td>Narbonne</td><td>recursive</td><td>4</td><td>SATHE</td><td></td><td></td><td></td><td>0.838</td><td></td><td></td><td></td><td>0.599</td></tr><tr><td>Narbonne</td><td>recursive</td><td>4</td><td>Full</td><td></td><td></td><td></td><td>1.241</td><td></td><td></td><td></td><td>0.880</td></tr><tr><td rowspan="2">Fig. 14</td><td>Milan</td><td>recursive</td><td>4</td><td>SATHE</td><td></td><td></td><td></td><td>1.075</td><td></td><td></td><td></td><td>0.717</td></tr><tr><td>Milan</td><td>recursive</td><td>4</td><td>Full</td><td></td><td></td><td></td><td>1.460</td><td></td><td></td><td></td><td>0.982</td></tr></tbody></table>

We examined the model's performance for various seasons and compared it to a more complex model (full-model). Recalibrating the model's parameters for the different seasons consistently yields higher accuracy compared to complex models for temperature predictions. Notably, complex models such as those proposed by Casagli [16] and Béchet [4] have the advantage of not needing any recalibration, but they exhibit reduced accuracy, particularly during winter. The resulting SATHE model stands out in delivering accurate predictions following parameter recalibration.

We also demonstrate that the SATHE model outperforms the full-model in online applications. This aspect is particularly significant in showing the practical utility of the SATHE model for real-time model-based control. The improved prediction accuracy achieved by the SATHE model translates directly into better control outcomes. The introduction of the self-adaptive SATHE model represents a significant breakthrough in the field, providing a powerful tool for accurately forecasting temperature across different reactor setups.

In conclusion, this paper provides a valuable contribution to the field of microalgae cultivation by presenting a flexible model capable of accurately predicting reactor temperature.

## Acknowledgements

This work has received funding from the Digitalgaesation project within the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No. 955520.

## References

[1] E. Jacob-Lopes, M. M. Maroneze, M. I. Queiroz, L. Q. Zepka, Handbook of microalgae-based processes and products: Fundamentals and advances in energy, food, feed, fertilizer, and bioactive compounds, Academic Press, 2020.

[2] M. Kumar, Y. Sun, R. Rathour, A. Pandey, I. S. Thakur, D. C. Tsang, Algae as potential feedstock for the production of biofuels and value-added products: Opportunities and challenges, Science of the Total Environment 716 (2020) 137116, publisher: Elsevier.

[3] M. Ras, J.-P. Steyer, O. Bernard, Temperature effect on microalgae: a crucial factor for outdoor production, Reviews in Environmental Science and Bio/Technology 12 (2) (2013) 153–164. doi:10.1007/s11157-013-9310-6.

[4] Q. Béchet, A. Shilton, J. B. Park, R. J. Craggs, B. Guieysse, Universal temperature model for shallow algal ponds provides improved accuracy, Environmental science & technology 45 (8) (2011) 3702–3709, number: 8 Publisher: ACS Publications.

[5] F. Casagli, O. Bernard, Simulating biotechnological processes affected by meteorology: Application to algae–bacteria systems, Journal of Cleaner Production 377 (2022) 134190, publisher: Elsevier.

[6] E. Rodríguez-Miranda, F. G. Acién, J. L. Guzmán, M. Berenguel, A. Visoli, A new model to analyze the temperature effect on the microalgae performance at large scale raceway reactors, Biotechnology and Bioengineering 118 (2) (2021) 877–889.

[7] J. Pruvost, V. Goetz, A. Artu, P. Das, H. Al Jabri, Thermal modeling and optimization of microalgal biomass production in the harsh desert conditions of State of Qatar, Algal Research 38 (2019) 101381. doi:10.1016/j.algal.2018.12.006.

[8] A. Alonso-Bastida, M. Franco-Nava, M. Adam-Medina, B. López-Zapata, P. Álvarez Gutiérrez, J. García-Morales, Mathematical modeling of thermal interactions in a self-cooling pilot-scale photobioreactor, Case Studies in Thermal Engineering 31 (2022) 101825, publisher: Elsevier.

[9] F. Casagli, S. Rossi, J. P. Steyer, O. Bernard, E. Ficara, Balancing microalgae and nitrifiers for wastewater treatment: can inorganic carbon limitation cause an environmental threat?, Environmental science & technology 55 (6) (2021) 3940–3955, number: 6 Publisher: ACS Publications.

[10] V. Goetz, F. Le Borgne, J. Pruvost, G. Plantard, J. Legrand, A generic temperature model for solar photobioreactors, Chemical Engineering Journal 175 (2011) 443–449. doi:10.1016/j.cej.2011.09.052.
[11] A. Gharib, W. Djema, F. Casagli, O. Bernard, An auto-adaptive heat transfer model predicting the temperature evolution in photobioreactors, IFAC-PapersOnLine 56 (2) (2023) 9715-9720. doi:10.1016/j.ifacol.2023.10.284.

[12] O. Bernard, B. Rémond, Validation of a simple model accounting for light and temperature effect on microalgal growth, Bioresource Technology 123 (2012) 520-527. doi:10.1016/j.biortech.2012.07.022.

[13] E. F. Camacho, C. Bordons, Model Predictive control, Advanced Textbooks in Control and Signal Processing, Springer London, London, 2007. doi:10.1007/978-0-85729-398-5.

[14] J. Huisman, F. J. Weissing, Light-Limited Growth and Competition for Light in Well-Mixed Aquatic Environments: An Elementary Model, Ecology 75 (2) (1994) 507-520. doi:10.2307/1939554.

[15] Q. Béchet, A. Shilton, B. Guieysse, Modeling the effects of light and temperature on algae growth: state of the art and critical assessment for productivity prediction during outdoor cultivation, Biotechnology advances 31 (8) (2013) 1648-1663, number: 8 Publisher: Elsevier.

[16] F. Casagli, O. Bernard, Coupling heat transfer modelling to ALBA model for full predictions from meteorology, IFAC-PapersOnLine 55 (20) (2022) 558-563, number: 20 Publisher: Elsevier.

[17] D. Dochain (Ed.), Bioprocess control, Control systems, robotics and manufacturing series, Wiley [u.a.], Hoboken, NJ, 2008.

[18] E. Walter, L. Pronzato, Identification of parametric models from experimental data, Communications and control engineering, Springer; Masson, Berlin; New York: Paris, 1997.

[19] S. Stanhope, J. E. Rubin, D. Swigon, Identifiability of linear and linear-in-parameters dynamical systems from a single trajectory, SIAM Journal on Applied Dynamical Systems 13 (4) (2014) 1792-1815, number: 4 Publisher: SIAM.

[20] H. L. Smith, Monotone dynamical systems: an introduction to the theory of competitive and cooperative systems: an introduction to the theory of competitive and cooperative systems, no. 41, American Mathematical Soc., 1995.

[21] B. Chachuat, O. Bernard, Probabilistic observers for a class of uncertain biological processes, International Journal of Robust and Nonlinear Control: IFAC-Affiliated Journal 16 (3) (2006) 157-171.

[22] A. Pizzera, D. Scaglione, M. Bellucci, F. Marazzi, V. Mezzanotte, K. Parati, E. Ficara, Digestate treatment with algae-bacteria consortia: A field pilot-scale experiment in a sub-optimal climate area, Bioresource technology 274 (2019) 232-243, publisher: Elsevier.

[23] J.-B. Caillau, R. Ferretti, E. Trélat, H. Zidani, An algorithmic guide for finite-dimensional optimal control problems, in: Handbook of Numerical Analysis, Vol. 24, Elsevier, 2023, pp. 559-626. doi:10.1016/bs.hna.2022.11.006.

## Appendix A. Identifiability

### Appendix A.1. Practical identifiability

The sensitivity functions $\frac{\partial y}{\partial \theta} = (\sigma_{T_p}, \sigma_h)^T$ in (31) can be calculated by solving the ODE

$$ \frac{d}{dt} \left( \frac{\partial y}{\partial \theta} \right) = \begin{bmatrix} \dot{\sigma}_{T_p} \\ \dot{\sigma}_h \end{bmatrix} = \begin{bmatrix} \frac{\partial \dot{T}_p}{\partial T_p} & \frac{\partial \dot{T}_p}{\partial h} \\ \frac{\partial \dot{h}}{\partial T_p} & \frac{\partial \dot{h}}{\partial h} \end{bmatrix} \begin{bmatrix} \sigma_{T_p} \\ \sigma_h \end{bmatrix} + \begin{bmatrix} \frac{\partial \dot{T}_p}{\partial \theta} \\ \frac{\partial \dot{h}}{\partial \theta} \end{bmatrix} \quad (A.1) $$

where the Jacobian matrix entries are calculated by

$$ \frac{\partial \dot{T}_p}{\partial T_p} = \frac{1}{h} \left( \theta_1 \lambda_1 + 4\theta_2 \lambda_2 T_p^3 - \frac{q^{in} + q^{rain}}{S} + \theta_6 \lambda_6 \left( 1 - \lambda_7 T_p \right) \left( -\frac{P(T_p)}{T_p} + \frac{0.5 \lambda_9 P(T_p)}{(\lambda_{10} + T_p)^{0.5}} \right) v_0^{0.8} - \lambda_7 \theta_6 F(T_p, T_a, v_0, \eta) \right) \quad (A.2a) $$

$$ \frac{\partial \dot{T}_p}{\partial h} = -\frac{1}{h^2} \left( \theta_1 A(T_p) + \theta_2 B(T_p) + \theta_3 C(T_a) + \theta_4 D(T_a) + \theta_5 E(I_0) + \theta_6 F(T_p, T_a, v_0, \eta) + \frac{Q_{flow}(T_p, q^{in})}{S} \right) \quad (A.2b) $$

$$ \frac{\partial \dot{h}}{\partial T_p} = \theta_6 \lambda_6 \lambda_7 \left( -\frac{P(T_p)}{T_p} + \frac{0.5 \lambda_9 P(T_p)}{(\lambda_{10} + T_p)^{0.5}} \right) v_0^{0.8} \quad (A.2c) $$

$$ \frac{\partial \dot{h}}{\partial h} = 0 \quad (A.2d) $$

and the partial derivatives of the states regarding the parameters are:

$$ \left[ \begin{array}{c} \frac{\partial \dot{T}_p}{\partial \theta} \\ \frac{\partial \dot{h}}{\partial \theta} \end{array} \right] = \left[ \begin{array}{ccccc} \frac{\lambda_1 T_p}{h} & \frac{\lambda_2 T_p^4}{h} & \frac{\lambda_3 T_a}{h} & \frac{\lambda_4 T_a^4}{h} & \frac{\lambda_5 I_0}{h} \\ 0 & 0 & 0 & 0 & 0 \end{array} \right] \frac{(1 - \lambda_7 T_p) F(T_p, T_a, v_0, \eta)}{\lambda_7 F(T_p, T_a, v_0, \eta)} \quad (A.3) $$

## Appendix B. Derivating the heat transfer from the Energy balance

The change of energy inside a system is defined by

$$ \frac{dE}{dt} = \frac{d(E_{kin} + E_{pot} + E_{in})}{dt} = \frac{dE_{in}}{dt} \quad (B.1) $$

since we don't care about the external energies, we assume the kinetic and potential energies to be 0 ($E_{kin} = E_{pot} = 0$) and we focus only on the inner energy $E_{in}$ of our system.

$$ \begin{aligned} \frac{dE}{dt} &= \rho_w c_{p_w} \frac{d(V_{T_p})}{dt} = \rho_w c_{p_w} \left( V \frac{d(T_p)}{dt} + T_p \frac{d(V)}{dt} \right) \\ &= \rho_w c_{p_w} (h S \dot{T}_p + T_p S \dot{h}) \end{aligned} \quad (B.2) $$

with

$$ S \dot{h}(t) = q^{in} - q^{out} + S \theta_6 q_{flow, ev} + q^{rain} \quad (B.3) $$

where $q_{flow, ev}$ includes the evaporation term $F(T_p, T_a, v_0, \eta)$ (see (23). The general form of the energy balance equation is given by:

$$ \begin{aligned} \frac{dE}{dt} &= \sum(\dot{m}_{i,o} \cdot h_{i,o}) + \tilde{Q}_g - W^{\approx 0} + \dot{m}_{gen} \cdot h_{gen} \approx 0 \\ \sum(\dot{m}_{i,o} \cdot h_{i,o}) &= \rho_w c_{p_w} ((q^{in} + q^{rain}) \cdot T_a - q^{out} \cdot T_p) \\ \tilde{Q}_g &= \theta_1 A(T_p) + \theta_2 B(T_p^4) + \theta_3 C(T_a) + \theta_4 D(T_a^4) \\ &\quad + \theta_5 E(I_0) + \theta_6 F(T_p, T_a, v_0, \eta) \\ \tilde{Q}_g &= \tilde{Q}_g S \rho_w c_{p_w} \end{aligned} $$

where $\dot{m}_{i,o}$ is the inflow and outflow and $\dot{m}_{gen}$ the generated mass/volume flow rate (here also rain), $h_{i,o}$, $h_{gen}$ the specific enthalpy of the inflowing, outflowing, and generated substance, $\tilde{Q}_g$ all heat fluxes into the system and $W$ the work done on the system (e.g. paddle work of the raceway pond).

Setting (B.2) = (B.4):

$$ \begin{aligned} h S \dot{T}_p + T_p S \dot{h} &= \left( (q^{in} + q^{rain}) \cdot T_a - q^{out} \cdot T_p \right) + \frac{\tilde{Q}_g}{\rho_w c_{p_w}} \\ \dot{T}_p &= \frac{\tilde{Q}_g + \rho_w c_{p_w} ((q^{in} + q^{rain}) \cdot (T_a - T_p) - T_p S \lambda_7 F(T_p, T_a, v_0, \eta))}{\rho_w c_{p_w} h S} \\ \dot{T}_p &= \frac{\tilde{Q}_g}{h} + \frac{Q_{flow}(T_p, q^{in})}{h S} \end{aligned} \quad (B.4) $$