Lluis F Marsal

Gael Gautier

Stéphane Bastide
email: bastide@icmpe.cnrs.fr

Elisa Pinna

Sylvain Le Gall

Encarnacion Torralba

Guido Mula

Christine Cachet-Vivier

Mesopore Formation and Silicon Surface Nanostructuration by Metal-Assisted Chemical Etching With Silver Nanoparticles

Keywords: MACE, silver nanoparticle, black silicon, reflectivity, texturization, nanostructuration

This article presents a study on Metal-Assisted Chemical Etching (MACE) of silicon in HF-H 2 O 2 using silver nanoparticles as catalysts. Our aim is a better understanding of the process to elaborate new 3D submicrometric surface structures useful for light management. We investigated MACE over the whole range of silicon doping, i.e., p ++ , p + , p, p -, n, n + , and n ++ . We discovered that, instead of the well-defined and straight mesopores obtained in p and n-type silicon, in p ++ and n ++ silicon MACE leads to the formation of cone-shaped macropores filled with porous silicon. We account for the transition between these two pore-formation regimes (straight and cone-shaped pores) by modeling (at equilibrium and under polarization) the Ag/Si/electrolyte (HF) system. The model simulates the system as two nanodiodes in series. We show that delocalized MACE is explained by a large tunnel current contribution for the p-Si/Ag and n-Si/HF diodes under reverse polarization, which increases with the doping level and when the size of the nanocontacts (Ag, HF) decreases. By analogy with the results obtained on heavily doped silicon, we finally present a method to form size-controlled cone-shaped macropores in p silicon with silver nanoparticles. This shape, instead of the usual straight mesopores, is obtained by applying an external anodic polarization during MACE. Two methods are shown to be effective for the control of the macropore cone angle: one by adjusting the potential applied during MACE, the other by changing the H 2 O 2 concentration. Under appropriate etching conditions, the obtained macropores exhibit optical properties (reflectivity ∼3 %) similar to that of black silicon.

INTRODUCTION

Metal-assisted chemical etching (MACE) of silicon is a powerful technique to produce surface nanostructures with high aspect ratios. Networks of nanowires, nanopores or nanocones and more complex patterns synthesized in this way find applications in various fields like microelectronics, solar energy conversion or chemical/biological sensing.

Recently, the texturing at the submicrometer scale of silicon into so called "black silicon" has been the object of a renewed interest for crystalline silicon solar cells [START_REF] Otto | Black silicon photovoltaics[END_REF]. This is largely due to advances in surface passivation, notably with the optimization of the emitter doping to avoid Auger recombination [START_REF] Oh | An 18.2%-efficient black-silicon solar cell achieved through control of carrier recombination in nanostructures[END_REF] and the use of thin passivating dielectric layers (e.g., Al 2 O 3 ) to overcome surface recombination involved with large surface areas [START_REF] Savin | Black silicon solar cells with interdigitated back-contacts achieve 22.1% efficiency[END_REF]. Thus, "black silicon" is being considered for practical photovoltaic applications since it offers a much lower reflectivity than conventional surfaces (e.g., pyramids with sizes of several micrometers), without the need of antireflection coating, and is also very efficient for light trapping in the cell. However, an adequate three-dimensional (3D) control of the texture is essential to achieve the required efficiency of the cells.

Within this framework, MACE has been shown to be a method of choice to produce efficient nanostructures for light management in silicon solar cells [START_REF] Koynov | Black nonreflecting silicon surfaces for solar cells[END_REF]. It continues to be studied as it is a relatively easy chemical method to implement and uses simple reagents. Its potential has not yet been fully explored, far from it, as many parameters have a relevant impact on the etching process. As a result, current research is devoted both to understanding the mechanisms of MACE and to determining etching conditions and modalities to fabricate tailored surface structures for efficient light coupling.

The nature of the metal is of primary importance in MACE. Silver nanoparticles are known to dig well-defined straight mesopores whose walls are smooth and with diameters set by the nanoparticle sizes [START_REF] Tsujino | Boring deep cylindrical nanoholes in silicon using silver nanoparticles as a catalyst[END_REF]. Silver is therefore used for highly localized etching and the obtention of well-resolved nanostructures. On the contrary, gold or platinum nanoparticles give rise to the formation of mesopores surrounded by a cone-shaped volume of porous silicon [START_REF] Lee | Pore formation in silicon by wet etching using micrometre-sized metal particles as catalysts[END_REF]. The reason for this difference is the nature of the Schottky metal/Si junction, which can be rectifying (silver) or ohmic (gold, platinum). In the first case (i.e., rectifying junction), the injection of holes is not possible while in the second it causes the polarization of the bulk and thus leads to the formation of porous silicon at the mesopore Si/electrolyte (HF) interface [START_REF] Torralba | Tunable surface structuration of silicon by metal assisted chemical etching with Pt nanoparticles under electrochemical bias[END_REF].

Other factors are also important in MACE, such as the doping of the silicon substrate, the metal shape (e.g., nanoparticle or mesh) or the composition of the etching solution. Under certain etching conditions, silver nanoparticles or silver-meshes have been shown to lead to the formation of porous silicon, i.e., to a delocalized rather than localized etching, in contrast to what is usually observed [START_REF] Chartier | Metal-assisted chemical etching of silicon in HF-H 2 O 2[END_REF][START_REF] Zhang | Preparation of large-area uniform silicon nanowires arrays through metal-assisted chemical etching[END_REF][START_REF] Geyer | Ag-mediated charge transport during metal-assisted chemical etching of silicon nanowires[END_REF]. The influence of substrate doping on the formation of mesopores by MACE with silver nanoparticles has therefore not been fully established yet.

In this work, we have studied the MACE process in view of an improved understanding for the control of the surface structures 3D shapes. In particular, the formation of mesopores in silicon with silver nanoparticles in HF/H 2 O 2 was studied on a wide range of substrate dopings, namely p ++ , p + , p, p -, n, n + , and n ++ . Our most significant discovery is that, instead of the well-defined nanometer-sized mesopores with smooth walls obtained in p-type and n-type silicon, a MACE process based on silver nanoparticles in highly doped p-and n-type silicon leads to the formation of micrometer-sized cone-shaped pores filled with porous silicon. We try to account for the transition between these two pore formation regimes by modeling the Ag/Si/electrolyte interface at the nanoscale, at equilibrium and under etching conditions (electrostatic polarization).

Finally, we present results on the effect of adding an external polarization during MACE of p-type silicon with silver nanoparticles as a method to form cone-shaped macropores rather than the usual straight mesopores. We test how the cone angle can be controlled by tuning the applied potential or the concentration of H 2 O 2 and the effect it has on the surface reflectivity.

MATERIALS AND METHODS

Silicon Substrates

Polished crystalline (100) p-type (boron) and n-type (phosphorus) silicon wafers (300-500 µm in thickness) from Sil'tronix and ITME were used, with three resistivities for each doping: n ++ (∼2.5 10 -3 Ω cm), n + (∼3.0 10 -2 Ω cm), n (∼2.0 Ω cm), p (∼2.0 Ω cm), p + (8.0 10 -2 Ω cm), p ++ (∼3.5 10 -3 Ω cm), and in addition p -doped silicon with a resistivity of ∼10 Ω cm.

Chemicals

Analytical grade (VWR chemicals) 30 wt.% H 2 O 2 , 96 wt.% H 2 SO 4 , 40 wt.% HF, 65 wt.% HNO 3 and ultra-pure water (18.2 MΩ cm, Millipore) were used in the all experiments. The AgNO 3 salt (VWR) was 99 wt.%.

Silver Deposition

Before silver deposition and etching, the silicon wafers were cleaved into 1.2 × 1.2 cm 2 pieces, cleaned in H 2 SO 4 -H 2 O 2 (3:1) and rinsed with ultra-pure water.

In order to obtain silver nanoparticles at the surface of the silicon substrate, a 10 nm thick silver layer was deposited using a Cressinton 208HR sputter coater and subsequently annealed 20 min at 275 • C under Ar atmosphere. The nanoparticle size distribution was established from SEM observations and ImageJ analysis (cf. Figure A in Supplementary Information). In some cases, silver nanoparticles were deposited chemically by dipping the silicon samples in a solution of AgNO 3 (1 mmol L -1 ) -HF (0.14 mol L -1 ) during 1 min (cf. Supplementary Information).

MACE and MACE Under Electrochemical Polarization

MACE was performed by dipping the samples in an aqueous solution of 1.21 mol L -1 HF and 0.21 mol L -1 H 2 O 2 for 20 min in all cases. The molar ratio ρ of the solution, i.e., [HF]

/([HF] + [H 2 O 2 ]), is equal to 0.85.
All MACE experiments under electrochemical polarization were performed in a home-made O-ring three-electrode cell. The silicon substrates (with or without silver nanoparticles) were used as working electrode and a Pt wire as counter electrode. A K 2 SO 4 saturated Hg/Hg 2 SO 4 electrode (SME), protected from the HF solution by a KCl saturated agar-agar bridge, was used as reference. The geometrical area of the Si/electrolyte contact was 0.38 cm 2 (delimited by an O-ring). The silicon back contact was

Instrumentation

Scanning electron microscopy (SEM) images and Energy Dispersive X-ray Spectroscopy (EDS) were obtained with a Merlin FEG microscope from Zeiss equipped with AZtec systems (EDS Advanced, HKL Advanced Nordlys Nano, Oxford Instruments).

Cyclic voltammetry and chronoamperometry were performed with a PGSTAT 20 potentiostat/galvanostat Metrohm Autolab, equipped with Nova software.

The optical reflectivity measurements were performed in the 350-1500 nm range using a Hitachi UV-VIS-NIR 4001 spectrophotometer fitted with an integrating sphere.

