\section{Introduction}

Numerical design of microalgae production technologies has been for one decade a source of many interesting challenges not only in engineering but also in the area of scientific computing~\cite{Eilers1993,Masojidek2003,Yoo2012,Lamare2018}. 
The potential of these emerging photosynthetic organisms finds interests for the cosmetics, pharmaceutical fields, feed, food and - in the longer term - green chemistry and energy applications~\cite{Wijffels2010}. Outdoor production is mainly carried out in open bioreactors with a raceway shape. Algae grow while exposed to solar radiation in such circular basins, where the water is set in motion by a paddle wheel~\cite{Chiaramonti2013} which homogenizes the medium for ensuring an equidistribution of the nutrients
 and guarantees that each cell will have regularly access to light~\cite{Demory2018}. 
The algae are periodically harvested, and their concentration is maintained around an optimal value~\cite{Munoz-Tamayo2013,Posten2016}. Light penetration is strongly reduced by the algal biomass, resulting in limited algal concentrations, which stay  below 1\%~\cite{Bernard2015} of the liquid mass. Above this value, the light extinction is so high that a large fraction of the population evolves in the dark and does not grow anymore. At low biomass density, a fraction of the solar light is not used by the algae and the productivity is suboptimal. Theoretical works have determined the optimal biomass for maximising productivity \cite{Masci2010,Grognard2014}. Remark that other possibilities exist to improve the entire photoproduction process in order to optimize its potential. For instance, a non flat topography shape has been observed to improve the growth rate of the algae in a specific case in~\cite{Bernard02970776}. Here, we are going to prove this result, considering different cases. For this we develop a coupled model to describe the growth of algae in a raceway pond, accounting for the light that they receive. More precisely, extending the study of~\cite{Bernard2013}
 by combining the Han photosynthesis equations with an hydrodynamic law based on the shallow water equations. 

This approach enables us to formulate an optimization problem where the raceway topography is designed to maximize the productivity. For this problem, we present an adjoint-based optimization scheme which includes the constraints associated to the shallow water regime. On the contrary to a widespread belief, we prove that the flat topography is optimal in a periodic case for productivity in laminar regime. However, non-trivial topographies can be obtained in other contexts, e.g., when an extra mixing strategy is included in the model. Note that in the examples considered in our numerical tests such topographies only slightly improved the biomass production.


The outline of the paper is as follows. In Section 2, we present the biological and hydrodynamic models underlying our coupled model. In Section 3, we describe the optimization problem and a corresponding numerical optimization procedure. Section 4 is devoted to the numerical results obtained with our approach. We conclude in Section 5 with some perspectives opened by this work.

\section{Coupling hydrodynamic and biological models}

Our approach is based on a coupling between the hydrodynamic behavior of the particles and the evolution of the photosystems driven by the light intensity they receive when traveling across the raceway pond.





\subsection{Modeling the photosystems dynamics}

We consider the Han model~\cite{Han2001} which describes the dynamics of the reaction centers. These subunits of the photosynthetic process harvest photons and transfer their energy to the cell to fix CO$_2$.  In this compartmental model, the photosystems can be described by three different states: open and ready to harvest a photon ($A$), closed while processing the absorbed photon energy ($B$), or inhibited if several photons have been absorbed simultaneously ($C$). The relation of these three states are presented in Fig.~\ref{fig:han}.
\begin{figure}[thpb]
\begin{center}
\usetikzlibrary{shapes.geometric}
\tikzset{block/.style = {draw,rectangle,minimum height = 2.5em,minimum width = 2.5em}}
\begin{tikzpicture}
\node at (0,0)[circle,draw,scale=1.5] (A){$A$};
\node at (3,0)[circle,draw,scale=1.5] (B){$B$};
\node at (6,0)[circle,draw,scale=1.5] (C){$C$};
\node at (1.5,0) [block,draw,scale=0.8](sigI){$\sigma I$};
\node at (4.5,0) [block,draw,scale=0.8](kdsigI){$k_d\sigma I$};
\node at (1.5,1.5) [block,draw,scale=0.8](tau1){$\tau^{-1}$};
\node at (4.5,1.5) [block,draw,scale=0.8](kr){$k_r$};
\node at (0,2)(Pho1){Photon $I$};
\node at (3,2)(Pho2){Photon $I$};
\draw (A)--(sigI)
[->](sigI)--(B);
\draw (B)--(kdsigI)
[->](kdsigI)--(C);
\draw (C) edge[bend right=40] (kr)
(B) edge[<-,bend left=40] (kr);
\draw (B) edge[bend right=40] (tau1)
(A) edge[<-,bend left=40] (tau1);
\draw [red,->](Pho1)--(A);
\draw [red,->](Pho2)--(B);
\end{tikzpicture}
\end{center}
\caption{Scheme of the Han model, representing the probability of state transition, as a function of the photon flux density.} 
\label{fig:han}
\end{figure}

