The light/dark cycle of microalgae in a thin-layer photobioreactor

A numerical study of the motion of algal cells in a representative thin-layer-cascade (TLC) photobioreactor is presented. The goal is to determine the time scale associated with the light/dark (L/D) cycle seen by the cells during their turbulent motion in the liquid culture. Owing to the limited reliability of the available numerical results which deal with time-averaged quantities and thus lack time-resolved information, the present study is based upon the Direct Numerical Simulation of the Navier-Stokes equations, a reliable but consequently expensive numerical approach which does not incur in turbulence modelling errors. Indeed, the simulation is successfully validated in terms of averaged velocity with experimental data. The availability of full temporal information allows algae cells to be followed in time along their trajectories. A large number (up to a million) of tracers is placed in the flow to mimic the algae cell. Their trajectories are statistically studied and linked to the turbulent mixing. Results indicate that, in a typical TLC reactor designed to mimic an experimental setup, cells undergo an L/D cycle with a time scale in the range 0.1–2 s. Such time scale, albeit much longer than the typical time scale of the photosynthesis, significantly benefits the productivity of the algae compared to a steady illumination.


Introduction
A strong interest in cultivation of phototrophic microalgae is driven by their several envisaged application areas: a non-exhaustive list includes biofuels, chemicals, medicine and nutrition (Posten and Walter 2012). Microalgae hold potential in food production and in purification of contaminated water and wastewater (Acién et al. 2016;Bernaerts et al. 2017). They are important in the food chain as they are the main feeding source for rotifers, fishes, etc.; microalgae possess high nutritional value and contain lipids, proteins, carotene and other essential minerals (Hakalin et al. 2014). Microalgae are also used as fertilisers and in products for skin health. Over the last years, owing to the decline of fossil fuel supply and to the increasing concerns on global warming from Alessandro Chiarini alessandro.chiarini@polimi.it Maurizio Quadrio maurizio.quadrio@polimi.it 1 Department of Aerospace Sciences and Technologies, Politecnico di Milano, Milan, Italy excessive carbon dioxide emissions, microalgae became interesting also as an alternative source of renewable energy (Makarevičiene et al. 2011;Gour et al. 2016).
To date, the most widely used photobioreactor (PBR) for large-scale cultivation of microalgae is of the openpond type (Grima et al. 1999). An open-pond PBR is simple and enjoys low construction costs and low power requirements. Furthermore, it can be easily scaled up in size; cleaning and maintenance are simple and affordable. Unfortunately, however, the open-pond PBR has a low production efficiency. This is mainly due (Davies et al. 2016) to unpredictable and/or uncontrollable algal growth conditions with respect to nutrient quality and quantity, pH, salinity, temperature and turbulence in the flow. Above all, the conventional open-pond PBR suffers from poor light-to-biomass conversion efficiency. How light is delivered to microalgae cultivation is of paramount importance (Richmond 2004;Posten 2009). Conventional PBR typically employ natural illumination. At daytime, the sunlight intensity is usually excessive near the irradiated liquid surface and causes photoinibition (see Han et al. (2000) and references therein), whereas a large fluid volume beneath the surface receives insufficient light for photosynthesis, leading of the so-called respiration phenomenon, where net photosynthesis cannot take place (Iluz and Dubinsky 2012).
Strategies do exist to increase the light-to-biomass conversion efficiency. A popular approach is for example to deliver light to the algae cells according to a certain temporal light/dark (L/D) cycle that best suits their internal characteristics. The photosynthetic efficiency of microalgae under intermittent illumination is known to be higher than that under continuous illumination, provided the parameters of the L/D cycle are tuned correctly; see for example (Kok 1953;Janssen et al. 2000a, b;Grobbelaar 2006;Xue et al. 2011). This is because photosynthesis itself is a cyclic process, where almost instantaneous photic reactions are followed by slower thermochemical ones. The time scale of the photosynthetic cycle is thus determined by the slower thermochemical reaction and is set at approximately 1-15 ms (Richmond 2004), implying that the optimal L/D frequency is around 70-1000 Hz. Nedbal et al. (1996) experimentally investigated the effect of intermittent illumination in dilute cultures of different species of algae using arrays of red LEDs. They found that for periods of the L/D cycle shorter than 1 ms, the growth rate is close to that found in the continuous light regime. For periods in the range 1-100 ms, instead, the growth rate is 10-20% higher than in the continuous case. Moreover, further experiments carried out in short-light-path PBRs have found that L/D cycles with frequencies much lower than the typical frequency of the thermochemical reactions, i.e. in the range of 0.5-1 Hz, are still very beneficial to the algal growth rate and light conversion efficiency and therefore to the productivity (Miller et al. 1964;Laws et al. 1983;Vejrazka et al. 2011;Liao et al. 2014).
Several groups are thus trying to exploit such a flashinglight principle to enhance productivity. Some strategies involve the optical modification of the irradiated surface, as done for example by Liao et al. (2014), who applied a spatially alternate shading to a tubular PBR. Most of the proposed solutions aim at creating L/D cycles through the design of a PBR with suitable mixing properties, such that algal cells are periodically moved to and from the illuminated regions of the reactor. Clearly, the spatial extension of the illuminated and dark regions depends on the concentration of the algae suspension (Richmond 2003): the thickness of the illuminated region becomes smaller when the culture density increases, as the algae cells absorb more of the incident light.
Among the mixing-type PBR designed for enhanced efficiency, the Taylor-Couette and the thin-layer cascade reactors are worth citing. In the former, the algae cultivation is confined to the annular region between two concentric cylinders, with the inner cylinder rotating (Miller et al. 1964). By a fluid instability mechanism (Taylor 1923;Diprima 1960), this configuration efficiently produces large-scale toroidal vortices in the fluid, known as Taylor vortices, which shuttle the fluid between the dark region close to the inner cylinder (if the PBR receives illumination from the outer ambient) and the illuminated region close to the outer cylinder. The process can be energetically efficient because of the instability mechanism that creates and sustains the Taylor vortices. It has been demonstrated by several works (Kong and Vigil 2013a, b;Gao et al. 2015aGao et al. , b, 2017a that in a Taylor-Couette PBR, the biomass yield can be significantly higher than that obtained without Taylor vortices (i.e. the same setup without cylinder rotation).
In principle, the thin-layer cascade (TLC) PBR consists of a microalgae suspension flowing down a sloped wall, forming a thin liquid layer with thickness typically less than 1 cm. At the bottom of the wall, the suspension is pumped up back to the starting point. Experimental information (Doucha and Lavanský 1995, 2006, 2009Masojídek et al. 2011;Apel et al. 2017) provides evidence of high biomass productivity by using a TLC PBR. The cost of a TLC is estimated to be significantly lower (approximately 5 times) than the cost of a conventional open-pond raceway (Doucha and Lívanský 2009). Other significant advantages of the TLC PBR derive from the possibility of employing a much denser suspension, which leads to significant savings in the final phase of cell harvesting. In the TLC PBR, the mixing process is believed to be driven by turbulence fluctuations, which shuttle the algae cells from the dark near-wall region to the illuminated region close to the irradiated surface, owing to the turbulent nature of the flow. However, no evidence is available to link the turbulent motions in the suspension with the properties of the mixing-induced L/D cycles in the TLC and its productivity features. As a result, the choice of the best mixing mode for a TLC, and the detailed design of an efficient PBR based on its mixing properties, still are open problems.
To address such issues one needs to resort to a detailed understanding of the mixing properties of the system. This can be achieved via either experimental or numerical studies. Within the first class, we mention in particular the excellent study by Apel et al. (2017), which is numerically replicated in the present work: in that study, the productivity characteristics of a TLC are studied under very controlled and realistic environmental conditions. Numerical studies, albeit few, are becoming progressively more popular, as they allow a better insight in the physics of the suspension. Among them we recall the work of Severin et al. (2018), where the vertical mixing in the main cultivation channel of the same TLC is investigated for different operating conditions and fluids using Computational Fluid Dynamics (CFD) and turbulence modelling. They also introduced a measure of the mixing capability of a TLC, although because of limitations of their model, it does not account for the typical meandering features of each and every trajectory of the algal cells. More recently, this work has been extended by Akhtar et al. (2020) to the full TLC. They employ CFD to investigate separately the hydrodynamic performance and the mixing of algae cells in each compartment of the same reactor and address the effect of the channel width, water depth and reactor slope. However, they are still limited by the use of turbulence modelling and by the lack of time-resolved information. Chen et al. (2016) investigated the L/D cycle in raceway ponds via CFD and addressed the influence of the paddlewheel rotational speed and the effect of flow-deflection baffles. They found that in a typical open pond, the mean L/D cycle is about 20 s and that it decreases as the paddlewheel rotational speed increases. Perner-Nochta and Posten (2007) studied the L/D cycle in a tubular PBR with and without static mixers. They found that in such tubular PBR, a static mixer enhances the mixing and optimises the L/D cycle. By using CFD modelling Huang et al. (2014) developed a novel mixer to improve the biomass productivity of an air-lift flat-plate PBR. Sato et al. (2010) developed a virtual PBR which combines a numerical simulation of two-phase turbulent flow and a photosynthesis model to characterise the light distribution in the PBR. Regarding the Taylor-Couette PBR, we have cited above the studies of Gao and coworkers: they performed both experiments and numerical simulations.
For other examples, see Pires et al. (2017) and references therein.
Despite the several attempts, all these works are limited by (and dependent upon) their modelling choices, in particular turbulence modelling. Most of the existing turbulence models, in fact, are derived from high-Reynoldsnumber flows and designed for high-Reynolds-number industrial or aeronautical applications; their use to model the low-Reynolds-number flow in a PBR requires extra care. Moreover, for the TLC PBR, the typical computational approaches are based on a temporal average and do not resolve the flow in time, preventing the investigation of the mixing. Therefore, to the best of our knowledge, to date none has been capable to discuss the typical L/D cycle of the algae cells in a specific TLC PBR: the lack of fully reliable numerical simulations prevents a complete understanding of the physics of such flows and of the mechanisms that may enhance the microalgae productivity.
In this work, we study through numerical simulations one specific TLC photobioreactor, and aim at understanding its working principle, by characterising the L/D cycle the algal cells are typically subjected to. To this end we resort to the most reliable numerical approach, the Direct Numerical Simulation of the equations of motion, which-at the cost of a large computational cost-can handle turbulence properly without the need for models of any sort and completely resolves the flow in time. We will reproduce the operating conditions of the cultivation channels of the TLC PBR described in the reference work of Apel et al. (2017) and complement their results and those of Severin et al. (2018) and Akhtar et al. (2020). In particular, by inserting in the flow a large number of algae cells and by tracking them in time, we will be able to statistically characterise their trajectories and to describe the properties of their L/D cycle.

Methods
The TLC experimentally described by Apel et al. (2017) is considered in the present work. Our focus is on the core part of one channel of the PBR; hence, turns, steps and the return circuit are neglected. These are important factors determining the global energetic efficiency of the PBR, but do not affect the light conversion mechanism. The most important part of the PBR, in the context of the present work, is the thin layer of fluid that flows down the inclined plate under the action of gravity. The fluid dynamical characterisation of the flow is carried out by resorting to the most sophisticated numerical approach available, i.e. the Direct Numerical Simulation (DNS) of the Navier-Stokes equations. In contrast to alternative approaches for the numerical simulation of fluid flow, DNS is an expensive three-dimensional and timedependent approach which dispenses with the need of turbulence modelling, at the cost of resolving the smallest dynamically significant spatial and temporal scales in the flow. As such, it is particularly well suited to fundamental studies aimed at understanding the physical processes involved.
The experimental setup by Apel et al. (2017) is replicated. The numerical values of the key physical and geometrical parameters are provided in Table 1. In particular, as in the original reference, density ρ and kinematic viscosity ν of the algae suspension are those of seawater at 20 • C. The goal is to investigate the properties of the L/D cycle of the algae cells, which are considered in the simulations as particle tracers whose meandering motion is characterised individually and in a statistical sense.

The numerical simulation
The velocity field in the TLC is computed via a Direct Numerical Simulation (DNS) of a possibly turbulent flow in an indefinite open channel (see Fig. 1), where a fluid layer of uniform thickness h flows down an inclined flat plate placed at an angle α with the horizontal under the driving action of gravity. Since the extent of the inclined plate is much larger than the dimensions of the largest turbulent structure, the flow can be safely considered statistically homogeneous in the wall-parallel directions, and only a portion of the whole layer needs to be simulated. It is well known (Pope 2000) that in this configuration, the flow statistics only depend upon the distance from the wall and upon the Reynolds number. The latter is a dimensionless group that can be interpreted as the ratio of the inertia and viscous forces, or (as in the present case) as the ratio of the molecular and turbulent diffusive time scales for a fixed lengthscale. Its generic definition is: where U and L are the reference velocity and length scales of the problem and ν is the kinematic viscosity of the fluid. Simply put, at low values of Re, the flow is laminar, while at higher Re, the flow becomes turbulent. Two Reynolds numbers are appropriate here, depending on the choice of the reference velocity, with the layer thickness h setting the Fig. 1 Sketch of the open-channel geometry and the reference system. The blue/coloured surface represents the liquid surface, parallel to the wall (black line) and inclined by an angle α with respect to the horizontal (dashed surface) lengthscale L. A possible velocity scale is the bulk velocity where u is the mean velocity profile, a function of the sole wall-normal coordinate y, and the operator · denotes averaging in the homogeneous streamwise and spanwise directions and over time. It leads to the definition of a bulk Reynolds number: An alternate velocity scale is the friction velocity u τ : where τ w is the mean wall shear stress, and ρ is the fluid density. It leads to the definition of a friction Reynolds number: Re b can be seen as a dimensionless measure of the flow rate, whereas Re τ expresses in a dimensionless way the frictional force exerted by the flow upon the wall.
The downward inclination α of the thin-layer flat plate translates the force of gravity into an uniform acceleration volume force, which is equivalent to a uniform streamwise pressure gradient along the channel. At equilibrium, it is balanced by the friction force at the wall. Hence: implying that the friction velocity is constant along the plate. From the numerical values for h and α reported in Table 1, one obtains: and this sets Re τ = 163 as the Reynolds number of the present simulations. This value of Re τ is thus provided as a simulation input, whereas the flow field, its statistics, including the mean velocity profile, are simulation outputs. The evolution of the flow is governed by the incompressible Navier-Stokes equations, written here in their non-dimensional form: where repeated indices imply summation; x = x 1 (u = u 1 ), y = x 2 (v = u 2 ) and z = x 3 (w = u 3 ) are the streamwise, wall-normal and spanwise directions (velocity components), and p is the pressure. The Kronecker symbol δ i1 implies that the external pressure gradient only affects the streamwise component. All quantities are made dimensionless with u τ , h and ν; hence, the reference Reynolds number is Re τ . The in-house solver employed in the present study is a state-of-the-art pseudo-spectral DNS code, originally introduced and validated by Luchini and Quadrio (2006), and since then used in a large number of scientific works. The code solves the incompressible Navier-Stokes equations projected in the divergence-free space of the wallnormal components of the velocity and vorticity vectors, as suggested by Kim et al. (1987). The solver has parallelcomputing capabilities, uses a mixed spatial discretisation, with Fourier series for the homogeneous directions, where periodic boundary conditions are used, and fourth-order, explicit compact difference schemes in the wall-normal direction. At the wall, no-slip u = w = 0 and nopenetration v = 0 conditions are imposed, whereas at the fixed, planar interface, a free-slip condition is used, with v = 0 and ∂u/∂y = ∂w/∂y = 0.
The size of the computational domain is L x = 4πh and L z = 2πh (corresponding to a portion of 7.02 cm and 3.51 cm) in the streamwise and spanwise directions, where the discretisation employs N x = N z = 256 Fourier modes (the 3/2 rule is used to remove aliasing error in the nonlinear operations). The domain size is less than the actual length and width of the PBR; since the pioneering DNS work by Kim et al. (1987), we know that for statistically homogeneous flows, these dimensions are more than adequate to provide the periodic boundary conditions in the homogeneous directions with physical significance, and to properly compute statistical quantities to characterise the turbulent flow. In the wall-normal direction, the differential operators are discretised via fourth-order compact finite differences using 50 points on a non-uniform grid refined near the wall via an hyperbolic tangent distribution: where y 0 = 0 denotes the channel wall and a = 1.6. The spatial resolution fully satisfies the requirements for a wellresolved DNS that have been determined in the seminal work of Kim et al. (1987) cited above. They found that the grids spacing in the streamwise and spanwise direction need to be no larger than δx + ≤ 10 in the streamwise direction and δz + ≤ 5 (the + superscript denotes quantities expressed in wall units, i.e. made non-dimensional with u τ and ν); in the present work, we have δx + ≈ 8 and δz + ≈ 4. As for the wall-normal discretisation, the optimum is to have the cell closest to the wall not larger than δy + ≤ 1, but not much smaller owing to the penalty in computing costs; in the present work, the wall-normal size of the first near-wall cell is δy The solution is advanced in time by a partially implicit method, where the viscous terms are computed with a second-order-accurate Crank-Nicolson method, whereas a third-order Runge-Kutta scheme deals with the convective terms. The simulation is run for a time interval of 30h/u τ (corresponding to a physical time of 5.5 s) after reaching statistical equilibrium. The time step dynamically varies during the simulation in order to satisfy the CFL condition, i.e. the maximum Courant-Friederichs-Lewy CFL number is kept at CFL = 1 in the whole domain during the simulations. The resulting time step fluctuates slightly over time and averages at 0.000209h/u τ , corresponding to 0.038 ms.

Particles and particle tracking
The motion of the algae cells suspended in water is simulated by initially placing a large number of particles in the flow domain and then following their evolution in time. The strategy to track the algae cells within the suspension in time depends on the particle details. In particular, whether particles simply follow the fluid stream without affecting the flow field depends on a dimensionless quantity called the Stokes number Stk, defined as: where τ p , ρ p and d p are the relaxation time, the density and the diameter of the algae cell. If Stk < 1 particles follow the fluid stream without affecting the flow field (see McKeon et al. (2007) pp. [288][289]. In this study, we consider singlecell spherical-shaped cells of Chlorella with diameter d p = 7μm and density ρ p = 864 kg m −3 as reported in Williams and Laurens (2010), Ali et al. (2015), and Akhtar et al. (2020). For this algae cell species Stk 1, implying that in this flow conditions the Chlorella cells may be safely modelled as lagrangian tracers.
Eventually, the trajectories of the whole particles set are brought together for statistical characterisation. To assess the statistical convergence of the trajectory statistics, the total number N of particles is varied from 10 4 to 10 6 . The effect of the number of particles on the computed statistics is discussed later in the 'Discussion' section. The trajectory of each particle is computed via the numerical integration of the differential equation with the initial condition x(x 0 , 0) = x 0 , x 0 is the initial position of the particle and v(x 0 , t) is the particle velocity at time t. The lagrangian particle velocity v is related to the eulerian velocity u of the fluid by: v(x 0 , t) = u(x(x 0 , t), t). and the temporal integration of Eq. 7 above uses the same explicit Runge-Kutta integration scheme of the DNS code, as well as the same discretisation. During the DNS, the Eulerian velocity field u(x, t) is available at each timestep at discrete locations in the threedimensional computational grid. Since the instantaneous position of a particle does not normally coincide with a grid point, the velocity of every particle must be computed via a three-dimensional interpolation of the DNS velocity field, and it is known that such interpolation procedure may critically affect the accuracy of the procedure (Yeung and Pope 1988;Bernardini 2014). To this aim, the interpolation is carried out by using Lagrange polynomials of sixth order. The interpolation is locally reduced to fourth order for particles positioned between the wall and the first internal grid point.

The light/dark cycle and its statistical characterisation
Each particle can be followed along its trajectory to ascertain the amount of time spent in the illuminated volume of the fluid or in the dark one. The position of the boundary separating the dark and illuminated regions will be discussed below. After following each particle along its trajectory, every ith particle can be associated with its temporally averaged frequency f i of L/D switching, as well as the duty-cycle φ i , i.e. the fraction of time spent by the particle in the illuminated part of the L/D cycle. Accumulating such information for the whole particle set leads to a statistical characterisation of the light/dark cycle.
The average frequency f av and the average duty cycle φ av of the L/D cycle are the main quantities of interest for a statistical characterisation. They are computed as the arithmetic average of f i and φ i for all the N particles inserted in the flow, i.e.:

Velocity field
The turbulent velocity field is decomposed according to the classic Reynolds decomposition as: where U ≡ u is the mean field and u is the fluctuating field. The present flow is statistically stationary and has only one non-homogeneous direction, i.e. the wall-normal direction y. Therefore, the mean velocity vector field Under the present working conditions, namely at Re τ = 163, the flow is certainly turbulent; hence, a threedimensional and unsteady fluctuating velocity field u = 0 exists. This is easily visualised in Fig. 3 where the wallnormal profiles of the auto-and cross-covariances of the velocity fluctuations are plotted, which would be all zero in Fig. 3 Wall-normal distribution of the four non-zero components of the Reynolds stresses tensor u i u j a laminar flow. Four profiles are drawn, corresponding to the non-zero components of the symmetric tensor u i u j , also known as the tensor of Reynolds stresses. The components u w and v w are zero owing to the statistical symmetry of the flow with respect to an inversion of the z-axis.
Intense turbulent activity, corresponding to large u i u j , is observed near the wall at y + ≈ 15. Indeed, production of turbulent fluctuations is associated with the flow anisotropy (Pope 2000), which it is stronger in the vicinity of the wall. Moreover, one can notice that u u is larger than v v and w w . It is known that the mean flow, aligned to the x direction, triggers the largest fluctuations in the same direction. Since u w = v w = 0, the only off-diagonal component of the tensor is u v . It is proportional to the turbulent shear stress and is responsible for the deviation of the mean velocity profile from the laminar parabolic Poiseuille profile. A further insight into the complex turbulent motions can be obtained from   (1995). It is based on the second largest eigenvalue λ 2 of the tensor: are the symmetric and antisymmetric part of the velocity gradient tensor ∂u i /∂x j . This method accurately extracts the three-dimensional pattern of the vortical structures in the flow (Jeong and Hussain 1995;Cucitore et al. 1999) whose internal core is identified as the isosurface where the intermediate eigenvalue λ 2 becomes negative. The turbulent structures are observed to preferentially populate the nearwall region, confirming that the main turbulent activity in wall-bounded flows occurs in the vicinity of the solid wall. The presented pattern agrees with the picture of elongated vortices aligned with the streamwise direction typical of near-wall turbulence, which are intermixed with elongated stripes or streaks of low longitudinal velocity, whose extent is comparable to the length of the computational domain. The interaction between 'streaks' and 'quasistreamwise vortices' is essential to the near-wall turbulence regeneration cycle (Jiménez and Pinelli 1999) and cannot be accounted for by the simpler CFD simulations, such as RANS or URANS, based on the time average.

Lagrangian particles
Unlike in a laminar channel, where the velocity field is parallel to the streamwise direction, in a turbulent channel, the velocity field is unsteady and highly three-dimensional. As a consequence, particles released in a turbulent flow follow intricate three-dimensional trajectories, each moving with an apparently random pattern along the z and y directions while being convected along the streamwise direction.
In the present work, we are interested in the temporal evolution of the wall-normal coordinate of each particle, which identifies the particle position relative to the photic layer adjacent to the upper side of the computational domain. Given the density of the culture and the illumination conditions, one may assume that only the instantaneous wall-normal particle position determines the external radiation received by the algae cell. The nature of such evolution is qualitatively shown in Fig. 5, where the wall-normal position of a few randomly chosen particles is shown after projection onto a z − y plane at a certain time t. Only a few particles are shown for clarity and colour-coded with their initial distance from the wall. The figure demonstrates that particles loose memory of their initial wall-normal position. Indeed, it can be theoretically shown that, in a turbulent flow, ideal massless lagrangian tracers become eventually fully mixed by the velocity field in the limit t → ∞ (Richardson 1926;Pitton et al. 2012). This is clearly visible in video 2 and video 3 in Supplementary Material, which show the motion of the lagrangian tracers in time in the three-dimensional space and in a cross-stream plane.
To statistically investigate the properties of the L/D cycle faced by the algae cells owing to the turbulent fluctuations, which is the subject of the next section, one needs to evaluate for each particle the time history of its wall-normal position y p . Owing to the presence of non-zero vertical velocity fluctuations v , y p varies continuously in time. Of particular importance is the period of time that every algal cell spends near the free surface or in the bulk and in the near-wall regions, i.e. the light and dark regions of the culture, as well as the frequency by which the particle switches between those regions. Figure 6 shows an example of the time history of the wall-normal position of one particle. If the thickness of the euphotic layer is known, from such time history, one can immediately compute the fraction of time spent by the particle inside the euphotic layer and the corresponding duty cycle φ i , as well as the number of times the particle enters or exits the euphotic layer and the corresponding frequency f i (see video 3 in Supplementary Material where the motion of 8 particles in a cross-stream plane is shown and the euphotic layer is highlighted). From such information, and after properly averaging over the entire particle set, one arrives at the mean Fig. 6 Time history of the wall-normal distance y p of a lagrangian particle. The red dotted line identifies the euphotic layer corresponding to the low-density example described in the 'Low-density culture and/or high-intensity illumination' section values f av and φ av for the frequency and the duty cycle of the L/D cycle, as defined by Eq. 8.

Discussion
The trajectories of the particles within the turbulent flow contain information about the motion of the algae cells in the TLC PBR. Because of the vertical motion induced by the turbulent fluctuations, the algae cells continuously move from light to dark regions, and vice-versa. The average frequency f av and the average duty cycle φ av of their L/D cycle primarily depend upon the following factors: (i) the turbulent intensity of the flow, (ii) the intensity of the light source and (iii) the culture density. Factor (i) determines the intensity of the driving wall-normal turbulent fluctuations, whereas the combination of factors (ii) and (iii) sets the thickness of the euphotic layer in the suspension.
The dependence of the photosynthetically active radiation intensity I par on the depth δ of the optical path can be estimated through the Beer-Lambert law (Parnis and Oldham 2013): where I 0 , expressed in W m −2 , is the intensity of the incident light on the culture, η par is the percentage of the total radiation which is photosynthetically active, is the absorption coefficient expressed in L g −1 cm −1 and c is the concentration of the algae cells expressed in g L −1 . Following Luhtala and Tolvanen (2013), the thickness of the euphotic layer is conventionally defined as the distance from the surface where the light intensity drops to 1% of the surface value.
In the following, two different examples are considered. In the first, a relatively thick euphotic layer corresponds to a culture with low density and/or with high illumination intensity. At the other extreme, a thin euphotic layer corresponds to an high-density culture and/or low illumination intensity. In both cases, once the thickness δ, or equivalently the position y δ = h−δ of the euphotic layer is set, the frequency and duty-cycle of the L/D cycle seen by the particles are computed by following in time the instantaneous wall-normal position y p of each particle: when y p > y δ the particle is said to be within the euphotic layer, whereas when y p < y δ it is said to be in the dark region. Hence, y = y δ is the distance from the wall which defines the plane separating the illuminated and dark flow regions. The L/D frequency is computed as one half the number of times the particles cross the plane y = y δ during the simulation time, divided by the entire simulation time; the duty-cycle results from the ratio between the time spent in each of the two flow regions.

Low-density culture and/or high-intensity illumination
As an example of culture with thick euphotic layer, we consider a culture of Chlorella with density c of 5 g L −1 . The light source is sunlight with intensity of I 0 = 800 W m −2 . In this and the following example, the absorption coefficient has been set at = 0.0034 L mg −1 cm −1 and η par = 0.4 as in Apel et al. (2015Apel et al. ( , 2017. This results in a euphotic layer thickness of 2.7 mm, i.e. approximately one-half of the overall suspension thickness, as shown in Fig. 7, where the change of irradiance along the optical path is plotted. In this case, a combined analysis of the trajectories of the whole set of particles leads to a computed mean frequency of the L/D cycle of f av = 1.62 Hz. Moreover, the simulation indicates that, on average, the algae cells spend almost one half of their time exposed to light, i.e. φ av = 48.6%. This means that in such cultures, the algae cells move from dark to light regions and viceversa with a time scale of approximately 0.6 s; during their motion, the algae cells in average spend almost the same amount of time in illuminated and in dark regions. Average values only provide a partial description of the L/D cycle, but the rich DNS dataset, which contains full temporal information, can be perused to extract more complete information in terms of the Probability Density Function (PDF) of the quantities of interest. A more informative look at the data is thus obtained via Fig. 8. Its top panel presents an histogram representation of the PDF of the frequency of the L/D cycle, and the bottom panel plots a similar histogram for the duty cycle. The figure reveals for example that 90.7% of the cells have an L/D frequency larger than 1 Hz, a percentage that drops to 23.1% if the threshold of 2 Hz is considered. Analogously, 53.3% of the cells have a duty cycle between 40 and 60%, and 3.4% only have a duty cycle between 0 and 20%.

High-density culture
The opposite case of a culture with high density and/or less intense light source involves an euphotic layer of reduced thickness. In the example that follows, the same conditions as in the previous case are considered but with a density culture of 40 g L −1 . This results is an euphotic layer thickness of only 0.34 mm, i.e. 6% of the overall suspension thickness, as shown in Fig. 9, where the irradiance is plotted.
In this case, a combined analysis of the trajectories of the whole set of particles leads to computing for the L/D cycle a mean frequency of f av = 0.51 Hz. Moreover, data indicate that, on average, the algae cells spend less than 1/10 of their time exposed to light, namely φ av = 6.1%. This means that in such cultures, the algae cells move from dark to light regions and vice versa with a time scale of approximately 1.9 s; during their motion, the time spent by the algae cells in the illuminated region is more than one order of magnitude larger than the time spent in dark regions.
Again, a more informative way of looking at temporal data consists in looking at the PDF, shown by Fig. 10. Its top panel presents an histogram of the PDF of the L/D cycle frequency, and the bottom panel plots the histogram of the duty cycle. The figure reveals for example that 11.4% of the cells have an L/D frequency larger than 1 Hz, a percentage that drops to less than 0.1% if a threshold of 2 Hz is considered. Analogously, 69.3% of the cells have a duty cycle between 0 and 10%, and 20.9% have a duty cycle between 10 and 30%.
In conclusion, Fig. 11 shows that these results are robust with respect to the number of particles placed in the flow and used for computing the statistics. If the number of particles is reduced from 10 6 to 10 4 , a variation of less than 0.8% is observed for f av and a variation less than 1.7% for φ av . Moreover, the PDF of the frequency obtained using 10 5 particles shows negligible differences with respect to that obtained with the whole set of 10 6 particles and previously shown in Fig. 10.   Fig. 9 As for Fig. 7, but for a density of 40 g L −1 Fig. 10 As for Fig. 8, but for high-density culture (or low light intensity)

Comments
The main features and limitations implied by our model are recalled and further discussed here. The model is based on a numerical approach that is free from any turbulence modelling errors, and employs discretisation and computational procedures which are state of the art. However, it is still based on hypotheses and affected by limitations.
1. The flow dynamics within the TLC photobioreactor is studied numerically in an idealised setup of an openchannel gravity-driven flow with constant depth. The air-water interface is flat and does not deform, and the extent of the computational domain is smaller than the real flat plate, which is only a part of an actual TLC. The use of periodic boundary conditions for the homogeneous direction in the DNS of wall-bounded turbulent flows is well accepted in the literature and provides reliable turbulence statistics for a finite computational domain, provided that the domain is large enough to represent the largest flow structures.
Although not shown, this is verified in the present case, according to the standard check of verifying that all the correlation functions drop to zero within the separation of one-half of the computational domain length. Centre: duty cycle of the L/D cycle. Bottom: PDF of the frequency for 10 5 particles (to be compared with top panel in Fig. 10) 2. The algae cells are modelled as massless tracers which are passively driven around by the velocity field. This assumption is fully justified on consideration of the negligible mass of a single cell. The dimensionless Stokes number Stk, expressing the ratio between the particle relaxation time and a time scale of the flow, quantifies the tendency of a solid particle to follow the streamlines of the flow, with Stk = 0 corresponding to the limit of perfect advection. In the present case, the particle Reynolds number (computed by using the particle diameter as the reference length) is as low as 0.5, confirming that cells dynamics develops in the Stokes regime. The Stokes number Stk is at most of the order of 10 −3 , lending full support to the masslessparticle assumption. 3. A one-way coupling is assumed between particles and the velocity field: the motion of the particles is determined by the fluid velocity field, but no feedback exists from the particles to the fluid, and no particle-particle interaction is considered. Neglecting the complex interactions among algae cells and between them and the fluid is justified by the minute dimensions of the cells (7 μm for Chlorella), whose presence is unable to directly affect the flow field. Unless very high-density cultures are considered, direct cell-cell interactions are also rare events that are neglected here.

Conclusions
The work has considered a thin-layer cascade (TLC) photobioreactor, which reportedly offers increased productivity in comparison to more traditional types of reactors, and has focused on the light/dark (L/D) cycle which the algae cells undergo during their turbulent motion in the liquid culture.
The goal was to establish whether the time scale associated to L/D can be related to the degree of unsteadiness in the light supply that is associated with an increase of photosynthetic efficiency.
To this aim, we have studied the numerical replica of a recent laboratory work where a TLC setup was designed and closely monitored to provide accurate control of external conditions and detailed information about several process variables. Our numerical study is based upon a numerical technique, the Direct Numerical Simulation (DNS) of the Navier-Stokes partial differential equations governing the motion of fluids, which to our knowledge has never been applied to such problem, and which is essential to provide temporal information. The significant computational cost of the DNS is counterbalanced by the absence of any type of turbulence modelling and by the availability of fully resolved temporal information.
A large number of particles have been placed in the flow to mimic the algae cells, and the trajectories of every particle have been studied to provide a statistical description of the L/D cycle. Our results indicate that, in the considered TLC reactor, the L/D cycle involves an average time scale ranging in the interval 0.1-2 s, and the single cell may actually experience significantly faster or slower L/D cycles. The PDF of the process has been fully characterised, since over one million particles could be considered.
The identified time scale is fast enough to support the hypothesis that turbulent mixing and the consequent unsteady illumination of every algae cell might be involved in the increased productivity of the thin-layer cascade. Our computational approach and in particular the statistical characterisation of the L/D cycle can be useful to further our understanding of the TLC and to provide design criteria aimed at maximising algal productivity. However, it should be noted that the actual productivity depends on several factors beyond the L/D cycle, and that some of them, as the nutrient uptake, are affected by the turbulent motions.
Our study has concerned a specific TLC setup. Although extremely expensive, a parametric study to systematically consider how changes flow rate and channel slope would affect the L/D characteristics would certainly be most valuable. Increasing either parameter leads to an increase in the Reynolds number, and the frequency of the L/D cycle is expected to increase as the intensity of the turbulent fluctuations increases and their time scale shrinks. However, increasing the Reynolds number above a certain threshold may negatively impact the algal culture. Indeed, the resulting shear stress may become excessive and cause increased cell mortality and decreased growth rate, as discussed for example by Wang and Lan (2018). Moreover, increasing the Reynolds number leads to a rapid surge of the computational demand, which for a DNS is estimated to be proportional to Re 3 .
As a concluding comment, we mention that the present numerical procedure can be straightforwardly extended to deal with other types of PBR, as for example Taylor-Couette, regardless of the physical mechanism driving the mixing in the culture.
Funding Open access funding provided by Politecnico di Milano within the CRUI-CARE Agreement. This work has been funded by Nord Energia S.p.A. Computing time has been provided by the Italian supercomputing centre CINECA under the ISCRA C project ALGAE-TL.

Compliance with ethical standards
Conflict of interest The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creativecommonshorg/licenses/by/4.0/.