Modeling

Numerical simulations in 2D of the valence and conduction band modulation at the Ag/Si/electrolyte interfaces were performed using the commercial TCAD software (Atlas from Silvaco, cf. [START_REF] Torralba | Tunable surface structuration of silicon by metal assisted chemical etching with Pt nanoparticles under electrochemical bias[END_REF] based on a finite volume method. This simulator solves the physical equations governing the electrostatics (Poisson, electro-neutrality) and the transport of e - and h + (drift-diffusion) self-consistently on a 2D mesh.

The modeled structure is schemed in Figure B of Supplementary Information. It consists of a silicon substrate with a thickness of 100 µm and a width of 100 nm, with a 12 nm large silver pad surrounded by two electrolyte contacts. The electrolyte contacts are short-circuited (i.e., at same potential). The silver and electrolyte electrodes are separated by 1 nm of insulating vacuum to avoid charge transfer between them.

The work functions of silicon are taken at 4.07 eV < W Si < 4.29 eV depending on the doping level, W Ag = 4.64 eV [START_REF] Hölzl | Solid Surface Physics[END_REF], and W El = 4.5 eV (determined in our experimental conditions, cf. [START_REF] Torralba | Tunable surface structuration of silicon by metal assisted chemical etching with Pt nanoparticles under electrochemical bias[END_REF] and its Supplementary Information). The Fermi level is set at 0 eV at the equilibrium. To mimic the MACE process, a positive polarization can be applied between the silver contact and the electrolyte.

RESULTS AND DISCUSSION

Pore Morphology

Figure 1 presents SEM images at the same magnification of mesopores obtained under identical MACE conditions of silicon substrates with doping levels ranging from p ++ to n ++ .

Figure 2 shows the pore tips (left) and tops (right) at higher SEM magnification.

In the case of p, p -, n and n + -type silicon, the tip sections exhibit a constant diameter and pore walls free of porous silicon. However, at the level of the surface, they all exhibit a certain widening and rough inner surfaces. Most of the pores are perpendicular to the surface but a significant number is also found parallel to the surface [START_REF] Tsujino | Morphology of nanoholes formed in silicon by wet etching in solutions containing HF and H 2 O 2 at different concentrations using silver nanoparticles as catalysts[END_REF] have shown that silver nanoparticles create lateral pores because they can follow other preferential orientations equivalent to [100] [e.g., (010), (001). . . ] or when they have non-spherical shapes (pores are then dug in random or twisted directions). As a result, many silver nanoparticles are found close to the surface even after a long etching time. In the case of highly doped p ++ and n ++ -type silicon, pores with a surrounding cone-shaped mesoporous silicon volume are formed. For p + -type silicon, mesoporous silicon is formed around the mesopores as well, but to a much lower extent. The pore length decreases strongly as the doping level increases, with a ratio of ∼ 6-7 for p vs. p ++ and n vs. n ++ . FIGURE 2 | SEM images of mesopore tips (left) and tops (right), for the different substrate dopings. Same magnification for all images of each group. Note that the presence or absence of a nanoparticle at a pore tip is only a random result of cutting the sample for cross-sectional SEM imaging (depending on whether it remains stuck to the observed piece). Scale bar: 200 nm.

These results are specific to the HF and H 2 O 2 concentrations chosen for etching. For other HF and H 2 O 2 concentrations, the frontiers in terms of doping level between cone-shaped pores (i.e., with mesoporous silicon) and straight pores formation may change. However, MACE with silver nanoparticles in low-doped (n-and p-type) silicon always leads to straight pores, irrespective of the reactant concentrations (for 0.7 < ρ < 1); this is reported for instance in Chartier et al. for p-type silicon (1-2 ohm.cm) in 14 mol L -1 HF and 3 mol L -1 H 2 O 2 (i.e., ten times higher concentrations than in the present case).

In order to understand why and how the pore morphology changes with the type and level of doping, a modeling describing the silicon bands modulation around nanometer sized silver and electrolyte contacts has been carried out. This type of modeling has already provided insights into the photo/electrochemical behavior of metal decorated silicon electrodes [START_REF] Nakato | Effect of microscopic discontinuity of metal overlayers on the photovoltages in metal-coated semiconductor-liquid junction photoelectrochemical cells for efficient solar energy conversion[END_REF], on MACE [START_REF] Kolasinski | Electron transfer during metal-assisted and stain etching of silicon[END_REF] and electrochemically assisted MACE mechanisms [START_REF] Chourou | Metal-assisted etching of p-type silicon under anodic polarization in HF solution with and without H 2 O 2[END_REF][START_REF] Huang | Metal-assisted electrochemical etching of silicon[END_REF][START_REF] Torralba | Tunable surface structuration of silicon by metal assisted chemical etching with Pt nanoparticles under electrochemical bias[END_REF][START_REF] Bastide | 3D patterning of Si by contact etching with nanoporous metals[END_REF]. The results have also been interpreted in the light of data from the large literature existing on the electrochemistry of silicon in HF media, in particular the reference book by [START_REF] Lehmann | The Electrochemistry of Silicon: Instrumentation, Science, Materials and Applications[END_REF].

The common characteristic of MACE processes catalyzed by silver nanoparticles, irrespective of the silicon type and doping level, is the formation of a main pore having the dimensions and shape of the nanoparticle with nanometer-level accuracy. It is actually the basis for the formation of nanowires networks in AgNO 3 /HF solutions. This ultra-localized dissolution must correspond to a surface phenomenon where the silicon atoms located very close to the metal contact are subjected to a lateral electric field strong enough to attract their electrons directly to the metal.

The change in pore morphology is related to the additional presence of mesoporous silicon at the main pore walls. This can represent a thin layer of a few nanometers either along the pore (cf. p + silicon in Figure 2) or just at its apex (cf. n + silicon in Figure 2), or a cone of mesoporous silicon much wider than the main pore (cf. p ++ and n ++ silicon in Figure 2). This constitutes a delocalized dissolution implying a polarization of the Si/HF interface far away from the silver nanoparticle which is only possible through a polarization of the silicon bulk. It is this case that we can try to account the observed phenomena by a modeling of the process and from data in the literature.

Band Bending Modeling

The modeled Ag/Si/electrolyte (HF) system consists in a silver pad on flat silicon surrounded by two electrolytic contacts (short-circuited) with a total width of 100 nm (cf. Figure B of Supplementary Information for details on the modeled device structure). The silver and electrolyte electrodes are separated by a gap of 1 nm, hence charge transfers only occurs through the Si/electrolyte (HF) and Si/Ag interfaces.

Figure 3 presents the valence and conduction band diagrams for p-type and and n-type silicon at equilibrium (top row). They correspond to the band modulations along a cutline (ycoordinate) that connects the center of the silver nanoparticle to the electrolyte 50 nm away along the x-axis, within bulk silicon over a depth of 1 µm. A scheme of the device is given on the left side (diagrams for p-type silicon at equilibrium) and in more detail in Figure B of the Supplementary Information. Three doping level are represented: moderate (p, n; doping of 3 × 10 15 cm -3 ), high (p + , n + ; 3 × 10 17 cm -3 ) and very high doping level (p ++ , n ++ ; 1 × 10 19 cm -3 ). In all cases, E F is set at 0 eV. The two Schottky diodes Ag/Si and Si/HF are back-to-back. As expected, the space-charge region (SCR) decreases as the doping level increases. At the Ag/p-Si contact, the valence band describes a near-surface peak that may appear unusual in a band diagram. This is actually due to the influence of the Electrolyte/Si junction which has a more pronounced band bending (at equilibrium). It results in a lateral modulation of the band even under the silver contact (at 6 nm from the edge surrounding the contact). This illustrates the advantage of 2D modeling to reveal such lateral effects that cannot be depicted from simplified 1D band diagrams.

To mimic MACE, the oxidizing power of H 2 O 2 is accounted for by a positive polarization applied to silver vs. electrolyte, as shown in Figure 3 (bottom row) in the case of moderate doping (p, n). The potential drop is distributed spatially between these two junctions according to their electronic characteristics. With n-type silicon, the Ag/n-Si diode is under forward bias while the n-Si/HF diode is blocking. Most of the potential drop occurs at the Ag/Si interface. In the case of p-type silicon, the Ag/p-Si diode is blocking while the p-Si/HF diode is under forward bias and most of the potential drop occurs at the Ag/Si interface as well. Therefore, the systematic presence of a blocking diode should not allow delocalized MACE, i.e., electron transfer from surface silicon atoms (oxidized and dissolved in HF) to H 2 O 2 molecules (reduced on silver).

Results from the literature show, however, that different silicon doping conditions and enhanced electric field due to geometrical effects (contact size, pore tip) can radically alter the charge transfer for this type of diode under reverse polarization [START_REF] Lehmann | The Electrochemistry of Silicon: Instrumentation, Science, Materials and Applications[END_REF][START_REF] Smit | Scaling of nano-Schottkydiodes[END_REF][START_REF] Smit | Conductance distribution in nanometer-sized semiconductor devices due to dopant statistics[END_REF][START_REF] Donolato | Approximate analytical solution to the space charge problem in nanosized Schottky diodes[END_REF][START_REF] Vostokov | On the role of tunneling in metal-semiconductor nanocontacts[END_REF].

N-Type Si/HF Diode

While for n-type silicon under anodic (reverse) polarization in the dark, the absence of h + in the valence band prevents porous silicon formation, a mesoporous silicon layer is formed for n + and n ++ silicon. The anodic current is actually produced by tunneling of electrons through the SCR from silicon surface atoms to the conduction band. This occurs at relatively low potentials (especially for n ++ ) because the electric field at a FIGURE 4 | I-V characteristics of n and n + silicon in contact with an HF electrolyte (W = 4.5 eV). The electrolyte contact diameter is 100 nm (A) and 12 nm (B). The current is calculated either for thermionic emission (including image potential effect), and with or without tunneling current contribution. mesopore tip is much stronger (due to the radius of curvature) than at a flat surface (Lehmann, 2002, chapter 8). This effect becomes significant when the radius of curvature of the pore tip is smaller than the width of the SCR. Regarding MACE, to the best of our knowledge, n-type silicon etching via the conduction band has only been proposed by [START_REF] Yae | High catalytic activity of palladium for metal-enhanced HF etching of silicon[END_REF], in the particular case of palladium particles in HF medium without oxidizing agent. They demonstrated that the oxidation of surface silicon atoms was coupled to the reduction of water on palladium by electron injection into the conduction band.

I-V characteristics established with our n-type silicon samples in HF at the concentration used for MACE (cf. Figure D in Supplementary Information), clearly show a lack of anodic current for n-type silicon, a small current for n + -type silicon and a significant current for n ++ -type silicon. A porous silicon layer is rapidly built at the surface in the latter case.

We investigated whether it was possible to retrieve these results by establishing the I-V characteristics from electrostatic modeling (Silvaco) of an electrolytic contact (HF) on heavily, moderately and lightly doped silicon. To evaluate the size effect, circular silver pads with a diameter of 100 and 12 nm were tested (cf. Figure C of Supplementary Information). Figure 4 plots the I-V characteristics so obtained. The current is calculated for thermionic emission only (including image potential effect) or with taking into account tunneling transport through the SCR (field emission and thermionic field emission [START_REF] Sze | Physics of Semiconductor Devices[END_REF][START_REF] Rhoderick | Metal-semiconductor contacts[END_REF].

With an electrolytic contact of 100 nm, the contribution of tunneling to the reverse thermionic current is negligible for ntype silicon but becomes significant for n + -doped silicon (i.e., it increases by a factor of ∼3), as shown in Figure 4A. When the contact diameter is reduced from 100 to 12 nm, the reverse current density increases regardless of the doping level. This is due to a reduced SCR thickness [START_REF] Smit | Scaling of nano-Schottkydiodes[END_REF], but the (thermionic) forward current also increases (not shown in Figure 4B) which indicates a lowering of the barrier height as well [START_REF] Vostokov | On the role of tunneling in metal-semiconductor nanocontacts[END_REF].

The modeling thus makes it possible to account for the size effect in Schottky contacts (higher tunneling current, lower barrier height) that has been highlighted by [START_REF] Smit | Scaling of nano-Schottkydiodes[END_REF], [START_REF] Donolato | Approximate analytical solution to the space charge problem in nanosized Schottky diodes[END_REF], and [START_REF] Hugelmann | Schottky diode characteristics of electrodeposited Au/n-Si(111) nanocontacts[END_REF].

P-Type Si/Ag Diode

Under anodic polarization, the p-Si/Ag diode is blocking. Figure 5 shows the I-V characteristics of p-Si/Ag diodes established by modeling, for two contact sizes (100 and 12 nm) and three doping levels (p, p + and p ++ ).

The reverse-biased Si/Ag solid junction (V Ag-Si > 0) produces a non-negligible current if tunneling is taken into account, even when the doping level is low (p) and the contact size large (100 nm). The current density at +0.5 V is 220 and 360 times higher than that calculated for thermionic emission alone, in the case of p and p + -doped silicon, respectively. As a result, the diodes are less rectifying as the silicon doping level increases.

The size effect (reduction in diameter from 100 to 12 nm) is not pronounced for p and p ++ -doped silicon but important for p + -doped silicon. In this latter case, the reverse current density (tunneling) is almost equal to the forward current density (thermionic emission) for the 12 nm silver contact. Therefore, both the increase in doping level and the downsizing of Schottky diodes to nanometer sized silver contacts lead to significant tunneling currents.

Summary of the Modeling Results

To sum up, it appears that delocalized MACE is not possible with lightly doped p-or n-type silicon since in each case one of the two (back to back) diodes is blocking. Some porous silicon is actually visible at the very top of the pores though (cf. Figure 2), probably as a result of a very small current accumulated during the 20 min of etching (the surface experiences the longest etching time). This phenomenon is accentuated by the presence of numerous silver nanoparticles near the surface that etch the substrate laterally (see the large number of pore sections visible in the cross-section images of Figure 2), contributing to a higher leakage current density in this region and thus to some porosification of the surface. In addition, silver oxidation by H 2 O 2 with diffusion of Ag + in solution far from the initial location can also lead to non-local etching, as studied in detail by [START_REF] Chiappini | Biodegradable porous silicon barcode nanowires with defined geometry[END_REF] and [START_REF] Geyer | Ag-mediated charge transport during metal-assisted chemical etching of silicon nanowires[END_REF], and this phenomenon would be more important near the surface because of the high number of silver nanoparticles there.

On the contrary, MACE can be delocalized with heavily doped silicon (n ++ , p ++ ), for which the reverse currents are significant. In both cases, only majority carriers are involved. With highly doped silicon (n + , p + ), the situation is intermediate for the MACE conditions used in this work. With p + silicon, we observe the presence of porous silicon all along the main pores, whereas for n + doped silicon, porous silicon only appears at the top of main pores (i.e., for long exposure times to HF), which implies that the reverse current at n + -Si/HF exists but is low.

Control of the Pore Morphology

The experimental and modeling results discussed above acquire a special interest in the fabrication of submicrometric silicon surface textures, whose 3D aspect is a key to control the samples behavior when interacting with light, e.g., in the case of solar cells or photodetectors.

This type of texturization has shown its interest in solar cells with the so-called "black silicon." Black silicon can be obtained by various techniques, e.g., plasma etching, laser or chemical treatments [START_REF] Otto | Black silicon photovoltaics[END_REF]. Some studies have shown that MACE with silver nanoparticles followed by a light alkaline etching allows the formation of black silicon for very high efficiency solar cells [START_REF] Oh | An 18.2%-efficient black-silicon solar cell achieved through control of carrier recombination in nanostructures[END_REF]. In this case, the pores were limited in depth and width to ∼ 250 and 100 nm, respectively, without a well-defined morphology. More generally, and as we have shown, MACE with silver nanoparticles produces straight mesopores over a range of doping going from p + to n + . An increase in HF and H 2 O 2 concentrations or in the relative fraction of H 2 O 2 vs. HF results in straight mesopores as well [START_REF] Chartier | Metal-assisted chemical etching of silicon in HF-H 2 O 2[END_REF], hence the chemical control of the nanostructure morphology is not possible with silver catalysts.

On the other hand, MACE with other metals gives different results. Gold nanoparticles have also been investigated for MACE applied to solar cell processing [START_REF] Koynov | Black nonreflecting silicon surfaces for solar cells[END_REF][START_REF] Algasinger | Improved black silicon for photovoltaic applications[END_REF][START_REF] Otto | Black silicon photovoltaics[END_REF]. Silicon nanocones are obtained (∼ 500 nm in height and 250 nm in width). The resulting optical properties are very good, with a low effective reflectivity (∼1%) and efficient light trapping. However, no means of morphology control was reported. MACE with platinum nanoparticles leads to the formation of cone-shaped macropores (after dissolution of porous silicon) with efficient light coupling properties as well [START_REF] Torralba | Tunable surface structuration of silicon by metal assisted chemical etching with Pt nanoparticles under electrochemical bias[END_REF]. Obtaining these structures rather than straight pores with silver is linked to the nature of the Si/Metal contact, ohmic with gold and platinum, Schottky with silver (for silicon doped between p + and n + ).

However, if the objective is to precisely control the morphology of the etched structures, for instance for the design of specifically texturized devices [START_REF] Pinna | Electrochemical nanolithography on silicon: an easy and scalable method to control pore formation at the nanoscale[END_REF], even with these metals MACE needs to play on additional parameters. In the case of platinum, it was necessary to add an electrochemical polarization in order to adjust the rate of porous silicon formation around the mesopore etched by the nanoparticles and hence to control the opening of the coneshaped macropores [START_REF] Torralba | Tunable surface structuration of silicon by metal assisted chemical etching with Pt nanoparticles under electrochemical bias[END_REF]. From a practical point of view, the disadvantage of using gold or platinum as catalysts is their high cost and the need to use a strong oxidant (aqua regia) to ultimately remove the nanoparticles (at the bottom of the pores) to avoid subsequent contamination or formation of gold/platinum silicide during high temperature treatments.

Taking these problems into account, it would be advantageous to be able to use silver rather than gold or platinum in electrochemically assisted MACE since it is both much easier to dissolve and has a lower cost. For these reasons, a study on MACE of p-type silicon with silver nanoparticles under anodic polarization of the substrate has been carried out. Figure 6 (top line) shows the result of the etching under the same conditions as those of Figures 1, 2 but with the addition of an increasing anodic polarization ranging from +0.1 to +0.5 V. The OCP measured under MACE conditions being of the order of-0.5 V vs. SME, this represents potentials ranging from -0.4 to 0 V vs. SME.

In all cases, a porous silicon corolla is formed around the main mesopores, the amount of which increases with the potential. This is expected for p-type silicon under forward bias in HF medium. The anodic current density increases from 1.0 mA cm -2 to 6.1 mA cm -2 at -0.4 and 0 V vs. SME, respectively.

The SEM images of the middle and bottom lines in It should be noted that at -0.2 V vs. SME the angle shown in Figure 7A is 75 • (with a relatively large standard deviation) whereas in the corresponding SEM image in Figure 6, very open cones are also observed. Due to the density of the silver nanoparticles, there is a large overlap between the macropores and only those which are by chance isolated develop a porous silicon cone in accordance with the imposed anodic current. In the vast majority of cases, only the base of the macropores remains due to the superposition of the porous silicon cones.

The whole approach is summarized by the diagrams in Figure 8, which brings together the results obtained (on ptype silicon) by electrochemical anodization in HF medium (V A > OCP), by conventional MACE (OCP, V A = 0) and eventually by the two processes carried out simultaneously (MACE at V A > OCP).

To control the macropore dimensions, the applied potential can be adjusted. More opened macropores can be obtained by increasing the anodic potential, as shown for V 2 > V 1 .

We have also investigated another option that consist in changing the H 2 O 2 concentration since it controls the penetration rate of the silver nanoparticles. For that purpose, we tested the effect of doubling the H 2 O 2 concentration at constant HF concentration. The ρ-value is slightly modified (from 0.85 to 0.75) and the penetration rate is increased by ∼50 %. The open circuit potential shifts from -0.50 V to -0.46 V vs. SME when the H 2 O 2 concentration is increased from 0.22 to 0.42 mol L -1 . One sample was treated with MACE under anodic polarization at -0.37 V vs. SME (+0.1 V vs. OCP) and another at -0.26 V vs. SME (+0.2 V vs. OCP). Figure E in Supplementary Information gives the chronopotentiometry of these sample during etching. Figure 9 shows SEM images of the sample surface after these treatments and subsequent etching of the porous silicon layer.

The difference is striking at low applied potential (-0.37 V vs. SME), with cone shaped macropores much more developed and narrower. At -0.26 V vs. SME, the macropore depth with respect to the porous silicon surface is ∼1.95 µm (after 20 min).

The angle values of the cone-shaped macropores thus obtained are reported in Figure 7A (red dashed line and dots). Compared to those reported for a lower H 2 O 2 concentration (blue line and dots), it is clear that, at constant potential, the macropores are ∼16-18% narrower. This evolution is logical if we consider that the porous silicon cone shape results from the combination of two perpendicular etchings: (i) the etching normal to the surface of the nanoparticle penetrating the substrate (mesopore) and (ii) the formation of porous silicon normal to this mesopore (i.e., perpendicular to the mesopore wall). This is schemed in Figure 7B. Assuming that the anodic current depends only on the HF concentration, the raise of H 2 O 2 only increases the penetration rate of nanoparticles (i.e., the mesopore depth). For example, for an applied potential of -0.26 V vs. SME, the cone angle increases from 40 • to 55 • when the H 2 O 2 concentration is doubled, which would correspond to an increase in the penetration rate of the nanoparticles by a factor of 1.4, in relative agreement with the factor 1.6 measured from the SEM images of Figure 9 (1.9 µm/20 min at -0.26 V vs. SME) and Figure 6 (1.2 µm/20 min at -0.3 V vs. SME). Note that we exclude an influence of the anodic current on the penetration rate of silver nanoparticles due to a preferential collection of h + at the mesopore tips. This remarkable effect has been reported in the literature [START_REF] Chourou | Metal-assisted etching of p-type silicon under anodic polarization in HF solution with and without H 2 O 2[END_REF][START_REF] Huang | Metal-assisted electrochemical etching of silicon[END_REF], but it concerns the electrochemistry in HF medium only (i.e., without H 2 O 2 ). In the presence of H 2 O 2 , the lateral curvature of the bands is no longer favorable to the diffusion of h + from p-Si/electrolyte to p-Si/Ag areas.

The surface reflectivity was measured from 350 to 650 nm for textured surfaces with cone-shaped macropores of different opening angles, as shown in Figure 10, for a low (A) and high (B) concentration in H 2 O 2 .

The reflectivity spectra in Figure 10A are found to be similar to that of polished silicon, indicating that there is no optical interference due to the presence of an intermediate index layer, as would be the case with remains of porous silicon. The evolution of the reflectivity at 550 nm with the cone angle is given in the inset of Figure 10A. For the sample treated at -0.2 V vs. SME, a high reflectivity is measured, that corresponds to shallow macropores with large openings (cf. Figure 6). The least reflective surface (9 % at 550 nm) is obtained for cone-shaped macropores with an opening angle of ∼ 49 • (cf. Figure 6, -0.3 V vs. SME), which is comparable to the reflectivity of (100) oriented silicon with inverted square-based pyramids [START_REF] Magnin | Angle-dependent ray tracing simulations of reflections on pyramidal textures for silicon solar cells[END_REF] obtained by lithography and alkaline etching for high efficiency solar cells. At -0.4 V vs. SME, the macropores are smaller and narrower and the reflectivity slightly higher (13 %).

These results can be rationalized from the comprehensive analysis of the optical properties of submicrometer structures (black silicon) developed by [START_REF] Otto | Black silicon photovoltaics[END_REF]. The light trapping performances depend mainly on the correlation length L c , which corresponds approximatively to the lateral distance between adjacent peaks and valleys and influences the fraction of scattered light. The larger the correlation length, the better the light trapping. At the same time, the reflectivity also increases with L c due to backscattering. The second important parameter is the peak to valley height H of the structures, which is detrimental to the antireflection properties and light trapping if it is less than ∼500 nm, since the refractive index gradient is not sufficient to couple light efficiently.

At low H 2 O 2 concentration, macropores obtained at -0.4, -0.3 and -0.2 V vs. SME exhibit: L c = 0.23, 0.32, 0.52 µm, and H = 0.28, 0.52 µm, 0.28 µm, respectively. Therefore, the lowest reflectivity is obtained at-0.3 V because the macropores are narrower than those obtained at -0.2 V vs. SME and much deeper than those obtained at -0.4 V vs. SME.

At high H 2 O 2 concentration, the macropores get deeper and narrower, which translates for etching at -0.37 V and-0.26 V in: L c = 0.23 µm and 0.34 µm; H = 1.49 and 0.87 µm, respectively. The reflectivity spectra of the sample etched at -0.37 V vs. SME does not vary monotonically like the others but exhibit a minimum (5.9 %) around 410 nm (cf. Figure 10B). Despite suitable L c and H values, the reflectivity is not that low, probably because the refractive index profile in the macroporous layer is too steep (cf. Figure 9A). The macropores obtained at-0.26 V vs. SME exhibit cone angle of ∼ 40 • , and a reflectivity of 3% at 550 nm (H = 0.87 µm, L c = 0.34 µm). This number is comparable to some values reported in the literature for black silicon [START_REF] Oh | An 18.2%-efficient black-silicon solar cell achieved through control of carrier recombination in nanostructures[END_REF][START_REF] Hirsch | Optoelectronic properties of Black-Silicon generated through inductively coupled plasma (ICP) processing for crystalline silicon solar cells[END_REF][START_REF] Li | Hybrid black silicon solar cells textured with the interplay of copper-induced galvanic displacement[END_REF]. However, black silicon obtained by MACE with Au nanoparticles can lead to effective reflectivity as low as ∼ 0.6-1 % under certain etching conditions, as reported by [START_REF] Algasinger | Improved black silicon for photovoltaic applications[END_REF][START_REF] Otto | Black silicon photovoltaics[END_REF]. In this case, the difference lies in structures with a L c closer to the optimal value of 0.1 µm, together with a significant depth H of 0.65 µm. On the other hand, larger L c values, as obtained here by MACE with silver nanoparticles under anodic polarization, lead to strong light trapping. With this respect, the cone-shaped macropores should be as efficient as the best "black silicon" nanostructures obtained by photoelectrochemical etching that exhibit a L c value of 0.36 µm [START_REF] Otto | Black silicon photovoltaics[END_REF]. This can be advantageous in several cases, like for instance in ultrathin silicon solar cells.

CONCLUSION

In this work, we studied the mechanisms involved in MACE of silicon in HF-H 2 O 2 medium using silver nanoparticles with the aim to explain and control the process of 3D surface nanostructuring.

A systematic investigation over the whole range of silicon doping was conducted. We show for the first time that there is a delocalized MACE of p ++ and n ++ silicon with formation of a porous silicon cone around the main mesopores. Bibliographic data and 2D modeling have allowed to rationalize the Ag/Si/electrolyte (HF) system into two nanodiodes (Ag/Si and Si/HF) in series. During MACE, H 2 O 2 induces a positive polarization of silver with respect to the HF electrolyte, with one of the nanojunctions under reverse (blocking) polarization, either Ag/Si for p-type silicon or Si/HF for n-type silicon. Simulations of the I-V curves of these nanodiodes showed that the reverse current increases with the doping level and when the size of the metal contact (i.e., nanoparticle size) decreases, both factors leading to a more important contribution of tunneling current to the overall current (in the case of high dopings, the reverse and forward currents can even be of the same magnitude). Therefore, in the case of n-type silicon, and as already known for the electrochemical formation of mesoporous silicon in the dark, the delocalization of MACE can occur by injection of electrons from silicon surface atoms into the conduction band (majority carrier process). By analogy with these results, the formation of coneshaped macropores in p-type silicon was obtained with silver nanoparticles and the help of an external anodic polarization during MACE. The anodic current makes it possible to form a porous silicon cone around the mesopores etched by the nanoparticles. We have shown that the angle of the cone and the depth of the macropores can be controlled either by the anodic potential (i.e., the rate of porous silicon formation, parallel to the surface) or through the H 2 O 2 concentration which determines the penetration rate of the silver nanoparticles (perpendicular to the surface). Black silicon with reflectivity values of about 3% was obtained (after porous silicon removal) with a macropore depth of 0.8 µm and a correlation length of 0.34 µm. Gold or platinum nanoparticles are already known as MACE catalysts leading to the formation of nanocones or cone-shaped pores. As demonstrated here, silver allows a high level of control in nanostructuring but presents also the advantage of being less expensive and more easily removed, a crucial step to avoid contamination during subsequent processing.

FIGURE 1 |

 1 FIGURE 1 | SEM images of mesopore tips formed by MACE with silver nanoparticles and HF-H 2 O 2 in silicon substrates with doping levels: p ++ , p + , p, p -, n, n + , n ++ . Same magnification for all images. Scale bar: 1µm.

FIGURE 3 |

 3 FIGURE 3 | Band diagram (valence band and conduction band for p-type and n-type silicon, respectively) of the modeled Ag/Si/electrolyte (HF) system at equilibrium (top row) and under a positive polarization of silver vs. Electrolyte by +0.5 V (bottom row). x is the lateral position at the surface of the device and y the position along the cutline through bulk Si, as schematized in the upper-left band diagram and in Figure B of the Supplementary Information. Three doping levels for n and p-type silicon are represented (right and left column, respectively). The spatial distribution of the voltage drop is also presented (red dash-dot line).

FIGURE 5 |

 5 FIGURE 5 | I-V characteristics of p, p + and p ++ silicon in contact with silver (W = 4.64 eV). The silver contact diameter is 100 nm (A) and 12 nm (B). The current is calculated either for thermionic emission (including image potential effect), and with or without tunneling current contribution.

FIGURE 6 |

 6 FIGURE 6 | SEM images of cone-shaped pores formed in p-type silicon by MACE with silver nanoparticles in HF-H 2 O 2 with in addition an external anodic potential.Same magnification for all images. Top row: cross section views after etching. Middle and bottom rows: after etching and porous silicon removal in HF-HNO 3 in cross section and plane view, respectively. Scale bar: 500 nm.

  Figure 6 are obtained after dissolution of the porous silicon layer in HF-HNO 3 (1:99) in transverse and plan view, respectively. Cone-shaped macropores are then observed with an opening angle that increases as a function of the applied potential, as reported in Figure 7A (blue dashed line and squares, [H 2 O 2 ] = 0.22 mol L -1 ).

FIGURE 7 |

 7 FIGURE 7 | (A) Opening angle of the cone-shaped macropores obtained at different anodic potentials during MACE performed in 1.24 mol L -1 HF for two different H 2 O 2 concentrations corresponding to the experiments presented in Figure 6 (0.22 mol L -1 ) Figure 9 (0.42 mol L -1 ). The error bar represents the standard deviation for about ten measured macropores. The dashed lines are guides for the eyes. (B) Scheme of a cone-shaped pore formation by combined MACE and porous silicon anodic formation with etch rates V AgNP and V PSi , for low and high H 2 O 2 concentrations (blue and red half-cones, respectively).

FIGURE 8 |

 8 FIGURE 8 | Scheme of the electrochemical anodization and MACE of p-type silicon and their combination for two different applied potentials.

FIGURE 9 |

 9 FIGURE 9 | SEM images in cross section of cone shaped pores formed in p-type silicon by MACE with silver nanoparticles in HF-H 2 O 2 with in addition an external anodic polarization (A: -0.37 V vs. SME; B: -0.26 V vs. SME) before (left) and after (right) removing porous silicon in HF-HNO 3 . Same magnification for all images.

FIGURE 10 |

 10 FIGURE 10 | Surface reflectivity of p-type silicon covered with cone shaped macropores obtained by MACE under anodic polarization, followed by dissolution of porous silicon in HF/HNO 3 . MACE is performed in 1.24 mol L -1 [HF] and (A) 0.22 mol L -1 [H 2 O 2 ] or (B) 0.42 mol L -1 [H 2 O 2 ]. Anodic potentials applied during MACE are indicated in all cases. The inset of (A) gives the reflectivity at 550 nm as a function of the cone angle.

Frontiers in Chemistry | www.frontiersin.org

July 2020 | Volume 8 | Article 658

ACKNOWLEDGMENTS

The authors acknowledge the support of the Centre National de la Recherche Scientifique (CNRS) and of the Institut de Chimie (INC).

DATA AVAILABILITY STATEMENT

The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS

All authors participated in the design and planning of the research. EP and ET participated to the deposition of silver nanoparticles, MACE experiments and SEM analysis. SLG did the modeling and numerical simulations of the Ag/Si/HF system and the nanodiodes. ET and SB performed the electrochemical experiments. All authors participated to the data analysis discussion. SB wrote the first draft of the manuscript with SLG and ET providing feedback. GM and SB wrote the final version of the manuscript.

FUNDING

EP thanks the Regione Autonoma della Sardegna for its contribution to her PhD grant.

SUPPLEMENTARY MATERIAL

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem. 2020.00658/full#supplementary-material

Conflict of Interest:

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Copyright © 2020 Pinna, Le Gall, Torralba, Mula, Cachet-Vivier and Bastide. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.


Olivier Bernard
email: olivier.bernard@inria.fr

Liudi Lu
email: liudi.lu@inria.fr

Jacques Sainte Marie
email: jacques.sainte-marie@inria.fr

Julien Salomon
email: julien.salomon@inria.fr

Jacques Sainte

Shape optimization of a microalgal raceway to enhance productivity

Keywords: microalgae, hydrodynamic, optimization, han model, shallow water systems

We consider a coupled physical-biological model describing growth of microalgae in a raceway pond cultivation process, accounting for hydrodynamics. Our approach combines a biological model (based on the Han model) and shallow water dynamics equations that model the fluid into the raceway pond. We present an optimization procedure dealing with the topography to maximize the biomass production over one lap or multiple laps with a paddle wheel. On the contrary to a widespread belief in the microalgae field, the results show that a flat topography is optimal in a periodic regime. In other frameworks, non-trivial topographies can be obtained. We present some of them, e.g., when a mixing device is included in the model.

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 [START_REF] Eilers | Dynamic behaviour of a model for photosynthesis and photoinhibition[END_REF][START_REF] Masojídek | A closed solar photobioreactor for cultivation of microalgae under supra-high irradiance: Basic design and performance[END_REF][START_REF] Sung | Design of experiments and sensitivity analysis for microalgal bioreactor systems[END_REF][START_REF] Lamare | Gradientbased optimization of a rotating algal biofilm process[END_REF]. 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 [START_REF] René | An outlook on microalgal biofuels[END_REF]. 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 [START_REF] Chiaramonti | Review of energy balance in raceway ponds for microalgae cultivation: Re-thinking a traditional system is possible[END_REF] which homogenizes the medium for ensuring an equidistribution of the nutrients and guarantees that each cell will have regularly access to light [START_REF] Demory | How do microalgae perceive light in a high-rate pond? towards more realistic lagrangian experiments[END_REF]. The algae are periodically harvested, and their concentration is maintained around an optimal value [START_REF] Muñoz-Tamayo | Optimizing microalgal production in raceway systems[END_REF][START_REF] Posten | Microalgae Biotechnology[END_REF]. Light penetration is strongly reduced by the algal biomass, resulting in limited algal concentrations, which stay below 1% [START_REF] Bernard | Modelling of microalgae culture systems with applications to control and optimization[END_REF] 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 [START_REF] Masci | Microalgal biomass surface productivity optimization based on a photobioreactor model[END_REF][START_REF] Grognard | Optimal strategies for biomass productivity maximization in a photobioreactor using natural light[END_REF]. 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 [START_REF] Bernard | Controlling the bottom topography of a microalgal pond to optimize productivity[END_REF]. 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 [START_REF] Bernard | A 2d model for hydrodynamics and biology coupling applied to algae growth simulations[END_REF] by combining the Han photosynthesis equations with an hydrodynamic law based on the shallow water equations. 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.

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.

Modeling the photosystems dynamics

We consider the Han model [START_REF] Han | Photosynthesis-irradiance response at physiological level: A mechanistic model[END_REF] 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. 1.

A B C σI k d σI τ -1 k r
Photon I Photon I Their evolution satisfy the following dynamical system

       Ȧ = -σIA + B τ , Ḃ = σIA -B τ + k r C -k d σIB, Ċ = -k r C + k d σIB.
Here A, B and C are the relative frequencies of the three possible states with

A + B + C = 1, ( 1
)
and I is the photon flux density, a continuous time-varying signal. The other parameters are σ, that stands for the specific photon absorption, τ which is the turnover rate, k r which represents the photosystem repair rate and k d which is the damage rate. As shown in [START_REF] Lamare | Gradientbased optimization of a rotating algal biofilm process[END_REF], one can use (1), a fast-slow approximation and singular perturbation theory to reduce this system to a single evolution equation:

Ċ = -α(I)C + β(I), (2)
where

α(I) = β(I) + k r , with β(I) = k d τ (σI) 2 τ σI+1
. The net specific growth rate is obtained by balancing photosynthesis and respiration, which gives

µ(C, I) = -γ(I)C + ζ(I), (3)
where ζ(I) = γ(I) -R, with γ(I) = kσI τ σ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 µ:

Ẋ = µ(C, I)X -DX, (4)
where D is the dilution rate. The growth rate (3) associated with the steady state, for a constant I is given by

µ(I) = -γ(I)C * (I) + ζ(I), (5)
where

C * (I) = β(I) α(I) . ( 6
)

Steady 1D shallow water equations

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 [START_REF] Gerbeau | Derivation of viscous saint-venant system for laminar shallow water; numerical validation[END_REF]). 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:

∂ x (hu) = 0, (7)
∂ x (hu 2 + g h 2 2 ) = -gh∂ x z b , ( 8
)
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 η is given by η = h + z b and the averaged discharge Q = hu. This system is presented in Fig. 2. 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).