Their evolution satisfy the following dynamical system
\begin{equation*}
\left\{
\begin{array}{lr}
\dot{A} = -\sigma I A + \frac B{\tau},\\
\dot{B} =  \sigma I A - \frac B{\tau} + k_rC - k_d\sigma I B,\\
\dot{C} = -k_r C + k_d \sigma I B.
\end{array}
\right.
\end{equation*}
Here $A, B$ and $C$ are the relative frequencies of the three possible states with 
\begin{equation}\label{eq:abc}
A+B+C=1,
\end{equation}
and $I$ is the photon flux density, a continuous time-varying signal. The other parameters are $\sigma$, that stands for the specific photon absorption, $\tau$ which is the turnover rate, $k_r$ which represents the photosystem repair rate and $k_d$ which is the damage rate. As shown in~\cite{Lamare2018}, one can use~\eqref{eq:abc}, a fast-slow approximation and singular perturbation theory to reduce this system to a single evolution equation:
\begin{equation}\label{eq:ct}
\dot{C} = -\alpha(I) C + \beta(I),
\end{equation}
where
$\alpha(I) = \beta(I) + k_r,\text{with } \beta(I) =  k_d\tau \frac{(\sigma I)^2}{\tau \sigma I+1}.$
The net specific growth rate is obtained by balancing photosynthesis and respiration, which gives
\begin{equation}\label{eq:mu}
\mu(C,I) = -\gamma(I)C + \zeta(I),
\end{equation}
where
$\zeta(I) = \gamma(I)-R,\text{ with } \gamma(I) = \frac{k\sigma I}{\tau \sigma I+1}.$
Here, $k$ is a factor that relates received energy with growth rate. The term $R$ represents the respiration rate. 

In this framework, the dynamics of the biomass $X$ in an open system where a fraction of the biomass is continuously sampled and is derived from $\mu$:
\begin{equation}\label{eq:Xt}
\dot{X} = \mu(C,I)X-D X,
\end{equation}
where $D$ is the dilution rate. The growth rate~\eqref{eq:mu} associated with the steady state, for a constant $I$ is given by
\begin{equation}\label{eq:muss}
\mu(I) = -\gamma(I)C^*(I) + \zeta(I),
\end{equation}
where
\begin{equation}\label{eq:css}
C^*(I) =\frac{\beta(I)}{\alpha(I)}.
\end{equation}





\subsection{Steady 1D shallow water equations}\label{sec:model_hydro}

The shallow water equations are one of the most popular model for describing geophysical flows, which is derived from the free surface incompressible Navier-Stokes equations (see for instance~\cite{Gerbeau2001}). In the current study, we focus on the smooth steady state solutions of the shallow water equations in a laminar regime. Such steady states are governed by the following partial differential equations:
\begin{align}
&\partial_x(h u) = 0, \label{eq:sv1} \\ 
&\partial_x(h u^2+g\frac{h^2}2) = -g h \partial_x z_b, \label{eq:sv2}
\end{align}
where $h$ is the water elevation, $u$ is the horizontal averaged velocity of the water, the constant $g$ stands for the gravitational acceleration, and $z_b$ defines the topography. The free surface $\eta$ is given by $\eta = h + z_b$ and the averaged discharge $Q=hu$. This system is presented in Fig.~\ref{fig:shallow}.
\begin{figure}[thpb]
\begin{center}
\begin{tikzpicture}
\draw [thick](0,-4) -- (6,-4);
\draw (0,-4) node[anchor=north] {0};
\draw (6,-4) node[anchor=north] {$L$};
\draw[thick, ->] (0,-4) -- (0,1) node[anchor=east] {$z$};
\draw (0,0) node[anchor=east] {$0$};
\draw[dashed,->] (0,0) -- (6,0) node[anchor=west]  {$x$};
\draw [dashed, <-](1.5,-0.48) -- (1.5,0);
\draw (1.5,0.3) node{$\eta(x)$};
\draw[thick,red,<-] (2.5,-0.25) -- (2.5,1);
\draw[thick,red,<-] (3,0) -- (3,1);
\draw[thick,red,<-] (3.5,0.2) -- (3.5,1);
\draw[thick,red,<-] (4,0.35) -- (4,1);
\draw[thick,red,<-] (4.5,0.4) -- (4.5,1);
\draw[thick,red,<-] (5,0.3) -- (5,1);
\draw[thick,red,<-] (5.5,0.1) -- (5.5,1);
\draw (4,1.5) node {$I_s$};
\draw (0,0.2) sin (1.5,-0.5) cos (3,0) sin (4.5,0.4) cos (6,-0.2);
\draw[fill, gray] (0,-4) -- (0,-3.6) sin (1.5,-3.2) cos (3,-3.4) sin (4.5,-3.6) cos (6,-3.2) -- (6,-4);
\draw [dashed,->](5.3,0) -- (5.3,-3.5);
\draw (5.4,-1.8) node[anchor=west] {$z_b(x)$};
\draw[dashed,<-] (3.5,0.15) -- (3.5,-3.5);
\draw (3.6,-1.8) node[anchor=west]{$h(x)$};
\draw [->](1.3,-2) -- (1.6,-2) ;
\node at (2,-1.8){$u(x)$};
\draw [->](1.3,-1.5) -- (1.6,-1.5);
\draw [->](1.3,-2.5) -- (1.6,-2.5);
\draw [->](1.3,-1) -- (1.6,-1);
\draw [->](1.3,-3) -- (1.6,-3);
\end{tikzpicture}
\end{center}
\caption{Representation of the hydrodynamic model.}
\label{fig:shallow}
\end{figure}
The $z$ axis represents the vertical direction and the $x$ axis represents the horizontal direction. Besides, $I_s$ represents the light intensity at the free surface (assumed to be constant).

Integrating~\eqref{eq:sv1}, we get
\begin{equation}\label{eq:Q0}
h u=Q_0, 
\end{equation}
for a fixed positive constant $Q_0$, which implies a constant discharge in space. Then~\eqref{eq:sv2} can be rewritten as
\begin{equation}\label{eq:sv2mod}
h u\partial_x u + h\partial_x g h + h\partial_x g z_b = 0.
\end{equation}
Let us assume that $h$ is strictly positive. Dividing~\eqref{eq:sv2mod} by $h$ and using~\eqref{eq:Q0} to eliminate $u$, we get
\begin{equation*}
\partial_x \Big( \frac{Q_0^2}{2h^2} + g (h + z_b)\Big) = 0.
\end{equation*}
This equation corresponds to the Bernoulli's principle. Now let us consider two fixed constants $h(0), z_b(0)\in \R$. For all $x\in [0,L]$, we obtain
\begin{equation*}
\frac{Q_0^2}{2h^2} + g (h + z_b) = \frac{Q_0^2}{2h^2(0)} + g (h(0) + z_b(0)) =:M_0,
\end{equation*}
meaning that the topography $z_b$ satisfies
\begin{equation}\label{eq:zb}
z_b = \frac{M_0}g - \frac{Q_0^2}{2g h^2} - h.
\end{equation}
\begin{rem}\label{rem:flu}
Define the Froude number for the steady state by $Fr=u/\sqrt{g h}$. The situation $Fr<1$ corresponds to the subcritical case (i.e. the flow regime is fluvial) while $Fr>1$ is to the supercritical case (i.e. the flow regime is torrential). In particular, the threshold value of $h$ for $Fr=1$ is given by
\begin{equation*}
h_c := (\frac{Q^2_0}{g})^{\frac 13}.
\end{equation*}
\end{rem}

Because of~\eqref{eq:zb}, $h$ is the solution of a third order polynomial equation. Given a smooth topography $z_b$, if 
\begin{equation*}
h_c + z_b + \frac{Q_0^2}{2g h^2_c} - \frac{M_0}g < 0.
\end{equation*}
there exists a unique positive smooth solution of~\eqref{eq:zb} which satisfies the subcritical flow condition (see~\cite[Lemma 1]{Michel-Dansac2016}). The velocity $u$ can then be obtained thanks to~\eqref{eq:Q0}.







\subsection{Lagrangian trajectories of the algae and captured light intensity}

Let $z(t)$ be the depth of a particle at time $t$ in the raceway pond. We first determine the Lagrangian trajectory of an algal cell that starts at a given position $z(0)$ at time 0. 
        
From the incompressibility of the flow, we have $\nabla \cdot \underline{\textbf{u}} = 0$ with $\underline{\textbf{u}} = (u(x),w(x,z))$. Here,
$w(x,z)$ is the vertical velocity. This implies that
\begin{equation}\label{eq:incompress}
\partial_x u + \partial_z w = 0.
\end{equation}
Integrating~\eqref{eq:incompress} from $z_b$ to $z$ gives:
\begin{equation*}
\begin{split}
0 &= \int_{z_b}^z \big(\partial_x u(x) + \partial_z w(x, z) \big)\D z \\
&= \partial_x \int_{z_b}^z  u(x) \D z +\int_{z_b}^z \partial_z w(x, z) \D z\\
& =  \partial_x \big((z-z_b)u(x)\big) + w(x,z) - w(x, z_b) \\
& = (z - z_b) \partial_x u(x) - u(x)\partial_x z_b + w(x, z),
\end{split}
\end{equation*}
where we have used the kinematic condition at the bottom (i.e. $w(x,z_b) = u(x)\partial_x z_b$). It then follows from~\eqref{eq:zb} that
\begin{equation*}
w(x , z) = (\frac{M_0}g - \frac{3u^2(x)}{2g} - z)u'(x).
\end{equation*}
The Lagrangian trajectory is characterized by the system
\begin{equation}\label{eq:xtzt}
\begin{pmatrix}
\dot{x}(t)\\
\dot{z}(t)
\end{pmatrix}
= 
\begin{pmatrix}
u(x(t))\\
w(x(t),z(t))  
\end{pmatrix}.
\end{equation}
\begin{rem}
Note that a Brownian motion can be included in these equations to take into account the fluid viscosity. In this case, the optimisation procedure required a large set of simulations and an averaging strategy. Such questions are out of the range of this paper.
\end{rem}
Up to now  and for the sake of simplicity, we denote by $z(x)$ the depth of a particle at $x$. From~\eqref{eq:xtzt}, we get
\begin{equation}\label{eq:z'x}
z' := \frac{\dot{z}}{\dot{x}} = (\frac {M_0}g -\frac{3u^2}{2g} - z)\frac{u'}{u}.
\end{equation}
Note that from~\eqref{eq:Q0} and~\eqref{eq:zb}, we have $\eta=h+z_b = \frac{M_0}g - \frac{u^2}{2g}$, which implies $\eta'=-u u'/g$. Multiplying~\eqref{eq:z'x} both sides by $u$ and using the formulation of $\eta, \eta'$ gives
\begin{equation*}
\begin{split}
z' u + z u'&= (\eta - \frac{u^2}g)u'\\
&=\eta u' + \eta' u,
\end{split}
\end{equation*}
which implies that $\big(u(z-\eta)\big)'=0$. We then obtain


\begin{equation}\label{eq:zx}
z(x) = \eta(x) + \frac{u(0)}{u(x)}(z(0) - \eta(0)).
\end{equation}
The computation of trajectories in the shallow water system can be carried out with this formula, where $u(0)$, $z(0)$ and $\eta(0)=h(0)+z_b(0)$ are given as boundary conditions, whereas $u(x)$ can be obtained by~\eqref{eq:Q0} and~\eqref{eq:zb}. On the contrary, $h(x)$, hence $u(x)$ will be considered  in section~\ref{sec:opti} as a parameter to be optimized, and~\eqref{eq:zb} will be used to (post-)compute $z_b(x)$.  


\begin{rem}
Since $Q_0$ is chosen to be positive, $h$ is necessarily positive and so does $u$ from~\eqref{eq:Q0}. Moreover, if $z(0)$ belongs to  $[z_b(0),\eta(0)]$, then $z(x)$ belongs to $[z_b(x),\eta(x)]$. In particular, choosing $z(0)=z_b(0)$ in~\eqref{eq:zx} and using~\eqref{eq:Q0} gives $z(x) = z_b(x)$. In the same way, we find that $z(x)=\eta(x)$ when $z(0) = \eta(0)$. 
\end{rem}

To evaluate the light intensity observed on the trajectory $z$, we assume that the growth process occurs at a much slower time scale than those of hydrodynamics and is, as such, negligible for one lap over the raceway. As a consequence, the turbidity is supposed to be constant over the considered time scale. In this framework, the Beer-Lambert law describes how light is attenuated with depth:
\begin{equation}\label{eq:Beer}
I(x,z) = I_s\exp\Big(-\varepsilon (\eta(x) - z)\Big).
\end{equation}
Here $\varepsilon$ is the light extinction coefficient. We assume that the system is perfectly mixed so that the concentration of the biomass $X$ in~\eqref{eq:Xt} is homogeneous and $\varepsilon$ is constant. Combining~\eqref{eq:Beer} with~\eqref{eq:zx}, we get the following  expression for the captured light intensity along the trajectory $z$
\begin{equation}\label{eq:Ixz0}
I(x,z(x)) =I_s \exp\Big(-\varepsilon  \frac{u(0)}{u(x)} (\eta(0) -z(0) ) \Big).
\end{equation}
It follows that in this approach, the light intensity $I$ couples the hydrodynamic model and Han model: the trajectories of the algae define the received light intensity, which is used in the photosystems dynamics.

We can then derive the equation satisfied by $C$. Indeed, repeating the reasoning done to get~\eqref{eq:z'x} with \eqref{eq:ct}, we find a time-free reformulation, namely
\begin{equation}\label{eq:c'x}
C':=\frac{\dot{C}}{\dot{x}} = -\frac{\alpha(I)}{u}C +  \frac{\beta(I)}{u},
\end{equation}
where all the functions on the right-hand side only depend on $x$. 

\section{Optimization problem}\label{sec:opti}

In this section, we define the optimization problem associated with our biological-hydrodynamic model. We first introduce our procedure for a simple problem, and then extend this to other variant cases.





\subsection{Optimization functional}

The average net specific growth rate over the domain is defined by
\begin{equation}\label{eq:mucontinous}
\bar{\mu} := \frac 1L\int_0^L\frac 1{h(x)}\int_{z_b(x)}^{\eta(x)} \mu\big(C(x, z), I(x, z)\big) \D z \D x.
\end{equation}
In order to tackle numerically this optimization problem, let us consider a vertical discretization. Let $N_z$ denotes the number of particles, and consider a uniform vertical discretization of their initial position:
\begin{equation*}
z_i(0) = \eta(0) - \frac {i-\frac12}{N_z}h(0), \quad i=1,\ldots,N_z. 
\end{equation*}
From~\eqref{eq:zx}, we obtain
\begin{equation*}
z_{i}(x) - z_{i+1}(x) = \frac 1{N_z}h(x), \quad i=1,\ldots,N_z, 
\end{equation*}
meaning that the distribution of particles remains uniform along the trajectories. To simplify notations, we write $I_i(x)$ instead of $I(x,z_i(x))$ hereafter. 

Let $C_i(x)$ (resp. $I_i(x)$) the photo-inhibition state (resp. the light intensity) associated with the trajectories $z_i(x)$. Then the semi-discrete average net specific growth rate in the raceway pond can be defined from~\eqref{eq:mucontinous} by
\begin{equation}\label{eq:mudiscret}
\begin{split}
\bar \mu_{\Delta} :&= \frac 1L \int_0^L \frac 1{h(x)}\frac 1{N_z}h(x)\sum_{i=1}^{N_z}\mu( C_i(x),I_i(x) ) \D x\\
&= \frac 1{L N_z}\sum_{i=1}^{N_z} \int_0^L\mu( C_i(x),I_i(x) ) \D x.
\end{split}
\end{equation}

From now on, we focus only on the subcritical case, i.e. $Fr<1$, see Remark~\ref{rem:flu}. As we have mentioned in the previous section, in this regime, a given topography $z_b$ corresponds to a unique water height $h$ which verifies this assumption. Note that the volume of our system is simply given by
\begin{equation}\label{eq:vol}
V = \int_0^L h(x) \D x.
\end{equation}
Therefore, we choose to parameterize $h$ by a vector $a:=[a_1,\cdots,a_N]\in \R^{N}$ in order to handle the volume of our system. This vector will be the variable to be optimized. Once we find the optimum parameter $a^*$, we can determinate the associated topography by means of~\eqref{eq:zb}.
\begin{rem}
In usual shallow water solver, equations of type~\eqref{eq:zb}  are usually consider to compute $h$ in the simulations. Here, we use this equation in the opposite way, i.e., to recover the topography $z_b$ from $h$.
\end{rem}

Our goal is to optimize the topography to maximize $\bar \mu_{\Delta}$. Up to now, we denote explicitly that the functions depend on $a$.
In the case of one trajectory, we consider the functional
\begin{equation*}
J(a) = \frac1L\int_0^L \frac{-\gamma(I(x;a))C(x) + \zeta(I(x;a))}{u(x;a)} \D x ,
\end{equation*}
where $C$ satisfies the following parameterized version of~\eqref{eq:c'x}
\begin{equation}\label{eq:cond}
C' + \frac{\alpha(I(x;a))}{u(x;a)}C =  \frac{\beta(I(x;a))}{u(x;a)}.
\end{equation}
The optimization problem then reads:

\textit{Find $a^*$ solving the maximization problem:}
\begin{equation}\label{eq:Jopt}
\max_{a\in \R^{N}} J(a).
\end{equation}





\subsection{Optimality System for maximizing growth rate}


In this section, we introduce the optimality system derived from the optimization problem~\eqref{eq:Jopt} and some of its variants. For the sake of simplicity, we first remain in the setting where one single trajectory is considered, and then extend it to a multiple trajectories system. Next, we analyse the case where $C$ is periodic, and finally, we present a more general case where the volume can also be optimized. For simplicity of notation, we omit $x$ in the notation.



\subsubsection{One trajectory problem}

Define the Lagrangian of Problem~\eqref{eq:Jopt} by
\begin{equation*}
\begin{split}
\mathcal{L}(C,p,a) = \frac1L&\int_0^L \frac{-\gamma(I(a))C +\zeta(I(a))}{u(a)} \D x\\
- &\int_0^L p\big(C' + \frac{\alpha(I(a))C-\beta(I(a))}{u(a)}\big) \D x,
\end{split}
\end{equation*}
where $p$ is the Lagrange multiplier associated with the constraint~\eqref{eq:cond}.

The optimality system is obtained by cancelling all the partial derivatives of $\mathcal{L}$. Differentiating $\mathcal{L}$ with respect to $p$ and equating the resulting expression to zero gives~\eqref{eq:cond}. Integrating the terms $\int p C'\D x$ on the interval $[0,L]$ by parts and differentiating $\L$ with respect to $C$ and $C(L)$ gives rise to
\begin{equation}\label{eq:adjoint}
\left\{
\begin{array}{lr}
\partial_{C}\mathcal{L} = p' - p \frac{\alpha(I(a))}{u(a)} - \frac1L\frac{\gamma(I(a))}{u(a)} \\
\partial_{C(L)}\mathcal{L}  = p(L).\\
\end{array}
\right.
\end{equation}
Given a vector $a$, let us still denote by $C, p$ the corresponding solutions of~\eqref{eq:cond} and~\eqref{eq:adjoint}. The gradient $\nabla J(a)$ is obtained by
\begin{equation*}
\nabla J(a) = \partial_{a}\mathcal{L},
\end{equation*}
where
\begin{equation}\label{eq:grad1}
\begin{split}
\partial_{a}\mathcal{L} = \frac1L&\int_0^L \frac{-\gamma'(I(a))C+\zeta'(I(a))}{u(a)}\partial_a I(a)\D x \\
-\frac1L&\int_0^L \frac{-\gamma (I(a))C + \zeta(I(a))}{u^2(a)}\partial_a u(a)\D x \\
+&\int_0^L p\frac{-\alpha'(I(a))C + \beta'(I(a))}{u(a)}\partial_a I(a)\D x \\
-&\int_0^L p\frac{-\alpha(I(a))C + \beta(I(a))}{u^2(a)}\partial_a u(a) \D x.
\end{split}
\end{equation}





\subsubsection{Multiple trajectories problem}

We now extend the previous framework to deal with multiple trajectories. From~\eqref{eq:mudiscret}, the objective functional is defined by:
\begin{equation}\label{eq:muNza}
\bar \mu_{\Delta}(a) = \frac 1{L N_z}\sum_{i=1}^{N_z}\int_0^L\frac{-\gamma(I_i(a))C_i + \zeta(I_i(a))}{u(a)}   \D x,
\end{equation}
where $C_i$ satisfies~\eqref{eq:cond} for the light intensity $I_i$ with $i=1,\ldots,N_z$. Let us still denote by $\L$ the associated Lagrangian. Similar computations as in the previous section give rise to $N_z$ systems similar to~\eqref{eq:adjoint}, where $C$ is replaced by $C_i$. Denoting by $p_{i}$ the associated Lagrange multipliers, the partial derivative $\partial_a\mathcal{L}$ (hence the gradient $\nabla \bar \mu_{\Delta}(a)$) is then given by
\begin{equation*}
\begin{split}
\partial_{a}\mathcal{L} = \frac 1{L N_z}&\sum_{i=1}^{N_z}\int_0^L \frac{-\gamma'(I_i(a))C_i+\zeta'(I_i(a))}{u(a)}\partial_a I_i(a)\D x \\
-\frac 1{L N_z}&\sum_{i=1}^{N_z}\int_0^L \frac{-\gamma (I_i(a))C_i + \zeta(I_i(a))}{u^2(a)}\partial_a u(a)\D x \\
+&\sum_{i=1}^{N_z}\int_0^L p_{i}\frac{-\alpha'(I_i(a))C_i + \beta'(I_i(a))}{u(a)}\partial_a I_i(a)\D x \\
-&\sum_{i=1}^{N_z}\int_0^L p_{i}\frac{-\alpha(I_i(a))C_i + \beta(I_i(a))}{u^2(a)}\partial_a u(a) \D x.
\end{split}
\end{equation*}





\subsubsection{Periodic problem}

We now consider a variant of our problem, where the photo-inhibition state $C$ is periodic, with a period corresponding to a raceway lap. This situation
occurs, e.g., when an appropriate harvest is performed after each lap. To describe the corresponding model, let us first consider a variant of the usual Cauchy problem for~\eqref{eq:c'x}:


\textit{Given $I\in L^2(0,L;R)$, find $(C_0,C)\in \R\times \C(0,L;\R)$ such that}
\begin{equation}\label{eq:periodiceq}
\left\{
\begin{array}{lr}
C'(x) = \frac{-\alpha(I(x))C(x) + \beta(I(x))}{u(x)},\quad x\in[0,L]&\\
C(L) = C(0)= C_0.&
\end{array}
\right.
\end{equation}


The following theorem gives the existence and uniqueness of a (weak) solution of the problem~\eqref{eq:periodiceq}
\begin{thm}\label{thm:perC}
Given $I\in L^2(0,L;R)$, there exists a unique couple $(C_0,C) \in \R\times \C(0,L;\R)$ satisfying
\begin{equation}\label{periodicC}
\left\{
\begin{array}{lr}
C(x) = C_0 + \int_0^x\frac{-\alpha(I(s)) C(s) + \beta(I(s))}{u(s)}\D s&,\\
C(L)  = C_0&
\end{array}
\right.
\end{equation}
for all $x\in [0,L]$.
\end{thm}

The proof is given in Appendix A. The optimality system is obtained as in the previous sections. 

\begin{rem}\label{rem:pper}
Note that the periodicity of $C$ implies that $p$ is also periodic. Indeed, since $C(L)=C(0)$, then by differentiating $\mathcal{L}$ with respect to $C(L)$, we have
\begin{equation*}
\partial_{C(L)} \mathcal{L} = p(L) - p(0).
\end{equation*}
so that equating the above equation to zero gives the periodicity for $p$. Similar computations apply for multiple trajectories problem.
\end{rem}

Next, we state a result for Problem~\eqref{eq:Jopt} in the case where $C$ is periodic.

\begin{prop}\label{prop:flattopo}
Assume $h$ is of the form
\begin{equation}\label{eq:Fourier}
h(x;a) = a_0 + \sum_{n=1}^N a_n \sin(2n\pi \frac x L),
\end{equation}
with $a:=[a_1,\ldots,a_N]$ chosen such that the subcritical condition is satisfied (see Remark~\ref{rem:flu}). Then $\nabla J(0)=0$.
\end{prop}
As a consequence, the flat topography is critical point for Problem~\eqref{eq:Jopt}. Note that the larger $N$ is, the less valid is our hydrodynamic model valid, see Section~\ref{sec:model_hydro}, where a smooth topography is assumed. Hence, limit situations where $N\rightarrow+\infty$ are not considered in what follows.

\begin{pf}
Let $q \in[0,1]$, consider an initial position by
\begin{equation}\label{eq:z0}
z(0;a) =  \eta(0;a) - q h(0;a).
\end{equation}
Replacing $z(0;a)$ by the previous expression in~\eqref{eq:Ixz0} and using~\eqref{eq:Q0}, gives
\begin{equation}\label{eq:Ia}
I(a) = I_s \exp(-\varepsilon q\frac{Q_0}{u(a)}) = I_s \exp(-\varepsilon q h(a)).
\end{equation}
Introduce $a^* = [0, \ldots, 0]$. Replacing $a$ by $a^*$ in~\eqref{eq:Q0} and~\eqref{eq:Ia}, we find that $u(a^*)$ and $I(a^*)$ are both constants with respect to $x$. For a given $C(0)$,~\eqref{eq:cond} gives
\begin{equation}\label{eq:Cx}
C(x) = e^{-\frac{\alpha(I(a^*))}{u(a^*)}x}C(0) + \frac{\beta(I(a^*))}{\alpha(I(a^*))}(1 - e^{-\frac{\alpha(I(a^*))}{u(a^*)}x}).
\end{equation}
Since $C$ is periodic (i.e. $C(L)=C(0)$), we get from the previous equation that $ C(0) = \frac{\beta(I(a^*))}{\alpha(I(a^*))}$. Inserting this value in~\eqref{eq:Cx}, we find
\begin{equation*}
C(x) = \frac{\beta(I(a^*))}{\alpha(I(a^*))}, \quad \forall x\in[0,L],
\end{equation*}
which corresponds to the steady state~\eqref{eq:css}, that we denote by $C^*$ hereafter. We then compute all the partial derivatives appearing in~\eqref{eq:grad1}. We get:
\begin{equation*}
\begin{split}
&\partial_{a_n} u(a) = -\frac{Q_0}{h^2(a)}\partial_{a_n} h(a),\\
&\partial_{a_n} I(a) = -\varepsilon q I(a)\partial_{a_n} h(a),
\end{split}
\end{equation*}
but $\partial_{a_n} h(a) = \sin(2n\pi \frac x L)$, so that the integral from $0$ to $L$ of $\partial_{a_n} h(a)$ is zero. Replacing $a$ by $a^*$ in~\eqref{eq:grad1} , the first integral can be computed as
\begin{equation*}
\begin{split}
&\int_0^L \frac {-\gamma'(I(a^*))C^* + \zeta'(I(a^*))}{u(a^*)} \partial_{a_n} I(a^*) \D x\\
= &\frac {-\gamma'(I(a^*))C^* + \zeta'(I(a^*))}{u(a^*)} \int_0^L\partial_{a_n} I(a^*) \D x \\
= &0.
\end{split}
\end{equation*}
and the second integral can be computed as
\begin{equation*}
\begin{split}
 &\int_0^L\frac {-\gamma(I(a^*))C^*+\zeta(I(a^*))}{u^2(a^*)} \partial_{a_n} u(a^*)\D x \\
=&\frac {-\gamma(I(a^*))C^*+\zeta(I(a^*))}{u^2(a^*)} \int_0^L\partial_{a_n} u(a^*)\D x \\
= &0.
\end{split}
\end{equation*}
Since $p$ is also periodic (see Remark~\ref{rem:pper}), it can be proved in much the same way that the two last terms in~\eqref{eq:grad1} also cancels when $a=a^*$. The result follows. 
\end{pf}

\begin{rem}
Numerically, we observe that the flat topography is actually optimal in the periodic case. 
\end{rem}

\begin{rem}
Note that if $C$ is not periodic,~\eqref{eq:Cx} implies that $C$ depends also on the space variable $x$ and the computations for the first-two integrals in the proof above will no longer hold. For the same reason, $p$ shall be not constant and the last two integrals in~\eqref{eq:grad1}, hence the gradient will not cancel in general. In other words, the flat topography is not an optimum for the case $C$ non periodic, which is confirmed by our numerical tests  (see Section~\ref{sec:testC0}).
\end{rem}






\subsubsection{Non-constant volume problem for maximizing areal productivity}\label{sec:nonconstant}

In general, the extinction coefficient is an affine function with respect to the biomass $X$~\cite{Martinez201811}
\begin{equation}\label{eq:eps}
\varepsilon(X) = \alpha_0 X + \alpha_1,
\end{equation}
where $\alpha_0>0$ the specific light extinction coefficient of the microalgae specie and $\alpha_1$ the background turbidity that summarizes the light absorption and diffusion due to all non-microalgae components. Recall that according to~\eqref{eq:Xt}, the biomass concentration $X$ depends on the growth rate $\mu$ and on the harvesting of the raceway $D$. 

\begin{rem}\label{rem:optXcond}
A standard criterion to determine a relevant value of $X$~\cite{Masci2010,Grognard2014}, and then an harvesting frequency, consists in regulating it to a value such that the steady-state value of the growth rate $\mu$ at the (average) bottom depth $\bar z_b$ is 0, i.e., 
\begin{equation}\label{eq:muIzb}
\mu(I_{\bar z_b})=0,
\end{equation}
where $\mu(I)$ is given by~\eqref{eq:muss}. Solving~\eqref{eq:muIzb} provides a value of $I_{\bar z_b}$ thus of the extinction coefficient $\varepsilon(X)$, and finally of the biomass concentration $X$. This assumption is usually called the \textit{compensation condition}, for which the growth at the bottom compensates exactly the respiration.
\end{rem}

In this framework, maximizing areal productivity is a relevant target. Productivity per unit of surface for a given biomass concentration $X$ is given by:
\begin{equation}\label{eq:Pi}
\Pi = \bar \mu X \frac{V}{S},
\end{equation}
where $S$ is the ground surface of the raceway pond and $\bar \mu$ is defined in~\eqref{eq:mucontinous}. To optimize this objective function, a more general setting where the reactor volume varies can be considered. The total biomass $X$ in a given volume is supposed to satisfy~\eqref{eq:muIzb} and we assume that it is maintained constant to achieve the
compensation conditions, see Remark~\ref{rem:optXcond}.

Let us keep the parameterization $h$ given in~\eqref{eq:Fourier}. The volume of the system satisfies $V=a_0L$ (see~\eqref{eq:vol}). Moreover, since we consider a 1D framework, we have $S=L$. The biomass concentration $X$ is chosen to be solution of~\eqref{eq:muIzb}. We have 
\begin{equation*}
X(a_0)\frac{V(a_0)}{S}= \alpha_2 - \alpha_3 a_0,
\end{equation*}
where $\alpha_2=\frac{1}{\alpha_0}\ln(\frac{I_s}{I_{\bar z_b}})$ and $\alpha_3=\frac{\alpha_1}{\alpha_0}$, and $\alpha_0, \alpha_1$ are given in~\eqref{eq:eps}. The detail of the computation is given in Appendix C. Let us then consider the extended parameter vector $\tilde a:=[a_0,a]$. For simplicity, we again consider one single trajectory and optimize the objective function associated with~\eqref{eq:Pi} defined by  
\begin{equation*}
\tilde J(\tilde a) = \frac{\alpha_2-\alpha_3a_0}{L}\int_0^L \frac{-\gamma(I(\tilde a))C +\zeta(I(\tilde a)) }{u(\tilde a)} \D x,
\end{equation*}
where $C$ satisfies~\eqref{eq:cond} with $a$ replacing by $\tilde a$.
The corresponding optimization problem reads: 

\textit{Find $\tilde a^*$ solving the maximization problem:}
\begin{equation}\label{eq:Joptilde}
\max_{\tilde a\in \R^{N+1}} \tilde J(\tilde a).
\end{equation}
The Lagrangian of Problem~\eqref{eq:Joptilde} is 
\begin{equation*}
\begin{split}
\tilde{\mathcal{L}}(C,p,\tilde a) = &\frac{\alpha_2-\alpha_3a_0}{L}\int_0^L \frac{-\gamma(I(\tilde a))C+\zeta(I( \tilde a))}{u(\tilde a)} \D x\\
- &\int_0^L p\big(C' + \frac{\alpha(I(\tilde a))C-\beta(I(\tilde a))}{u(\tilde a)}\big) \D x,
\end{split}
\end{equation*}
where $p$ is the Lagrangian multiplier associated with the constraint~\eqref{eq:cond} for $C$. By computations similar to that of the previous section, we find as optimality system
\begin{equation*}
\left\{
\begin{array}{lr}
\partial_{C} \tilde{\mathcal{L}} = p' - p \frac{\alpha(I(\tilde a))}{u(\tilde a)} -\frac{\gamma(I(\tilde a))}{u(\tilde a)}\frac{\alpha_2-\alpha_3a_0}{L} \\
\partial_{C(L)}\tilde{\mathcal{L}}  = p(L).
\end{array}
\right.
\end{equation*}
However, there is an extra element in the gradient $\tilde J(\tilde a)$, namely the derivative with respect to $a_0$, so that $\nabla \tilde J(\tilde a): = [\partial_{a_0}\tilde{\mathcal{L}},  \partial_{a}\tilde{\mathcal{L}}]$, where
\begin{equation}\label{eq:grad1tilde1}
\begin{split}
\partial_{a_0}\tilde{\mathcal{L}} =  \frac{\alpha_2-\alpha_3a_0}{L}&\int_0^L \frac{-\gamma'(I(\tilde a))C+\zeta'(I(\tilde a))}{u(\tilde a)}\partial_{a_0} I(\tilde a)\D x\\
-\frac{\alpha_2-\alpha_3a_0}{L}&\int_0^L \frac{-\gamma (I(\tilde a))C + \zeta(I(\tilde a))}{u^2(\tilde a)}\partial_{a_0} u(\tilde a)\D x \\
-\frac{\alpha_3}{L}&\int_0^L \frac{-\gamma (I(\tilde a))C + \zeta(I(\tilde a))}{u(\tilde a)}\D x \\
+&\int_0^L p\frac{-\alpha'(I(\tilde a))C + \beta'(I(\tilde a))}{u(\tilde a)}\partial_{a_0} I(\tilde a)\D x \\
-&\int_0^L p\frac{-\alpha(I(\tilde a))C + \beta(I(\tilde a))}{u^2(\tilde a)}\partial_{a_0} u(\tilde a) \D x,
\end{split}
\end{equation}
\begin{equation}\label{eq:grad1tilde2}
\begin{split}
\partial_{a}\tilde{\mathcal{L}} =  \frac{\alpha_2-\alpha_3a_0}{L}&\int_0^L \frac{-\gamma'(I(\tilde a))C+\zeta'(I(\tilde a))}{u(\tilde a)}\partial_a I(\tilde a)\D x \\
-\frac{\alpha_2-\alpha_3a_0}{L}&\int_0^L \frac{-\gamma (I(\tilde a))C + \zeta(I(\tilde a))}{u^2(\tilde a)}\partial_a u(\tilde a)\D x \\
+&\int_0^L p\frac{-\alpha'(I(\tilde a))C + \beta'(I(\tilde a))}{u(\tilde a)}\partial_a I(\tilde a)\D x \\
-&\int_0^L p\frac{-\alpha(I(\tilde a))C + \beta(I(\tilde a))}{u^2(\tilde a)}\partial_a u(\tilde a) \D x.
\end{split}
\end{equation}
However, we still find that the flat topography is optimal in that case.
\begin{prop}
Assume $h$ is of the form~\eqref{eq:Fourier} and $C$ is periodic, then $\tilde J(\tilde a^*)=0$ with $\tilde a^*=[a^*_0,0,\cdots,0]$.
\end{prop}

\begin{pf}
The proof is similar to that of Proposition~\ref{prop:flattopo}. For a given $z(0;\tilde a)$ defined in~\eqref{eq:z0}, using the result in Appendix C to replace $\varepsilon$ in~\eqref{eq:Ia}, we have
\begin{equation*}
I(\tilde a) = I_s \exp\Big(- \frac 1{a_0}\ln(\frac{I_s}{I_{\bar z_b}})q h(\tilde a)\Big).
\end{equation*}
The two extra partial derivatives appearing in~\eqref{eq:grad1tilde1} are computed by
\begin{equation*}
\begin{split}
&\partial_{a_0} u(\tilde a) = -\frac{Q_0}{h^2(\tilde a)}\partial_{a_0} h(\tilde a),\\
&\partial_{a_0} I(\tilde a) = -( -\frac 1{a^2_0}h(\tilde a) + \frac 1{a_0}\partial_{a_0} h(\tilde a) )\ln(\frac{I_s}{I_{\bar z_b}})q I(\tilde a),
\end{split}
\end{equation*}
where $\partial_{a_0} h(\tilde a)=1$. The flat topography corresponds to the parameter $\tilde a^* = [a^*_0,0,\ldots,0]$. Then $\partial_{a_0} u(\tilde a^*)=-\frac{Q_0}{{a^*_0}^2}$ and $\partial_{a_0} I(\tilde a^*)=0$. Hence the first integral and the fourth integral in $\partial_{a_0}\tilde{\mathcal{L}}$ in~\eqref{eq:grad1tilde1} are zero. Following a reasoning similar to that of Proposition~\ref{prop:flattopo}, we find $C^*$ is constant. Hence, the last integral in $\partial_{a_0}\tilde{\mathcal{L}}$ is zero. Finally for 
\begin{equation}\label{eq:a0}
a^*_0=\frac{\alpha_2}{2\alpha_3},
\end{equation}
the sum of the second and the third integral in $\partial_{a_0}\tilde{\mathcal{L}}$ is zero. 

As for~\eqref{eq:grad1tilde2}, it follows from computations similar to that of the proof of Proposition~\ref{prop:flattopo} that this term equals zero, which completes the proof.
\end{pf}

\section{Numerical Experiments}


In this section, we show some optimal topographies obtained in the various previous frameworks.


\subsection{Numerical Methods}

We start with an algorithm to solve the optimization problem and the numerical solver that we use for our numerical tests.





\subsubsection{Gradient Algorithm}

We detail now a gradient-based optimization procedure and apply it to Problem~\eqref{eq:Jopt}. The procedure is given in Algorithm~\ref{algo}.
\begin{algorithm}
\begin{algorithmic}[1]
\State \textbf{Input}: Tol$>0$, $\rho>0$.
\State \textbf{Initial guess}: $a$.
\State \textbf{Output}: $a$
\State Set $err :=$ Tol$+1$ and define $h$ by~\eqref{eq:Fourier} using the input data.
\While{$err>$Tol and $\norm{h}_{\infty} > h_{c}$}
\State Compute $u$ by~\eqref{eq:Q0} and $I$ by~\eqref{eq:Ixz0}.
\State Set $C$ as the solution of~\eqref{eq:cond}.
\State Set $p$ as the solution of~\eqref{eq:adjoint}.
\State Compute the gradient $\nabla J$ by~\eqref{eq:grad1}.
\State $a = a + \rho \nabla J$,
\State Set $err := \| \nabla J\|$. 
\EndWhile
\end{algorithmic}
\caption{Gradient-based optimization algorithm}
\label{algo}
\end{algorithm}

Note that $\norm{h}_{\infty} := \max_{x\in[0,L]}h(x)$. In addition to a numerical tolerance criterion on the magnitude of the gradient, we need to take into account a constraint on the height $h$, to guarantee that the simulated flow remains in a subcritical regime, see Remark~\ref{rem:flu} (and in the range of industrial constraints, see~\cite{Chiaramonti2013}).


A similar algorithm can be considered to tackle the multiple trajectories Problem~\eqref{eq:muNza}. Remark that since no interaction between trajectories is considered, the gradient computation can be partially parallelized when computing $I_i$ (hence $C_i$) and $p_i$.

Finally, the optimality system presented in Section~\ref{sec:nonconstant} can still be solved with Algorithm~\ref{algo} by adapting the objective functional $\tilde J$ and the gradient $\nabla \tilde J$. 





\subsubsection{Numerical Solvers}

To solve our optimization problem numerically, we introduce a supplementary space discretization with respect to $x$. In this way, let us take a space increment $\Delta x$, set $N_x =[L/\Delta x]$ and $x^{n_x}=n_x\Delta x$ for $n_x=0,\ldots,N_x$. We choose to use the Heun's method for computing $C$ via~\eqref{eq:cond}. Following a first-discretize-then-optimize strategy, we get that the Lagrange multiplier $p$ is also computed by a
Heun's type scheme. 
Note that this scheme is still explicit since it solves a backward dynamics starting from $p(L) = 0$.





\subsection{Parameter settings}

We now detail the parameters used in our simulations.



\subsubsection{Parameterization of $h$}

As mentioned in Proposition~\ref{prop:flattopo}, we  parameterize $h$ by a truncated Fourier series~\eqref{eq:Fourier}.
The parameters to be optimized will be the Fourier coefficients $a:=[a_1,\ldots,a_N]$. Note that we choose to fix $a_0$, since it is related to the volume $V$ of the raceway with $V=a_0L$






\subsubsection{Parameter for the models}

The spatial increment is set to $\Delta x=\SI{0.01}{m}$ so that the convergence of the numerical scheme has been ensured, and we set the raceway length $L=\SI{100}{m}$, the averaged discharge $Q_0 = \SI{0.04}{m^2.s^{-1}}$, the initial value $a_0(=h(0;a))=\SI{0.4}{m}$ and $z_b(0)=-\SI{0.4}{m}$ to stay in standard ranges for a raceway. The free-fall acceleration $g=\SI{9.81}{m.s^{-2}}$. All the numerical parameters values for Han's model are taken from~\cite{Grenier2020} and given in Table~\ref{tab:parameter}.
\begin{table}[htpb]
\caption{Parameter values for Han Model.}
\begin{center}
\begin{tabular}{|c|c|c|}
\hline
$k_r$  & 6.8 $10^{-3}$ & s$^{-1}$\\
\hline
$k_d$ & 2.99 $10^{-4}$  & -\\
\hline
$\tau$ & 0.25 & s\\
\hline
$\sigma$ & 0.047 & m$^2\cdot(\mu$mol)$^{-1}$\\
\hline
$k$  & 8.7 $10^{-6}$ & -\\
\hline
$R$ &  1.389 $10^{-7}$ & s$^{-1}$\\
\hline
\end{tabular}
\end{center}
\label{tab:parameter}
\end{table}

In order to determinate the light extinction coefficient $\varepsilon$, let us assume that only $1\%$ of light can be captured by the cells at the bottom of the raceway, i.e. for $z=z_b$, meaning that $I_b=0.01I_s$, we choose $I_s=\SI{2000}{\mu mol.m^{-2}.s^{-1}}$ which approximates the maximum light intensity, e.g., at summer in the south of France. Then $\varepsilon$ can be computed by
\begin{equation*}
\varepsilon = (1/h(0;a))\ln(I_s/I_b).
\end{equation*}

Note that in the case where $a_0$ is also a parameter to be optimized, we take from~\cite{Martinez201811} the specific light extinction coefficient of the microalgae specie $\alpha_0=\SI{0.2}{m^2\cdot g}$C the background turbidity and $\alpha_1=\SI{10}{m^{-1}}$, under these settings, $a^*_0$ defined in~\eqref{eq:a0} is around $\SI{0.316}{m}$.





\subsection{Numerical results}

We now test the influence of various parameters on optimal topographies. In all our experiments, we always observe that the topographies satisfy $\norm{h}_\infty> h_c$.

\subsubsection{Influence of vertical discretization}

The first test consists in studying the influence of the vertical discretization number $N_z$. We choose $N=5$ and take 100 random initial guesses of $a$. Note that the choice of $a$ should respect the subcritical condition. Let $N_z$ varies from $1$ to $100$, and we compute the average value of $\bar \mu_{\Delta}$ for each $N_z$. The results are shown in Fig.~\ref{fig:cvNz}.
\begin{figure}[thpb]
\centering
\includegraphics[scale=0.4]{Is_2000_q_001_N_5_L_100_dx_001_convergence_Nz_optim_h_x}
\caption{The value of the functional $\bar \mu_{\Delta}$ for $N_z=[1,100]$.}
\label{fig:cvNz}
\end{figure}
We observe numerical convergence when $N_z$ grows, showing the convergence towards the model continuous in space. In view of these results, we take hereafter $N_z=50$.





\subsubsection{Influence of the initial condition}\label{sec:testC0}

The second test consists in studying the influence of the initial condition $C_0$ on the optimal shape of the raceway pond. We set the numerical tolerance Tol$=10^{-10}$, and consider $N=5$ terms in the truncated Fourier series as an example to show the shape of the optimal topography. As an initial guess, we consider the flat topography, meaning that $a$ is set to 0. Besides, we compare the optimal topographies obtained with $C_0=0.1$ and with $C_0=0.9$. The result is shown in Fig.~\ref{fig:C0depend}.
\begin{figure}[thpb]
\centering
\includegraphics[scale=0.4]{Is_2000_q_001_N_5_Nz_50_L_100_dx_001_C0_01_optim_topo_optim_h_x.pdf}
\includegraphics[scale=0.4]{Is_2000_q_001_N_5_Nz_50_L_100_dx_001_C0_09_optim_topo_optim_h_x.pdf}
\caption{The optimal topography for $C_0=0.1$ (top) and $C_0=0.9$ (bottom). The red thick line represents the topography ($z_b$), the blue thick line represents the free surface ($\eta$), and all the other curves between represent the different trajectories.}
\label{fig:C0depend}
\end{figure}
A slight difference between the two optimal topographies is observed. Note that this difference remains when the spatial increment $\Delta x$ goes to zero.





\subsubsection{Influence of Fourier series truncation}

The next test is dedicated to the influence of the order of the truncation of the Fourier series used to parameterize the water elevation $h$. 

Set $N=[0,5,10,15,20]$ and keep all the other parameter values. Table~\ref{tab:N} shows the optimal value of the objective function $\bar \mu_{\Delta}(a^*)$ for different order of truncation $N$ of the Fourier series.
\begin{table}[htpb]
\caption{The objective function for different order of truncation of the Fourier series.}
\begin{center}
\begin{tabular}{|c|c|c|c|}
\hline
$N$ & Iter & $\bar \mu_{\Delta}(\si{d^{-1}})$ & $\log_{10}(\|\nabla J\|)$\\
\hline
0 & 0 & 0.94764 & $-$ \\
\hline
5 & 16 & 0.96942 & -10.368944 \\
\hline
10 & 17 & 0.97509 & -10.262891 \\
\hline
15 & 17 & 0.97734 & -10.010322 \\
\hline
20 & 18 & 0.97851 & -10.196963\\
\hline
\end{tabular}
\end{center}
\label{tab:N}
\end{table}
As expected, the result shows a slight increase of the optimal value of the objective function $\bar \mu_{\Delta}(a^*)$ when $N$ becomes larger. However, the corresponding values of $\bar \mu_{\Delta}(a^*)$ remain close to the one associated with a flat topography.





\subsubsection{Simulation with paddle wheel}

The last test aims to simulate a more realistic raceway pond situation, where a paddle wheel considered in the system. More precisely, we simulate several laps, with a paddle wheel that mix up the algae after each laps. The paddle wheel is modelled by a permutation matrix $P$ that rearrange the trajectories at each lap. In our test, $P$ is an anti-diagonal matrix with the entries one. This choice actually corresponds to an optimum and has been shown in~\cite{Bernard02970756}, where other choices will be investigated. 

Meanwhile, this permutation matrix $P$ corresponds to the permutation
\begin{equation*}
\pi = (1 \ N_z)(2 \ N_z-1)(3 \ N_z-2)\cdots,
\end{equation*}
where we use the standard notation of cycles in the symmetric group. Note that $\pi$ is of order two. The photo-inhibition state $C$ is then set to be $2$-periodic (i.e $C^1(0)=PC^2(L)$). The details of the optimization procedure are given in Appendix B.

We choose a truncation of order $N=5$ in the Fourier series. The initial guess $a$ is still zero. Fig.~\ref{fig:2laps} presents the shape of the optimal topography and the evolution of the photo-inhibition state $C$ for two laps. In particular, we observe that the state $C$ is actually periodic for each laps. This result is actually proved in~\cite{Bernard02970756} in the case of a flat topography.
\begin{figure}[thpb]
\centering
\includegraphics[scale=0.4]{Is_2000_q_001_N_5_Nz_50_Nl_2_L_100_dx_001_C0_per_optim_topo_optim_h_matrix.pdf}
\includegraphics[scale=0.4]{Is_2000_q_001_N_5_Nz_50_Nl_2_L_100_dx_001_C0_per_C_optim_h_matrix.pdf}
\caption{The optimal topography (top) and the evolution of the photo-inhibition state $C$ (bottom) for two laps.}
\label{fig:2laps}
\end{figure}

The resulting optimal topography in this case is not flat. However, the increase in the optimal value of the objective function $\bar \mu_{\Delta}$ compared to a flat topography is around $0.217\%$, and compare to a flat topography and non permutation case is around $0.265\%$, in both cases, the increase remains small. 





\section{CONCLUSIONS AND FUTURE WORKS}

A non flat topography slightly enhances the average growth rate. However the gain remains very limited and it is not clear if the difficulty to design such a pattern could be compensated by the increase in the process productivity.

A flat topography cancels the gradient of the objective functional in many situations where $C$ is assumed to be periodic. However, including a mixing device gives rise to an optimal non flat topography with a slight gain of the average growth rate. 


\bibliographystyle{plain}       