0 L z 0 x η(x) I s z b (x) h(x) u(x)
Integrating [START_REF] Eilers | Dynamic behaviour of a model for photosynthesis and photoinhibition[END_REF], we get

hu = Q 0 , (9)
for a fixed positive constant Q 0 , which implies a constant discharge in space. Then (8) can be rewritten as

hu∂ x u + h∂ x gh + h∂ x gz b = 0. ( 10
)
Let us assume that h is strictly positive. Dividing (10) by h and using [START_REF] Grenier | Worldwide perspectives of rotating algal biofilm up-scaling[END_REF] to eliminate u, we get

∂ x Q 2 0 2h 2 + g(h + z b ) = 0.
This equation corresponds to the Bernoulli's principle. Now let us consider two fixed constants h(0), z b (0) ∈ R. For all x ∈ [0, L], we obtain

Q 2 0 2h 2 + g(h + z b ) = Q 2 0 2h 2 (0) + g(h(0) + z b (0)) =: M 0 ,
meaning that the topography z b satisfies

z b = M 0 g - Q 2 0 2gh 2 -h. ( 11
)
Remark 1 Define the Froude number for the steady state by F r = u/ √ gh. The situation F r < 1 corresponds to the subcritical case (i.e. the flow regime is fluvial) while F r > 1 is to the supercritical case (i.e. the flow regime is torrential). In particular, the threshold value of h for F r = 1 is given by

h c := ( Q 2 0 g ) 1 3
.

Because of [START_REF] Han | Photosynthesis-irradiance response at physiological level: A mechanistic model[END_REF], h is the solution of a third order polynomial equation. Given a smooth topography z b , if

h c + z b + Q 2 0 2gh 2 c - M 0 g < 0.
there exists a unique positive smooth solution of (11) which satisfies the subcritical flow condition (see [START_REF] Michel-Dansac | A well-balanced scheme for the shallow-water equations with topography[END_REF]Lemma 1]). The velocity u can then be obtained thanks to (9).

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 ∇ • u = 0 with u = (u(x), w(x, z)). Here, w(x, z) is the vertical velocity. This implies that ∂ x u + ∂ z w = 0. (12) Integrating [START_REF] Lamare | Gradientbased optimization of a rotating algal biofilm process[END_REF] from z b to z gives:

0 = z z b ∂ x u(x) + ∂ z w(x, z) dz = ∂ x z z b u(x)dz + z z b ∂ z w(x, z)dz = ∂ x (z -z b )u(x) + w(x, z) -w(x, z b ) = (z -z b )∂ x u(x) -u(x)∂ x z b + w(x, z),
where we have used the kinematic condition at the bottom (i.e. w(x, z b ) = u(x)∂ x z b ). It then follows from [START_REF] Han | Photosynthesis-irradiance response at physiological level: A mechanistic model[END_REF] that

w(x, z) = ( M 0 g - 3u 2 (x) 2g -z)u (x).
The Lagrangian trajectory is characterized by the system ẋ(t)

ż(t) = u(x(t)) w(x(t), z(t)) . ( 13
)
Remark 2 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.

Up to now and for the sake of simplicity, we denote by z(x) the depth of a particle at x. From (13), we get

z := ż ẋ = ( M 0 g - 3u 2 2g -z) u u . (14)
Note that from ( 9) and ( 11), we have

η = h + z b = M0 g -u 2 2g
, which implies η = -uu /g. Multiplying ( 14) both sides by u and using the formulation of η, η gives

z u + zu = (η - u 2 g )u = ηu + η u,
which implies that u(z -η) = 0. We then obtain

z(x) = η(x) + u(0) u(x) (z(0) -η(0)). ( 15
)
The computation of trajectories in the shallow water system can be carried out with this formula, where u(0), z(0) and η(0) = h(0) + z b (0) are given as boundary conditions, whereas u(x) can be obtained by ( 9) and ( 11). On the contrary, h(x), hence u(x) will be considered in section 3 as a parameter to be optimized, and (11) will be used to (post-)compute z b (x).

Remark 3 Since Q 0 is chosen to be positive, h is necessarily positive and so does u from [START_REF] Grenier | Worldwide perspectives of rotating algal biofilm up-scaling[END_REF].

Moreover, if z(0) belongs to [z b (0), η(0)], then z(x) belongs to [z b (x), η(x)].
In particular, choosing z(0) = z b (0) in (15) and using [START_REF] Grenier | Worldwide perspectives of rotating algal biofilm up-scaling[END_REF] gives z(x) = z b (x). In the same way, we find that z(x) = η(x) when z(0) = η(0).

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:

I(x, z) = I s exp -ε(η(x) -z) . ( 16
)
Here ε is the light extinction coefficient. We assume that the system is perfectly mixed so that the concentration of the biomass X in ( 4) is homogeneous and ε is constant. Combining ( 16) with ( 15), we get the following expression for the captured light intensity along the trajectory z

I(x, z(x)) = I s exp -ε u(0) u(x) (η(0) -z(0)) . (17)
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 ( 14) with (2), we find a time-free reformulation, namely

C := Ċ ẋ = - α(I) u C + β(I) u , (18)
where all the functions on the right-hand side only depend on x.

Optimization problem

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.

Optimization functional

The average net specific growth rate over the domain is defined by

μ := 1 L L 0 1 h(x) η(x) z b (x) µ C(x, z), I(x, z) dzdx. ( 19
)
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:

z i (0) = η(0) - i -1 2 N z h(0), i = 1, . . . , N z .
From [START_REF] Masojídek | A closed solar photobioreactor for cultivation of microalgae under supra-high irradiance: Basic design and performance[END_REF], we obtain

z i (x) -z i+1 (x) = 1 N z h(x), i = 1, . . . , N z ,
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 ( 19) by

μ∆ : = 1 L L 0 1 h(x) 1 N z h(x) Nz i=1 µ(C i (x), I i (x))dx = 1 LN z Nz i=1 L 0 µ(C i (x), I i (x))dx. (20)
From now on, we focus only on the subcritical case, i.e. F r < 1, see Remark 1. 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

V = L 0 h(x)dx. (21)
Therefore, we choose to parameterize h by a vector a := [a 1 , • • • , a N ] ∈ 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 [START_REF] Han | Photosynthesis-irradiance response at physiological level: A mechanistic model[END_REF].

Remark 4

In usual shallow water solver, equations of type [START_REF] Han | Photosynthesis-irradiance response at physiological level: A mechanistic model[END_REF] 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.

Our goal is to optimize the topography to maximize μ∆ . Up to now, we denote explicitly that the functions depend on a. In the case of one trajectory, we consider the functional

J(a) = 1 L L 0 -γ(I(x; a))C(x) + ζ(I(x; a)) u(x; a) dx,
where C satisfies the following parameterized version of ( 18)

C + α(I(x; a)) u(x; a) C = β(I(x; a)) u(x; a) . ( 22
)
The optimization problem then reads:

Find a * solving the maximization problem:

max a∈R N J(a). (23)

Optimality System for maximizing growth rate

In this section, we introduce the optimality system derived from the optimization problem (23) 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.

One trajectory problem

Define the Lagrangian of Problem (23) by

L(C, p, a) = 1 L L 0 -γ(I(a))C + ζ(I(a)) u(a) dx - L 0 p C + α(I(a))C -β(I(a)) u(a) dx,
where p is the Lagrange multiplier associated with the constraint (22).

The optimality system is obtained by cancelling all the partial derivatives of L. Differentiating L with respect to p and equating the resulting expression to zero gives (22). Integrating the terms pC dx on the interval [0, L] by parts and differentiating L with respect to C and C(L) gives rise to (25)

∂ C L = p -p α(I(a)) u(a) -1 L γ(I(a)) u(a) ∂ C(L) L = p(L). (24)

Multiple trajectories problem

We now extend the previous framework to deal with multiple trajectories. From [START_REF] Sung | Design of experiments and sensitivity analysis for microalgal bioreactor systems[END_REF], the objective functional is defined by:

μ∆ (a) = 1 LN z Nz i=1 L 0 -γ(I i (a))C i + ζ(I i (a)) u(a) dx, (26)
where C i satisfies (22) for the light intensity I i with i = 1, . . . , 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 (24), where C is replaced by C i .

Denoting by p i the associated Lagrange multipliers, the partial derivative ∂ a L (hence the gradient ∇μ ∆ (a)) is then given by

∂ a L = 1 LN z Nz i=1 L 0 -γ (I i (a))C i + ζ (I i (a)) u(a) ∂ a I i (a)dx - 1 LN z Nz i=1 L 0 -γ(I i (a))C i + ζ(I i (a)) u 2 (a) ∂ a u(a)dx + Nz i=1 L 0 p i -α (I i (a))C i + β (I i (a)) u(a) ∂ a I i (a)dx - Nz i=1 L 0 p i -α(I i (a))C i + β(I i (a)) u 2 (a)
∂ a u(a)dx.

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 [START_REF] Posten | Microalgae Biotechnology[END_REF]:

Given I ∈ L 2 (0, L; R), find (C 0 , C) ∈ R × C(0, L; R) such that C (x) = -α(I(x))C(x)+β(I(x)) u(x) , x ∈ [0, L] C(L) = C(0) = C 0 . (27)
The The proof is given in Appendix A. The optimality system is obtained as in the previous sections.

Remark 6 Note that the periodicity of C implies that p is also periodic. Indeed, since C(L) = C(0), then by differentiating L with respect to C(L), we have

∂ C(L) L = p(L) -p(0).
so that equating the above equation to zero gives the periodicity for p. Similar computations apply for multiple trajectories problem.

Next, we state a result for Problem (23) in the case where C is periodic.

Proposition 7 Assume h is of the form h(x; a) = a 0 + N n=1 a n sin(2nπ x L ), (29)
with a := [a 1 , . . . , a N ] chosen such that the subcritical condition is satisfied (see Remark 1). Then ∇J(0) = 0.

As a consequence, the flat topography is critical point for Problem (23). Note that the larger N is, the less valid is our hydrodynamic model valid, see Section 2.2, where a smooth topography is assumed. Hence, limit situations where N → +∞ are not considered in what follows.

PROOF. Let q ∈ [0, 1], consider an initial position by z(0; a) = η(0; a) -qh(0; a).

Replacing z(0; a) by the previous expression in [START_REF] Muñoz-Tamayo | Optimizing microalgal production in raceway systems[END_REF] and using (9), gives

I(a) = I s exp(-εq Q 0 u(a) ) = I s exp(-εqh(a)). ( 31
)
Introduce a * = [0, . . . , 0]. Replacing a by a * in ( 9) and (31), we find that u(a * ) and I(a * ) are both constants with respect to x. For a given C(0), (22) gives

C(x) = e -α(I(a * )) u(a * ) x C(0) + β(I(a * )) α(I(a * )) (1 -e -α(I(a * )) u(a * ) x ). ( 32
)
Since C is periodic (i.e. C(L) = C(0)), we get from the previous equation that C(0) = β(I(a * )) α(I(a * )) . Inserting this value in (32), we find

C(x) = β(I(a * )) α(I(a * )) , ∀x ∈ [0, L],
which corresponds to the steady state ( 6), that we denote by C * hereafter. We then compute all the partial derivatives appearing in (25). We get: and the second integral can be computed as

∂ an u(a) = - Q 0 h 2 (a) ∂ an h(a
L 0 -γ(I(a * ))C * + ζ(I(a * )) u 2 (a * ) ∂ an u(a * )dx = -γ(I(a * ))C * + ζ(I(a * )) u 2 (a * ) L 0 ∂ an u(a * )dx =0.
Since p is also periodic (see Remark 6), it can be proved in much the same way that the two last terms in (25) also cancels when a = a * . The result follows.

Remark 8 Numerically, we observe that the flat topography is actually optimal in the periodic case.

Remark 9 Note that if C is not periodic, (32) 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 (25), 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 4.3.2).

Non-constant volume problem for maximizing areal productivity

In general, the extinction coefficient is an affine function with respect to the biomass

X [13] ε(X) = α 0 X + α 1 , (33)
where α 0 > 0 the specific light extinction coefficient of the microalgae specie and α 1 the background turbidity that summarizes the light absorption and diffusion due to all non-microalgae components. Recall that according to (4), the biomass concentration X depends on the growth rate µ and on the harvesting of the raceway D.

Remark 10 A standard criterion to determine a relevant value of X [14,10], and then an harvesting frequency, consists in regulating it to a value such that the steady-state value of the growth rate µ at the (average) bottom depth zb is 0, i.e., µ(I zb ) = 0, (34) where µ(I) is given by [START_REF] Chiaramonti | Review of energy balance in raceway ponds for microalgae cultivation: Re-thinking a traditional system is possible[END_REF]. Solving (34) provides a value of I zb thus of the extinction coefficient ε(X), and finally of the biomass concentration X. This assumption is usually called the compensation condition, for which the growth at the bottom compensates exactly the respiration.

In this framework, maximizing areal productivity is a relevant target. Productivity per unit of surface for a given biomass concentration X is given by:

Π = μX V S , (35)
where S is the ground surface of the raceway pond and μ is defined in [START_REF] René | An outlook on microalgal biofuels[END_REF]. 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 (34) and we assume that it is maintained constant to achieve the compensation conditions, see Remark 10.

Let us keep the parameterization h given in (29). The volume of the system satisfies V = a 0 L (see ( 21)). Moreover, since we consider a 1D framework, we have S = L. The biomass concentration X is chosen to be solution of (34).

We have

X(a 0 ) V (a 0 ) S = α 2 -α 3 a 0 ,
where

α 2 = 1 α0 ln( Is Iz b
) and α 3 = α1 α0 , and α 0 , α 1 are given in (33). The detail of the computation is given in Appendix C. Let us then consider the extended parameter vector ã := [a 0 , a]. For simplicity, we again consider one single trajectory and optimize the objective function associated with (35) defined by

J(ã) = α 2 -α 3 a 0 L L 0 -γ(I(ã))C + ζ(I(ã)) u(ã) dx,
where C satisfies ( 22) with a replacing by ã. The corresponding optimization problem reads:

Find ã * solving the maximization problem:

max ã∈R N +1 J(ã). ( 36
)
The Lagrangian of Problem (36) is

L(C, p, ã) = α 2 -α 3 a 0 L L 0 -γ(I(ã))C + ζ(I(ã)) u(ã) dx - L 0 p C + α(I(ã))C -β(I(ã)) u(ã) dx,
where p is the Lagrangian multiplier associated with the constraint (22) for C. By computations similar to that of the previous section, we find as optimality system

∂ C L = p -p α(I(ã)) u(ã) -γ(I(ã)) u(ã) α2-α3a0 L ∂ C(L) L = p(L).
However, there is an extra element in the gradient J(ã), namely the derivative with respect to a 0 , so that ∇ J(ã

) := [∂ a0 L, ∂ a L],
where

∂ a0 L = α 2 -α 3 a 0 L L 0 -γ (I(ã))C + ζ (I(ã)) u(ã) ∂ a0 I(ã)dx - α 2 -α 3 a 0 L L 0 -γ(I(ã))C + ζ(I(ã)) u 2 (ã) ∂ a0 u(ã)dx - α 3 L L 0 -γ(I(ã))C + ζ(I(ã)) u(ã) dx + L 0 p -α (I(ã))C + β (I(ã)) u(ã) ∂ a0 I(ã)dx - L 0 p -α(I(ã))C + β(I(ã)) u 2 (ã) ∂ a0 u(ã)dx, (37)
∂ a L = α 2 -α 3 a 0 L L 0 -γ (I(ã))C + ζ (I(ã)) u(ã) ∂ a I(ã)dx - α 2 -α 3 a 0 L L 0 -γ(I(ã))C + ζ(I(ã)) u 2 (ã) ∂ a u(ã)dx + L 0 p -α (I(ã))C + β (I(ã)) u(ã) ∂ a I(ã)dx - L 0 p -α(I(ã))C + β(I(ã)) u 2 (ã) ∂ a u(ã)dx. (38)
However, we still find that the flat topography is optimal in that case.

Proposition 11 Assume h is of the form (29) and C is periodic, then

J(ã * ) = 0 with ã * = [a * 0 , 0, • • • , 0].
PROOF. The proof is similar to that of Proposition 7. For a given z(0; ã) defined in (30), using the result in Appendix C to replace ε in (31), we have

I(ã) = I s exp - 1 a 0 ln( I s I zb )qh(ã) .
The two extra partial derivatives appearing in (37) are computed by and ∂ a0 I(ã * ) = 0. Hence the first integral and the fourth integral in ∂ a0 L in (37) are zero. Following a reasoning similar to that of Proposition 7, we find C * is constant. Hence, the last integral in ∂ a0 L is zero. Finally for

∂ a0 u(ã) = - Q 0 h 2 (ã) ∂ a0 h(ã), ∂ a0 I(ã) = -(- 1 a 2 0 h(ã) + 1 a 0 ∂ a0 h(ã)) ln( I s I zb )qI(ã),
a * 0 = α 2 2α 3 , (39)
the sum of the second and the third integral in ∂ a0 L is zero.

As for (38), it follows from computations similar to that of the proof of Proposition 7 that this term equals zero, which completes the proof.

Numerical Experiments

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

Numerical Methods

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

Gradient Algorithm

We detail now a gradient-based optimization procedure and apply it to Problem (23). The procedure is given in Algorithm 1.

Algorithm 1 Gradient-based optimization algorithm 1: Input: Tol> 0, ρ > 0.

2: Initial guess: a. Compute u by ( 9) and I by [START_REF] Muñoz-Tamayo | Optimizing microalgal production in raceway systems[END_REF].

7:

Set C as the solution of (22).

8:

Set p as the solution of (24).

9:

Compute the gradient ∇J by (25). Set err := ∇J . 12: end while Note that h ∞ := max x∈[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 1 (and in the range of industrial constraints, see [START_REF] Chiaramonti | Review of energy balance in raceway ponds for microalgae cultivation: Re-thinking a traditional system is possible[END_REF]).

A similar algorithm can be considered to tackle the multiple trajectories Problem (26). 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 3.2.4 can still be solved with Algorithm 1 by adapting the objective functional J and the gradient ∇ J.

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 ∆x, set N x = [L/∆x] and x nx = n x ∆x for n x = 0, . . . , N x . We choose to use the Heun's method for computing C via (22). 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.

Parameter settings

We now detail the parameters used in our simulations.

Parameterization of h

As mentioned in Proposition 7, we parameterize h by a truncated Fourier series (29). The parameters to be optimized will be the Fourier coefficients a := [a 1 , . . . , a N ]. Note that we choose to fix a 0 , since it is related to the volume V of the raceway with V = a 0 L.

Parameter for the models

The spatial increment is set to ∆x = 0.01 m so that the convergence of the numerical scheme has been ensured, and we set the raceway length L = 100 m, the averaged discharge Q 0 = 0.04 m 2 s -1 , the initial value a 0 (= h(0; a)) = 0.4 m and z b (0) = -0.4 m to stay in standard ranges for a raceway. The free-fall acceleration g = 9.81 m s -2 . All the numerical parameters values for Han's model are taken from [START_REF] Grenier | Worldwide perspectives of rotating algal biofilm up-scaling[END_REF] and given in Table 1. In order to determinate the light extinction coefficient ε, 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 = 2000 µmol m -2 s -1 which approximates the maximum light intensity, e.g., at summer in the south of France. Then ε can be computed by ε = (1/h(0; a)) ln(I s /I b ).

Note that in the case where a 0 is also a parameter to be optimized, we take from [START_REF] Martínez | Theory of turbid microalgae cultures[END_REF] the specific light extinction coefficient of the microalgae specie α 0 = 0.2 m 2 • gC the background turbidity and α 1 = 10 m -1 , under these settings, a * 0 defined in (39) is around 0.316 m.

Numerical results

We now test the influence of various parameters on optimal topographies. In all our experiments, we always observe that the topographies satisfy h ∞ > h c .

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 μ∆ for each N z . The results are shown in Fig. 3. 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.

Influence of the initial condition

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. 4. A slight difference between the two optimal topographies is observed. Note that this difference remains when the spatial increment ∆x goes to zero.

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, [START_REF] Chiaramonti | Review of energy balance in raceway ponds for microalgae cultivation: Re-thinking a traditional system is possible[END_REF][START_REF] Grognard | Optimal strategies for biomass productivity maximization in a photobioreactor using natural light[END_REF][START_REF] Masojídek | A closed solar photobioreactor for cultivation of microalgae under supra-high irradiance: Basic design and performance[END_REF][START_REF] Sung | Design of experiments and sensitivity analysis for microalgal bioreactor systems[END_REF] and keep all the other parameter values. Table 2 shows the optimal value of the objective function μ∆ (a * ) for different order of truncation N of the Fourier series. As expected, the result shows a slight increase of the optimal value of the objective function μ∆ (a * ) when N becomes larger. However, the corresponding values of μ∆ (a * ) remain close to the one associated with a flat topography.

PROOF. Let x ∈ [0, L], as a weak solution of (A.1), the function Z is almost everywhere differentiable and a direct calculation gives d Z(x) x , we obtain: de

2 2 dx = 2 < Z(x), - α(I(x)) u(x) Z(x) > ≤ -2 k r u ∞ Z(x) 2
2 kr u ∞ x Z(x) 2 2 dx = e 2 kr u ∞ x f (x),
which gives by integration over [0, x] Z(x)

2 2 = e -2 kr u ∞ x Z 0 2 2 + x 0 e 2 kr u ∞ (s-x) f (s)ds.
Because of (A.3), f (x) ≤ 0, this concludes the proof.

Let us now give the proof for Theorem 5.

PROOF. Introduce the mapping Φ : R → R defined by

Φ(C 0 ) := C(L),
where C is the weak solution of [START_REF] Posten | Microalgae Biotechnology[END_REF]. Given C1 0 ∈ R and C2 0 ∈ R, define Z 0 = C2 0 -C1 0 and Z(x) = C2 (x) -C1 (x), where C1 and C2 are the weak solutions of (18) with initial conditions C1 0 and C2 0 , respectively. Subtracting the corresponding weak representations, we obtain that Z satisfies the assumptions of Lemma 12 so that (A.2) holds. As a consequence, C2 (L) -C1 (L) 2 ≤ e -kr u ∞ L C2 0 -C1 0 2 , which implies that Φ is a contraction. Applying Banach fixed-point theorem, it follows that there exists a unique C 0 ∈ R such that Φ(C 0 ) = C 0 . The corresponding weak solution C of ( 18) satisfies (28).

B System with a paddle-wheel

Let us denote by P the permutation matrix associated with π i.e., 1 as entries on the anti-diagonal. Let us denote by C 1 (resp. C 2 ) the photo-inhibition state for the first (resp. second) lap of the raceway. We then assume that the state C is 2-periodic, meaning that C 1 (0) = P C 2 (L). From (26), we define the objective function by

C Computation for X

Let us consider a 1D system with vertical light attenuation. We consider the simple configuration of a static growth rate where Han dynamics is considered at equilibrium. The optimal biomass shades light such that, at the pond bottom, we have (34). Replacing α, β, γ, ζ by their own definitions in ( 5) and ( 6), we find a second order equation for I k d τ R(σI) 2 + (k r τ σR -k r kσ)I + k r R = 0. Denote by I zb the suitable root of this equation (where 0 < I zb < I s ). From [START_REF] Masojídek | A closed solar photobioreactor for cultivation of microalgae under supra-high irradiance: Basic design and performance[END_REF] and [START_REF] Muñoz-Tamayo | Optimizing microalgal production in raceway systems[END_REF], the light intensity at the (averaged) bottom is given by I zb = I s exp(-ε(X)a 0 ). This implies ε(X) = 1 a 0 ln( I s I zb ).

Recall that ε(X) = α 0 X + α 1 , hence we have X(a 0 ) =

Fig. 1 .

 1 Fig. 1. Scheme of the Han model, representing the probability of state transition, as a function of the photon flux density.

Fig. 2 .

 2 Fig. 2. Representation of the hydrodynamic model.

0 p

 0 Given a vector a, let us still denote by C, p the corresponding solutions of (22) and (24). The gradient ∇J(a) is obtained by ∇J(a) = ∂ a L, where ∂ a L = 1 L L 0 -γ (I(a))C + ζ (I(a)) u(a) ∂ a I(a)dx -I(a))C + ζ(I(a)) u 2 (a) ∂ a u(a)dx + L -α (I(a))C + β (I(a)) u(a) ∂ a I(a)dx -I(a))C + β(I(a)) u 2 (a) ∂ a u(a)dx.

Theorem 5

 5 following theorem gives the existence and uniqueness of a (weak) solution of the problem (27) Given I ∈ L 2 (0, L; R), there exists a unique couple (C 0 , C) ∈ R × C(0, L; R) satisfying C(x) = C 0 + x 0 -α(I(s))C(s)+β(I(s)) u(s) ds , C(L) = C 0 (28) for all x ∈ [0, L].

0 ∂

 0 ), ∂ an I(a) = -εqI(a)∂ an h(a), but ∂ an h(a) = sin(2nπ x L ), so that the integral from 0 to L of ∂ an h(a) is zero. Replacing a by a * in (25) , the first integral can be computed as L 0 -γ (I(a * ))C * + ζ (I(a * )) u(a * ) ∂ an I(a * )dx = -γ (I(a * ))C * + ζ (I(a * )) u(a * ) L an I(a * )dx =0.

where ∂ a0 h(ã) = 1 .

 1 The flat topography corresponds to the parameter ã * = [a * 0 , 0, . . . , 0]. Then ∂ a0 u(ã * ) = -

3 :: a 4 :

 34 OutputSet err := Tol+1 and define h by (29) using the input data. 5: while err >Tol and h ∞ > h c do 6:

1 )Fig. 3 .Fig. 4 .

 134 Fig. 3. The value of the functional μ∆ for Nz = [1, 100].

2 ,(A. 3 )

 23 where we have used the fact that α(I) = β(I) + k r , β(I) ≥ 0 and 0 < u < u ∞ . Defining f (x) := d Z(x)

  I i (a))C j i + ζ(I i (a)) u(a)dx.Let us still denote by L the Lagrangian associated with this optimization problem. It can be writtenL(C, p,a) I i (a))C j i + ζ(I i (a)) i (a))C j i -β(I i (a)) u(a) dx,where p j i is the Lagrangian multiplier associated with the constraint (22) for C j i . Integrating the terms p j i C j i dx on the interval [0, L] by parts, we get i dx +p j i (L)C j i (L) -p j i (0)C j i (0).By definition of our paddle wheel model, C 2 (0) = P C 1 (L). Differentiating L with respect to C 1 (L) and C 2 (L), we find∂ C 1 (L) L = p 1 (L) -p 2 (0)P, ∂ C 2 (L) L = p 2 (L) -p 1 (0)P.Then differentiating L with respect to C j i gives the evolution of p j i by∂ C j i (L) L = p j i -p j i α(I i (a)) u(a) -1 2LN z γ(I i (a)) u(a) .Finally, the partial derivative ∂ a L (hence the gradient) is given by I i (a))C j i + ζ (I i (a)) u(a) ∂ a I i (a)dx -I i (a))C j i + ζ(I i (a)) u 2 (a)∂ a u(a)dx+ I i (a))C j i + β (I i (a)) u(a) ∂ a I i (a)dx -I i (a))C j i + β(I i (a)) u 2 (a) ∂ a u(a)dx.

Table 1

 1

	Parameter values for Han Model.
	kr	6.8 10 -3	s -1
	k d	2.99 10 -4	-
	τ	0.25	s
	σ	0.047	m 2 • (µmol) -1
	k	8.7 10 -6	-
	R	1.389 10 -7	s -1

Table 2

 2 The objective function for different order of truncation of the Fourier series.

	N	Iter μ∆(d -1 ) log 10 ( ∇J )
	0	0	0.94764	-
	5	16	0.96942	-10.368944
	10	17	0.97509	-10.262891
	15	17	0.97734	-10.010322
	20	18	0.97851	-10.196963

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 [START_REF] Bernard | Optimizing microalgal productivity in raceway ponds through a controlled mixing device[END_REF], where other choices will be investigated.

Meanwhile, this permutation matrix P corresponds to the permutation

where we use the standard notation of cycles in the symmetric group. Note that π is of order two. The photoinhibition state C is then set to be 2-periodic (i.e C 1 (0) = P C 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. 5 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 [START_REF] Bernard | Optimizing microalgal productivity in raceway ponds through a controlled mixing device[END_REF] in the case of a flat topography. The resulting optimal topography in this case is not flat. However, the increase in the optimal value of the objective function μ∆ 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.

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.

A Proof for Theorem 5

Let us give the details of the proof for Theorem 5. To prove the existence and uniqueness of (weak) a solution of this problem, we need the following result.