Bernard Tamba
email: bernard-tamba.sandouno@inria.fr

Sandouno Sogudo

Yamen Alsaba
email: yamen.alsaba@zoneadsl.com

Chadi Barakat
email: chadi.barakat@inria.fr

Walid Dabbous
email: walid.dabbous@inria.fr

Thierry Turletti
email: thierry.turletti@inria.fr

A Novel Approach to Mobile Outdoor QoS Map Generation

Keywords: Ray Tracing, QoS map, optimization

Ray Tracing is an electromagnetic wave propagation modeling approach used for accurate generation of Quality of Service (QoS) maps in mobile networks. Due to its complexity, current implementation of Ray Tracing fails to generate such maps in wide areas. In this paper, we propose an optimization to Ray Tracing able to accurately generate QoS maps in a reasonable time. Using a site-specific ray launching technique and an alternative to the reception test process, we divide by almost 1200 the execution time of Ray Tracing with less than 2% of memory usage as compared to baseline solutions.

I. INTRODUCTION

Ray Tracing (RT) is an electromagnetic wave propagation modeling approach that accurately estimates the signal power received at a given location. RT's accuracy comes from its ability to consider fine-grained information about the environment of interest. Indeed, according to RT, the interaction between waves and buildings is taken into account by different multipath propagation mechanisms [START_REF] Yun | Ray Tracing for Radio Propagation Modeling: Principles and Applications[END_REF], [START_REF] Erunkulu | Cellular Communications Coverage Prediction Techniques: A Survey and Comparison[END_REF]. With these mechanisms, the path followed by waves from a transmitter to a receiver is first traced and the corresponding power and bitrate are then computed [START_REF] Liu | Modelling radio wave propagation in tunnels with ray-tracing method[END_REF]. Due to the high level of accuracy obtained with RT, it is considered to be an appropriate tool for Quality-of-Service (QoS) maps generation in mobile wireless networks. Examples of QoS maps include the download bitrate and the signal power. Nevertheless, due to hardware limitations, to the big size of the area of interest and to the huge number of antennas and buildings present in the area, it is often impossible to generate those maps using RT. This is mainly due to the discretization process that precedes the map generation. Indeed, during this process, the area is represented as a large set of receiving points, where the bitrate is calculated to produce the map. To obtain a precise map, the distance between those points must be as small as possible. This may lead to thousands and even millions of points generated. The obtained points are considered as receivers and a reception test This work is supported by French ANRT agency under contract CIFRE 2020/1392. is performed between them and the rays to determine which rays are being received by the receivers. This process leads to high computational load, high execution time and often memory limit exceeding errors, hence causing the inability to produce the maps. For instance, after 16 hours of execution time on an Intel Core i7 machine with 16GB memory, Matlab's RT implementation crashed due to a memory limit exceeding error and was unable to generate the map in an 1072m x 730m area. This unfortunately makes it impossible for raw RT implementation to generate accurate maps at large scale.

On the other hand, the different RT acceleration techniques such as the space division techniques [START_REF] Zhengqing Yun | A fast ray tracing procedure using space division with uniform rectangular grid[END_REF], although very efficient for estimating the signal power of a single receiver, may fail when it comes to generating QoS maps since none of them directly tackles the overhead due to the huge number of receiving points in the area. These techniques mostly care about reducing the complexity related to the ray-object intersection test [START_REF] Zhengqing Yun | A fast ray tracing procedure using space division with uniform rectangular grid[END_REF]. Therefore, as the size of the environment and the number of receiving points get bigger, existing techniques become inefficient due to the big matrix operations performed during the reception test process which highly increases the execution time and consumes a lot of memory as in the case of Matlab.

To overcome the inability of the current models and techniques, we introduce in this paper an optimized RT solution able to accurately generate QoS maps in a reasonable time. Where current models as the one of Matlab take 16 hours to execute and often fail due to memory limitations, our solution produces the QoS map within 50 seconds. This is made possible on one hand by the use of a site-specific ray generation technique that we introduced in [START_REF] Tamba Sandouno | Site-specific Ray generation for Accurate Estimation of Signal Power[END_REF]. This technique helps to launch the optimal number of rays in order to fully cover the area. On the other hand, and instead of discretizing the area, our contribution in this paper is to merely capture the footprint of rays on a given plane rather than performing reception tests on a set of points. By doing so, we remove the discretization and the reception test processes which account for most of the overhead in RT when generating QoS maps. The main contribution of this paper explained above helps to generate accurate and precise QoS maps at a city or country level in a reasonable time without any memory constraint. Furthermore, by removing the discretization process, our method also removes the need for performing any further interpolation in the discrete case. To the best of our knowledge, this is the fastest tool able to generate QoS maps using RT. As we will see later, we can achieve this goal without compromising the accuracy of the obtained QoS maps.

The rest of the paper is organized as follows. In Section II, some RT fundamentals are explained. The related work is presented in Section III. The system overview is shown in Section IV. We present the validation and the performance evaluation of our solution in Section V. Finally, conclusions and perspectives of our work are discussed in Section VI.

II. RAY TRACING FUNDAMENTALS

Ray Tracing (RT) consists of simulating the different paths followed by rays from a transmitter to a receiver. During this simulation, multiple rays are launched from the transmitter and each of them is tracked by considering propagation mechanisms such as reflection and diffraction. To determine if a ray is received by a receiver, a reception test is performed. Rays are modeled as cones and the reception test consists in checking if the receiver is inside the cones' cross-sections. The radius R i of a ray cone's cross-section is given by (1) with α i being the maximum angular separation between adjacent rays and d i the distance travelled by the i th ray.

R i = α i d i √ 3 (1)
The i th ray is received if it passes through the reception sphere of radius R i centered at the receiver [START_REF] Zhengqing Yun | Development of a new shooting-and-bouncing ray (SBR) tracing method that avoids ray double counting[END_REF]. From Fig. 1, the reception sphere is centered at point P (the receiver). Ray i is received by point P if the distance between P and its orthogonal projection on the ray is less than or equal to the radius of the ray, i.e if ∥ --→ HP ∥ ≤ R i . Point H, the projection of point P is found using (2):

AH = -→ AP • ⃗ u H = AH * ⃗ u + A (2)
Now, when it comes to generating QoS maps, the area is discretized into multiple receivers. With M receivers, N rays and K reflections, the matrix calculation in Formula (2) must be conducted N × M × K times. For each reflection, the condition ∥ --→ HP ∥ ≤ R i has to be checked to identify for each receiver the associated rays that meet the reception test condition. The corresponding received powers are then summed to compute the bitrate for each receiver and to produce the QoS map. The bitrate value is given by the Shannon capacity formula in [START_REF] Erunkulu | Cellular Communications Coverage Prediction Techniques: A Survey and Comparison[END_REF], where B is the channel frequency bandwidth. The signal source (S) is the antenna with the highest signal power whereas other antennas at the same frequency act as sources of interference (I). N denotes the total noise power at the receiver [START_REF]Matlab, SINR[END_REF].

Bitrate = B * log 2 (1 + S I + N ) (3)
III. RELATED WORK Different techniques proposed in the literature are used to generate QoS maps. These techniques differ from each other by their level of accuracy and their complexity. The empirical propagation modeling approach for instance, allows QoS maps to be generated after a small execution time using offline calibrated models on real data. These models are fast because they do not consider thorough details about the environment of interest such as buildings and roads. Only with parameters such as the frequency and the distance, these models are able to estimate the received signal power [START_REF] Sharma | Comparative Analysis of Propagation Path loss Models with Field Measured Data[END_REF]. Although very fast for generating QoS maps, these models suffer from a low accuracy and a weak capacity to account for the details of the propagation environment.

In the same manner, another category of models called the stochastic models is even faster than the empirical ones when generating QoS maps. With stochastic models, the environment is modeled as a set of random variables. A probability density function is used afterwards to estimate the signal power. These models can be very useful for generating QoS maps in wide areas since they do not need a lot of information about the propagation environment for generating the maps. Thus, the maps can be generated with a small execution time [START_REF] Sharma | Comparative Analysis of Propagation Path loss Models with Field Measured Data[END_REF]. Nevertheless, due to the randomness of this technique and its neglect of environment details, it is the least accurate among all the propagation modeling approaches.

On the other hand, RT can produce accurate QoS maps due to its ability to consider the full 3D map details of the propagation environment. Raw implementation of RT as the one of Matlab allows for example to generate the Signal to Interference plus Noise Ratio (SINR) map of a given area. The SINR can then be used to produce maps of other performance metrics such as the bitrate. Different parameters such as the number of rays to be launched and the resolution (maximum distance between receiving points) are available and can be tuned. For a small area and a large distance between receivers, this model performs well by showing a heat map giving the SINR level of each location in the area of interest [START_REF]Matlab, SINR[END_REF]. However, when it comes to generating precise maps (small distance between receiving points) in wide areas, the execution time of Matlab's RT becomes very high. Moreover, due to memory limitations, it crashes and gets unable to produce the SINR map. To avoid such issues, Matlab limits for example the coverage area around an antenna to 500 meters. Although this may be helpful in some cases to reduce the complexity of the model, severe accuracy problems occur. Indeed, in real life, some antennas can cover more than 10 kilometers around them and within this radius, hundreds of other antennas may be present. Furthermore, with this setting, potential sources of interference are ignored by Matlab, leading to the generation of QoS maps showing unrealistic results. The larger the number of antennas in the environment of propagation, the higher is the probability of having sources of interference and the worse will be the results from Matlab's RT implementation. This shows the inability of the state-of-the-art raw RT implementation to be used in real life to generate accurate QoS maps at large scale.

To the best of our knowledge, none of the tools in the literature tackles the trade-off between the accuracy and the computational complexity of RT in order to generate accurate QoS maps in a reasonable time. In order to overcome this issue, we introduce in this paper an optimization to RT that allows to accurately generate QoS maps with small execution time. To achieve this trade-off, we simply perform the intersection of ray cones with a fixed horizontal plane Z instead of performing the reception test explained in Section II on different receiving points of the map with different heights. We afterwards use a simple condition to perform the ray-plane intersection test. This helps to drastically reduce the execution time since the need of performing all the matrix operations as in the original RT is removed. Moreover, with our approach, it is no longer necessary to perform the discretization of the map at the reception. Rather, the receivers are considered to be all located on a same fixed horizontal plane. Thus, our method is always precise because all the receiving points in the area of interest are considered. We further improve our approach by coupling it with a site-specific ray generation technique we proposed in an earlier work [START_REF] Tamba Sandouno | Site-specific Ray generation for Accurate Estimation of Signal Power[END_REF]. With the latter, we are able to launch the minimum number of rays needed to obtain an accurate map. With less rays launched and the reception test process removed, we will demonstrate next how our new method is able to produce accurate maps in a reasonable time.

IV. SYSTEM OVERVIEW

According to standard Ray Tracing, all the rays launched by an antenna need to be tracked and their intersection test with buildings must be performed upon each reflection. Additionally, the reception test between all the rays and all the considered receiving points is performed until the maximum number of reflections is reached. The larger the number of rays and the number of receiving points, the higher will be the execution time of RT, and the better will be the accuracy and the precision of the produced maps.

To reduce the cost in terms of execution time, we solve the complexity related to two main processes of RT. On one hand, instead of launching rays in all possible directions, we use a site-specific ray generation technique. This technique consists of generating the minimum possible number of rays that fully cover the area of interest without any gaps. With this technique, bunch of rays are saved from being wasted, because only potential rays that can reach the receiver are launched.

With RT, the parameter α used in (1) determines the radius of the ray cone. The smaller it is, the more rays are launched and the more accurate the signal power estimate is, but at the cost of higher execution time. Using a larger separation angle α allows a smaller number of rays to be launched, which reduces the computation time needed during the intersection test with buildings and the reception test with receivers. The challenge is to find the larger value of α that does not compromise the accuracy of the maps. This is what we solved in [START_REF] Tamba Sandouno | Site-specific Ray generation for Accurate Estimation of Signal Power[END_REF], using a site-specific and iterative new ray generation approach. With this approach, we proved that the drop of performance between an angular separation α = 15°and α = 1°is small, i.e., one can reach almost the same level of accuracy with a large value of α than with a small one [START_REF] Tamba Sandouno | Site-specific Ray generation for Accurate Estimation of Signal Power[END_REF].

Fig. 2. Ray-plane intersection

The reception test is another cumbersome process to RT as explained in Section II. Indeed, even when smaller number of rays are launched, the large number of receiving points will still lead to an overhead with high memory consumption and complex matrix manipulation. To avoid this, we consider the receivers to be all located on a given horizontal plane Z plane . With this, all the receivers have the coordinates (X i , Y i , Z plane ). For example, we can consider the receivers to be at the same height Z plane = 1.5 meters. This value typically corresponds to the average height of a person holding a mobile phone.

Fixing the height of the plane has generally little impact as the small variations of the height in the terrain have little effect on the power computed. In future works, we are planning to generalize our method in order to take into account all single terrain variations in the given area. Moreover, fixing Fig. 3. LOS scenario without any building at z plane = 1.5 meters the plane avoids performing the reception test between all the rays and all the receiving points. Rather, an intersection test is performed between ray cones' cross-sections and the horizontal plane of reception to calculate the contribution of the different rays to the QoS map. Since rays are cones, this typically consists in the intersection between a cone and a plane. The plane in the area being horizontal, this intersection is either a circle or an ellipse depending on the incident angle of the ray. For the sake of simplicity, we consider the intersection with the plane to be circles since rays' cross sections are circular themselves with a well-known radius given by (1). The consideration of ellipses would add considerable computation overhead to our approach. However, as we will see in the validation of Section V, the loss in terms of accuracy due to this approximation remains small.

Finally, as shown in Fig. 2, a ray intersects with the plane if condition (4) is met, with Z i = Z A being the z-value of the ray previous reflection point A, c i the z component of the ray's unit direction vector and d i the distance travelled by the ray between its previous reflection point A and its current endpoint.

d i ≥ Z plane -Z i c i (4)
Fig. 3 is an example of intersection with an horizontal plane showing the continuous footprints of rays. By removing the discretization, our approach always produces a precise QoS map without the need of having the receiving points to be as close as possible as in the traditional approach. Furthermore, there is no need to interpolate the received signal strength between the receiving points to find the bitrate for unknown points. Indeed, with rays' footprints, every single receiving point in the area is considered during the generation of the QoS map, which leads to a constantly high precision regardless of the targeted resolution. Note that with our approach, rays may be superposed leading to some areas covered by more than one ray. In those cases, the received power is the sum of the powers of all the rays covering that area, except for double counted rays [8] [4]. In the latter case, the power of the overlapping area corresponds to the power of a single ray. Additionally, we check whether overlapping cones are coming from the same antenna and have the same carrier frequencies in order to consider them or not for the calculation of the interference caused by other antennas. After all these considerations, the bitrate is computed and the QoS map showing the download bitrate as a heat map of the area of interest is generated.

V. PERFORMANCE EVALUATION

When buildings are present in the propagation area, it is not possible for Matlab's RT implementation to have receiving points on the same plane. Indeed, after the discretization process, the heights of the receiving points depend on whether they fall on a building or not. If they fall on a building, their final height is their own height plus the one of the building. Due to receivers having different heights, we could not validate our model using Matlab's RT implementation. To validate our solution, we implemented a raw RT that allows receivers to have the same height regardless of the buildings. However, for the computation load evaluation, the execution time of our solution was compared to the one of Matlab. To do this, we performed different simulations on different terrains with different scenarios. Since our model is free from any discretization unlike the traditional RT, we made possible, for validation purposes, to use our model in such a way it provides the SINR for the individual receivers used in our discrete reference RT implementation.

A. Bitrate validation

In Fig. 4, we see that in the Line-of-Sight (LOS) case, our model has almost the same distribution of the bitrate as our ground truth model with a Root Mean Square Error (RMSE) of 0.01 Mbps. This negligible RMSE means that our model has exactly the same accuracy as the traditional RT in the LOS case. This is because in LOS, the received power is only a function of d 2 [1] (d is the distance of the ray) and since the distances in LOS are typically small, the difference of distances between the traditional RT and our model is small, hence having almost the same received power and consequently the same bitrate. On the other hand, Fig. 5 shows the comparison of the two distributions in case of reflections, with the maximum number of possible reflections sets to 2. We can see that the two models have close distributions with a slight difference in the RMSE. This is because, in Non-LOS, the received power is a function of exp(d) [START_REF] Tamba Sandouno | Site-specific Ray generation for Accurate Estimation of Signal Power[END_REF]. Since in NLOS the distances are typically larger, the difference of distances between the traditional RT and our model is not negligible, leading to this slight increase of the RMSE. This difference in the distributions means that upon generating QoS maps for a given terrain, our model makes an average error of about 2 Mbps for 2 reflections. However, we believe this is still an acceptable level of trade-off in terms of time savings as the following subsection will show.

This validation was performed on 3 terrains in the city center of Nice in France. For the sake of our approach generalization, more terrains with different variations in heights will be considered. For now, the simulations performed in Nice were done in LOS, 2 and 4 reflections cases. From Table I, we can see that our model is robust to the change of terrain, since the difference of RMSE in the 3 cases is negligible, i.e., our model keeps almost the same level accuracy regardless of the terrain. Moreover, we also see that our model reaches an acceptable level of accuracy only with 2 reflections.

Next, we show the effect of the angular separation α on the overall accuracy of our model. From Table II, we can see the variation of the RMSE with respect to the angular separation α for 3 different reflections. We can notice from this table that smaller values of α have smaller RMSEs. Nevertheless, the difference of RMSE between the highest value of α and its lowest value is 0.18 Mbps in LOS, 0.67 Mbps with 2 reflections, and 1.62 Mbps in the last case. This means that the drop of performance when going from a small to a large value of α is small. This implies that it is possible to generate a QoS map with the highest value of α and get almost the same level of accuracy as the smallest ones, with the benefit of a lower execution time.

B. Execution time

For the purpose of execution time evaluation, we compared the running time of our model with Matlab's solution, which is a pure RT implementation. Matlab has two fundamental parameters while generating the QoS map: the maximum angular separation (α) and the resolution. The former has to be chosen between Low, Medium and High. Low is the most accurate angular separation. It launches more rays compared to High where less rays are launched but at the expense of a lower accuracy. On the other hand, the resolution is the maximum distance between receiving points. The smaller it is, the higher is the number of receiving points generated and the more precise the QoS map is.

We performed our simulations on an HP ELITEBOOK 850 G7 laptop with 16GB of memory, and a Core i7 CPU @ On the other hand, since we proved that increasing the value of α within some range does not alter the model's accuracy, we simulated our model only when α = 15°. The latter is the maximum value that α can take due to the small angle approximation made in the site-specific ray generation technique we used [START_REF] Tamba Sandouno | Site-specific Ray generation for Accurate Estimation of Signal Power[END_REF].

Fig. 6 gives the execution time of Matlab as a function of the resolution in the cases where there are 2 and 4 reflections with Low angular separation. In the case of 4 reflections, one must wait almost 1.5 hours to obtain a less precise map and up to 15 hours to obtain a more precise map. In the case of 2 reflections, this value ranges from 45 minutes to 7.5 hours. Given the size of the terrain, this is a very high value. The trend of the graphs shows that this value will increase more and more, that is, the larger the terrain is, the longer the waiting time will be. On the other hand, our approach took 50 seconds and 2 minutes respectively for 2 and 4 reflections. This shows the ability of our model to drastically reduce the execution time of RT and to easily scale while producing the map in a reasonable time.

Table III gives a summary of the average execution time of Matlab in seconds as a function of the resolution compared to the time taken by our model. The values for Matlab are the averages over the results obtained for the 3 angular separations allowed in Matlab: Low, Medium and High. This table shows that where our model takes 5 seconds in LOS to generate the map, Matlab's model takes on average between 102-709 seconds. On the other hand, our model takes respectively 50 and 125 seconds with 2 and 4 reflections while the time of Matlab is way higher. For larger environments, at country level for instance, Matlab may take days or weeks to produce the QoS map making it not effective for such cases. Moreover, when taking a resolution of 2 meters, Matlab failed at producing the map in the case of reflections. This is due to the large memory consumption generated by RT processes. For example, with 2 reflections and after 16 hours of execution time, Matlab crashed due to this memory issue. This shows that Matlab struggles to generate a precise map for larger environments where the number of receiving points is higher. Whereas our model is not very affected by memory limitations. Hence it can be applied in large scale scenarios with a lower execution time as compared to raw RT implementation as the one of Matlab.

VI. CONCLUSION

In this paper, we have shown the inability of the current Ray Tracing based models to accurately generate QoS maps in a reasonable time. Indeed, these models fail at generating QoS maps at large scale. With the use of some optimization techniques, we demonstrated that our model could find a good trade-off between computational complexity and accuracy. It is thus able to generate accurate QoS maps within a reasonable time. The consideration of different terrain elevations in our approach will be the topic of a future research work.

Fig. 1 .

 1 Fig. 1. Reception test process

Fig. 4 .

 4 Fig. 4. CDFs of our approach and a traditional RT: LOS

Fig. 5 .Fig. 6 .

 56 Fig. 5. CDFs of our approach and a traditional RT: 2 reflections

TABLE III MEAN

 III EXECUTION TIME IN SECONDS OF MATLAB VS OUR SOLUTION

	Resolution (meters) 3	4	5	8	10	Our model
	LOS	709	422	282	129	102	5
	2 reflections	12646	7310	4985 1884	1280 50
	4 reflections	23198	13321 8545 3465	2234 125


domain: sde math

text:
Francesca Casagli

Fabrice Beline

Elena Ficara

Olivier Bernard
email: olivier.bernard@inria.fr

Optimizing resource recovery from wastewater with algaebacteria membrane reactors

Keywords: Wastewater remediation, membranes, process optimization, algaebacteria, nitrous oxide, emerging pollutants

Exploiting the combination of algae and bacteria in High Rate Algal/Bacterial Ponds (HRABP) is an emerging approach for wastewater remediation and resource recovery.

In this study, the advantage of adding a solid/liquid separation system to uncouple Hydraulic Retention Time (HRT) and Solid Retention Time (SRT) is explored and quantified. A long-term validated model for HRABP was run to simulate and optimize a system at large scale treating digestate. It is shown that by uncoupling HRT and SRT, adapting the liquid depth and the alkalinity content, the algae productivity increases from 9.0-14.5 g m -2 d -1 (for HRT=SRT in the range of 5 to 10 days) to 20.3 g m -2 d -1 (for HRT =0.2 d and SRT= 2 d). Simulations pointed out that maximizing the algal productivity or the fraction of recovered nitrogen in the algal biomass are conflicting goals that are achieved under different operating conditions. Conditions maximising the algal productivity favour algae and heterotrophic bacteria while algae and nitrifying bacteria dominate the system under those conditions optimizing the efficiency of nitrogen recycling. Finally, increasing the influent alkalinity and adapting the water depth can boost the algal productivity without meeting conditions favourable to N2O emission, opening new perspectives for resource recovery through algal biomass valorisation.

Introduction

Exploiting the combination of algae and bacteria in High Rate Algal/Bacterial Ponds (HRABP) is an emerging technology for wastewater treatment addressing the most critical aspects of wastewater treatment [START_REF] Gonçalves | A review on the use of microalgal consortia for wastewater treatment[END_REF], [START_REF] Wang | Advances in algal-prokaryotic wastewater treatment: A review of nitrogen transformations, reactor configurations and molecular tools[END_REF]. Compared to conventional biological processes, HRABP do not need an external oxygen supply, which is known to account for more than 50% of the energy required for wastewater treatment [START_REF] Lee | Evaluation of aeration energy saving in two modified activated sludge processes[END_REF]. Recycling nitrogen and phosphorus into the algal biomass is another remarkable opportunity offered by this process [START_REF] Di Caprio | Microalgae cultivation by uncoupled nutrient supply in sequencing batch reactor (SBR) integrated with olive mill wastewater treatment[END_REF], since algae can be used as a source of energy or feedstocks [START_REF] Lage | Algal Biomass from Wastewater and Flue Gases as a Source of Bioenergy[END_REF], [START_REF] Shahid | Current advances in treatment technologies for removal of emerging contaminants from water -A critical review[END_REF], i.e. with their lipid fraction for biofuels or proteins for bioplastics [START_REF] Sutherland | 15 years of research on wastewater treatment high rate algal ponds in New Zealand: discoveries and future directions[END_REF], [START_REF] Sutherland | Microalgal bioremediation of emerging contaminants -Opportunities and challenges[END_REF]. Moreover, HRABP have been highlighted for their ability to remove emerging contaminants. These molecules are persistent in the environment, impact on the reproductive systems of aquatic organisms, and bioaccumulate in the food chain [START_REF] Yang | Occurrences and removal of pharmaceuticals and personal care products (PPCPs) in drinking water and water/sewage treatment plants: A review[END_REF]. The main emerging pollutants are pesticides, personal care products, pharmaceuticals and flame retardants [START_REF] Shahid | Current advances in treatment technologies for removal of emerging contaminants from water -A critical review[END_REF][START_REF] Sutherland | 15 years of research on wastewater treatment high rate algal ponds in New Zealand: discoveries and future directions[END_REF]. Conventional wastewater treatment technologies are not designed to remove or degrade these compounds [START_REF] Bai | Occurrence, distribution, and seasonality of emerging contaminants in urban watersheds[END_REF], while microalgae have been proven to be a powerful technology for bioremediation of emerging contaminants [START_REF] Lage | Algal Biomass from Wastewater and Flue Gases as a Source of Bioenergy[END_REF], [START_REF] Rempel | Cultivation of Microalgae in Media Added of Emergent Pollutants and Effect on Growth, Chemical Composition, and Use of Biomass to Enzymatic Hydrolysis[END_REF], [START_REF] Vale | Impact of parabens on microalgae bioremediation of wastewaters: A mechanistic study[END_REF]. Indeed, in the highly oxidizing environment that is achieved thanks to high oxygen levels, microalgae can catalytically degrade complex compounds, with efficiencies ranging from 30 to 80% for drugs such as ibuprofen, carbamazepine, and caffeine [START_REF] Ding | Biological removal of pharmaceuticals by Navicula sp. and biotransformation of bezafibrate[END_REF], [START_REF] Ding | Toxicity, degradation and metabolic fate of ibuprofen on freshwater diatom Navicula sp[END_REF], [START_REF] Matamoros | Assessment of the mechanisms involved in the removal of emerging contaminants by microalgae from wastewater: a laboratory scale study[END_REF], [START_REF] Matamoros | Capability of microalgaebased wastewater treatment systems to remove emerging organic contaminants: a pilot-scale study[END_REF]. Combined with bioadsorption and bioaccumulation, the algae process was found to achieve a remarkable removal efficiency for more than 50 emerging contaminants [START_REF] Rempel | Current advances in microalgae-based bioremediation and other technologies for emerging contaminants treatment[END_REF]. Adding a separation system, like membrane modules, is a technological breakthrough that leads to the uncoupling of Hydraulic Retention Time (HRT) and Solid Retention Time (SRT). On top of enhancing emerging contaminants removal [START_REF] Shahid | Current advances in treatment technologies for removal of emerging contaminants from water -A critical review[END_REF], HRT and SRT separation was shown to considerably enhance the overall process efficiency [START_REF] Robles | Performance of a membrane-coupled high-rate algal pond for urban wastewater treatment at demonstration scale[END_REF], [START_REF] Yunlong | and Henderson Simultaneous microalgae cultivation and wastewater treatment in submerged membrane photobioreactors: A review[END_REF] by increasing microalgal productivity by a factor 3.5 compared to a standard photobioreactor [START_REF] Marbelia | Membrane photobioreactors for integrated microalgae cultivation and nutrient remediation of membrane bioreactors effluent[END_REF] (see the synthesis of various literature studies in Table 1). The key idea is that algal biomass productivity results from the product of net growth rate and biomass concentration which is enhanced when increasing the biomass in the reactor while keeping growth rate at high levels. However, a separation system makes the process more flexible but also more complicated to operate and to optimize by increasing the number of parameters that can be tuned. Concentrating the algae in the system leads to a higher light attenuation and therefore, to a lower growth rate. In such a complex situation, a multi-parametric optimization must be carried out and the guidance of mathematical models has already proven very useful [START_REF] Sousa | Microalgae-based bioremediation of wastewaters -Influencing parameters and mathematical growth modelling[END_REF], [START_REF] Sánchez-Zurano | An Interactive Tool for Simulation of Biological Models Into the Wastewater Treatment With Microalgae[END_REF]. This is especially decisive in the case of complex dynamical bioprocesses involving a large range of interacting species, like for anaerobic digestion [START_REF] Pessoa | Numerical optimization based on generalized extremum seeking for fast methane production by a modified ADM1[END_REF].

Recently, the ALBA model (ALBA standing for ALgae-BActeria) was developed and validated. This mathematical model was built to simulate algal-bacterial interactions and competitions in HRABP, to evaluate bioremediation performances and to provide a better understanding of the interactions among the chemical, physical and biological processes. This model was validated and was proven able to accurately predict, under real outdoor conditions, the complex behaviour of two HRABP pilot-plants operated in two different locations and under different feeding conditions [START_REF] Casagli | A comprehensive growth model to optimize algae-bacteria wastewater treatment in raceway ponds[END_REF], [START_REF] Casagli | Balancing Microalgae and Nitrifiers for Wastewater Treatment: Can Inorganic Carbon Limitation Cause an Environmental Threat[END_REF]. The model turned out to be accurate in predicting the daily and seasonal responses, even during winter which was never effectively described by previously available models. In total, more than 630 days of pilot-scale experimental data validated the model performance [START_REF] Casagli | A comprehensive growth model to optimize algae-bacteria wastewater treatment in raceway ponds[END_REF], [START_REF] Casagli | Balancing Microalgae and Nitrifiers for Wastewater Treatment: Can Inorganic Carbon Limitation Cause an Environmental Threat[END_REF] under two different climatic conditions and with the same parameter set. The model confirmed, as experimentally observed [START_REF] Plouviez | Nitrous oxide emissions during microalgae-based wastewater treatment: current state of the art and implication for greenhouse gases budgeting[END_REF], that inorganic carbon could become limiting, as a result of the strong competition for inorganic carbon by algae and nitrifiers.

It was unexpected that this limitation still happened even under continuous CO2 injection; the ALBA model revealed that this was due to a shortage in alkalinity owing to the limited alkalinity/nitrogen ratio in the feed, insufficient to support full nitrification. At low alkalinity, inorganic carbon can no more be stored in the medium and becomes limiting both for algae and nitrifiers, leading to conditions favourable to N2O emission.

The high-fidelity ALBA model is therefore a powerful validated tool to assess the process productivity, the nitrogen removal rate, but also the risk of highly impacting N2O emission. In this work, the ALBA model was upgraded by implementing a solid/liquid separation system to investigate the complex interplay between SRT and HRT. The ALBA model was then validated in this new configuration using the experimental data of Robles et al. [START_REF] Robles | Performance of a membrane-coupled high-rate algal pond for urban wastewater treatment at demonstration scale[END_REF], carried out in an outdoor pilot reactor. The model was finally used to explore a large range of scenarios, varying HRT, SRT, and liquid depth. The trade-off between algal productivity and nitrogen recycling was explored, while focusing on operating modes that are not susceptible of leading to N2O emissions.

Finally, the simulation results are deeply discussed and compared with the experimental records reported in literature for outdoor algae-bacteria systems equipped with a separation system, comforting the model prediction and further quantifying the potential of separation systems for wastewater remediation, algae production and nitrogen recycling.

Material and methods

Brief recall of the ALBA model: structure, validation, and implementation

Simulations were run using the high fidelity ALBA model [START_REF] Casagli | A comprehensive growth model to optimize algae-bacteria wastewater treatment in raceway ponds[END_REF], simulating the dynamics of HRABP. The model considers a mixed culture of photoautotrophic algae, heterotrophic bacteria, Ammonium Oxidizing Bacteria (AOB) and Nitrite Oxidizing Bacteria (NOB). In total, it includes 19 biological processes and involves 17 state variables. The biokinetics equations are based on the Liebig's minimum law [START_REF] Baar | de Von Liebig's law of the minimum and plankton ecology (1899-1991)[END_REF] for limiting elements (carbon, nitrogen and phosphorus), meaning that the most limiting nutrient drives the overall kinetics. Moreover, biokinetics include dependences from light, temperature and pH, since they are the most influencing parameters on microbial dynamics (see SI.1).

The ALBA model embeds an in-depth description of the chemical reactions in the medium resulting in an accurate pH prediction. The pH sub-model is based on dissociation equilibria and mass balances of acids and bases, extending the approach proposed in the ADM1 (Anaerobic Digestion Model n.1 [START_REF] Batstone | The IWA Anaerobic Digestion Model No 1 (ADM1)[END_REF], [START_REF] Henze | van Activated Sludge Models ASM1, ASM2, ASM2d and ASM3[END_REF]).

The Total Alkalinity (TA) turns out to play an important role in the HRABP dynamics, since it determines the potential of inorganic carbon storage in the bulk. Its expression is defined in terms of molar quantities (mol m -3 ), as reported in Eq. 2.1.2 [START_REF] Casagli | Balancing Microalgae and Nitrifiers for Wastewater Treatment: Can Inorganic Carbon Limitation Cause an Environmental Threat[END_REF]:

TA = HCO 3 -+ 2CO 3 2-+ H 2 PO 4 -+ HPO 4 2-+ 2PO 4 3-+ OH -+NH 3 -H + -HNO 2 -HNO 3 -H 3 PO 4 2.1.2
The CO2, NH3, O2 gas/liquid transfer is also included, quantifying the rate through the global kLa and the diffusivity coefficients.

The model was previously calibrated and validated on two long-term datasets from two pilot-plants, operated under different climatic conditions. The first reactor was a demonstrative-scale raceway of 17 m 3 (with a 56 m 2 surface), located in the South of France (Narbonne area), fed on synthetic municipal wastewater [START_REF] Casagli | A comprehensive growth model to optimize algae-bacteria wastewater treatment in raceway ponds[END_REF]. The second reactor was a pilot-scale raceway of 1 m 3 (with a 3.8 m 2 surface), located in the North of Italy (Milan area), fed on the liquid fraction of piggery digestate [START_REF] Casagli | Balancing Microalgae and Nitrifiers for Wastewater Treatment: Can Inorganic Carbon Limitation Cause an Environmental Threat[END_REF]. In total, 30 days were used for the calibration phase, while 603 days were exploited for the model validation based on on-line measurements (dissolved oxygen, pH, temperature) and offline measurements (nitrogen compounds, algal biomass, total and volatile suspended solids, chemical oxygen demand). A more detailed description of the ALBA model, its calibration and validation procedure, can be found in Casagli et al. [START_REF] Casagli | A comprehensive growth model to optimize algae-bacteria wastewater treatment in raceway ponds[END_REF], and summarized in Supporting Information (SI.1 and SI.2).

Other models for describing algae-bacteria system are available in literature, as reviewed in Casagli et al. [START_REF] Casagli | Balancing Microalgae and Nitrifiers for Wastewater Treatment: Can Inorganic Carbon Limitation Cause an Environmental Threat[END_REF], they can also be modified similarly by introducing a solid/liquid separation term in the hydraulic balance.

The ALBA model was initially developed in AQUASIM and then implemented under MATLAB R2019b. Simulations were run on a PC with 8 i9 vPRO cores. Each simulation took approximately 3 minutes.

Case study

A theoretical case study was assumed for the optimization exercise, consisting in an industrial HRABP of 2000 m 2 treating digestate, located in the North of Italy (Milan area).

Simulations were run focusing on the spring season, where the most affecting environmental conditions (i.e., light and temperature) are closer to the optima for algal growth. The influent was assumed to be a diluted liquid fraction of digestate with the characteristics reported in Table 2 (taken from Pizzera et al., [START_REF] Pizzera | Digestate treatment with algae-bacteria consortia: A field pilot-scale experimentation in a suboptimal climate area[END_REF]). Note: all the biomasses (XALG, XAOB, XNOB, XH) were considered as zero in the influent.

Validated model parameters were taken from Casagli et al. [START_REF] Casagli | A comprehensive growth model to optimize algae-bacteria wastewater treatment in raceway ponds[END_REF]. Only the kLa value was assumed to be lower (10 day -1 ), which is more similar to an industrial scale application.

In all the simulations, the pH was controlled at 7.5 with CO2 bubbling (see Supp Info SI.3).

Implementation of HRT and SRT decoupling

The separation system was implemented by imposing a retention factor on all the particulate variables:

Ret = SRT-HRT SRT = 1- HRT SRT = 1-α 2.3.1
This factor can be positive or negative, depending on the value resulting from the ratio α = HRT/SRT. Consequently, the evolution in time of the particulate variables (Xj) that are retained by the separation system (dXj/dt [gXj m 3 d -1 ]) was implemented in the model as follows:

dX j dt = Q IN V •X j, IN -α•X j • ( Q IN -Q EVAP V ) ± ∑ v i,j ρ i i 2.3.2
where QIN is the imposed inflow rate [m 3 d -1 ]; QEVAP is the outflow leaving the system by evaporation, [m 3 d -1 ]; V is the reactor volume, [m 3 ]; Xj,IN is the concentration of Xj entering the system with the influent [g m -3 ]; Xj is the concentration in the system [g m - 3 ]; νi,j is the stoichiometric coefficient associated to the state variable Xj and the process i and ρi is the rate of process i [gXj m -3 d -1 ].

Optimization study

For the optimization study, 19200 simulations were run, testing the simultaneous influence of different factors on the HRABP performances: i) the HRT (0.2 -10 d); ii) the SRT (0.2 -10 d); iii) the liquid depth (δL 0.06 -1 m); and finally, the supplementation of Total Alkalinity (TA, 0 -12.5 mol m -3 added) through dosage of CaCO3.

Simulations were run considering that HRT and SRT could be independently tuned, simulating therefore cases with SRT lower than HRT (biomass harvesting with liquid recirculated in the reactor) or, reversely, HRT lower than SRT (biomass retention).

Algal productivity [gAlg m -2 d -1 ] was defined as:

ALG productivity ≔ X ALG • Q out • α S • φ 2.4.1
Where XALG is the algae concentration in the pond [gCOD m -3 ]; QOUT is the outflow, defined as the algebraic sum of QIN-QEVAP [m 3 d -1 ]; S is the reactor surface [m 2 ]; φ=0.64 is the conversion factor from CODALG to dry weight (DW) of algae (gDW-Alg.gCOD -1 ).

The efficiency of nitrogen recycling rate was defined as:

Eff N RECYC [%]= Flux N OUT ALG Flux N IN •100 2.4.2
Where: FluxNout,ALG is the mass flow of Nitrogen [gN d -1 ] in the harvested algal biomass and FluxN,IN is the mass flow of Nitrogen [gN d -1 ] fed to the HRABP:

Flux N OUT ALG [ gN d ] = X ALG •iN BM ALG •Q OUT • α 2.4.3 Flux N IN = Q IN (S NH IN +S NO2 IN +S NO3 IN +S S IN •IN SS +S I IN •IN SI +X S IN •iN XS +X I IN •iN XI ) 2.4.4
where iNBM,ALG is the N content in the algal biomass [gN gCOD -1 The N recycling efficiency, and the algal biomass productivity were used as Key Performance Indicators.

In addition, the N2O emission risk-factor was considered to spot unsuitable operation modes. Finally, a closer look at the biomass distribution among algae and bacteria biomasses was helpful to better interpret the evolution of the ecosystem.

Results and discussion

Model accuracy in predicting nitrogen removal rate and algal productivity

The model prediction capability for the most relevant monitored variables (pH, SO2, SNH, SNO3, SNO2 and XALG) can be found in Casagli et al. [START_REF] Casagli | A comprehensive growth model to optimize algae-bacteria wastewater treatment in raceway ponds[END_REF], [START_REF] Casagli | Balancing Microalgae and Nitrifiers for Wastewater Treatment: Can Inorganic Carbon Limitation Cause an Environmental Threat[END_REF] (see also a summary in Supporting Information SI.1 and SI.2) and it covers an uncommonly broad range of conditions, including all the seasons, and HRT (=SRT) ranging from 2.6 to 20 days (see

S.I.7).

The total ammoniacal nitrogen (TAN) removal rate (gN m -2 d -1 ) and the Total Suspended Solids (TSS) productivity (gTSS m -2 d -1 ) were also computed from the experimental measurements and compared to the model predictions (Fig. 2 The generally higher values obtained in Narbonne along all the seasons are mainly due to the lower HRT, i.e. 5 days, (not considering the period where it was varied between 2.6 and 10, for days 240-253), compared to Milan (HRT = 10 -20 days), corresponding to a lower biomass harvesting rate. The TSS productivities shown in Figure 1a,b are in agreement with the experimental values reported in other literature works, where the raceways were operated under similar conditions. For instance, 10±6.5 gTSS m -2 d -1

[32]; 15.2±2.1 gTSS m -2 d -1 [START_REF] Morillas-España | Improvement of wastewater treatment capacity using the microalga Scenedesmus sp. and membrane bioreactors[END_REF]; 5.5 ± 7.4 g TSS m -2 day -1 [START_REF] Mantovani | Outdoor pilot-scale raceway as a microalgae-bacteria sidestream treatment in a WWTP[END_REF] ; 12.3 gTSS m -2 d -1 as maximum value during winter and 32.4 gTSS m -2 d -1 as maximum value during summer [START_REF] Sánchez-Zurano | Annual assessment of the wastewater treatment capacity of the microalga Scenedesmus almeriensis and optimisation of operational conditions[END_REF].

The model was further used to estimate the contributions from the different biomasses (XALG, XAOB, XNOB, XH) and from inert and biodegradable compounds (XI and XS respectively) which all together make up the TSS.

Algae productivity in Milan represented 72% of the TSS productivity on a yearly basis (i.e. 7.03 gALG m -2 day -1 compared to 9.8 gTSS m -2 day -1 ), corresponding to a microbial biomass (XBIO = XALG+XAOB+ XNOB+XH) composed by 90.4% of algae on yearly average. For the Narbonne case study, algae productivity was 63% of the TSS productivity on a yearly bases (i.e. 9.5 gALG m -2 day -1 compared to 15.1 gTSS m -2 day -1 ), corresponding to a microbial biomass where algae represented 76.6% of the community on yearly average.

Looking at TAN removal rate in Narbonne (Fig. 3), it was on average 0.4 ± 0.1 gN m -2 d -1 in summer, with HRT of 5 days. Peaks of 0.9 gN m -2 d -1 were reached at HRT of 10 days, while the lowest value was achieved for a HRT of 2.6 days (0.2 gN m -2 d -1 ). In autumn and winter, the TAN removal rate was generally lower than 0.5 gN m -2 d -1 , due to the sub-optimal growth conditions for algae and bacteria. For the raceway located in Milan, the TAN removal rate was higher during summer, being in a range of 4.6-10.2 gN m -2 d -1 (average 6.9 ± 1.

Model accuracy when uncoupling HRT and SRT

The experimental campaign of Robles et al. [START_REF] Robles | Performance of a membrane-coupled high-rate algal pond for urban wastewater treatment at demonstration scale[END_REF] , was used for validating the ALBA model when HRT is uncoupled from SRT. The experimental campaign was performed on the same outdoor pilot-scale raceway located in Narbonne, fed with synthetic municipal wastewater [START_REF] Nopens | Stability analysis of a synthetic municipal wastewater[END_REF], previously used for the calibration and validation of the ALBA model [START_REF] Casagli | A comprehensive growth model to optimize algae-bacteria wastewater treatment in raceway ponds[END_REF]. The HRABP was operated with a liquid depth of 0.3 m and connected to a high-flow industrial-scale membrane compartment. The authors evaluated the process performances computing the average biomass productivity and nitrogen removal rate over a 20 days experimental period, for each operational condition tested. Uncoupling SRT and HRT by membrane filtration improved process efficiency, with higher biomass throughput and nutrient removal rate at lower HRT operation. At an SRT of 6 days, the biomass productivity increased up to 19.9 ± 0.5 and to 28.5 ± 0.5 gVSS m -2 d -1 when HRT was set to 4 and 2.5 days, respectively. The corresponding nitrogen removal rates were 2.3 ± 0.4 and 3.3 ± 0.4 gN m -2 d -1 , respectively. The authors ran also an experiment with HRT set to 2.5 and SRT set to 6 days, but under sub-optimal light and temperature conditions, obtaining lower values for both biomass productivity and nitrogen removal rate (19.6 gVSS m -2 d -1 and 2.4 g N m -3 d -1 , respectively), compared to the experiment carried out with the same operational conditions and under optimal light and temperature ranges.

The ALBA model was run under the same operational conditions, considering the climatology of the year 2018-2019 [START_REF] Casagli | A comprehensive growth model to optimize algae-bacteria wastewater treatment in raceway ponds[END_REF]. Biomass productivity and nitrogen removal rate were computed according to Robles et al., [START_REF] Robles | Performance of a membrane-coupled high-rate algal pond for urban wastewater treatment at demonstration scale[END_REF] and simulation results were averaged on a seasonal time scale. Spring 2018 presented similar average light and temperature compared to the one reported by Robles et al. [START_REF] Robles | Performance of a membrane-coupled high-rate algal pond for urban wastewater treatment at demonstration scale[END_REF] as being optimal. The autumn climatology of 2018 turned out to be comparable to the light and temperature conditions reported as suboptimal by Robles et al. [START_REF] Robles | Performance of a membrane-coupled high-rate algal pond for urban wastewater treatment at demonstration scale[END_REF]. The model predictions, for both biomass productivity and nitrogen removal rate, are perfectly in agreement with the experimental results (Figure 3), with a R 2 of 0.78 and 0.93, respectively. Discrepancies between simulated and experimental data are mainly due to the dynamic temperature and light dataset used for simulations, that are not belonging to the same year as in Robles et al. [START_REF] Robles | Performance of a membrane-coupled high-rate algal pond for urban wastewater treatment at demonstration scale[END_REF], since these data were not available. Moreover, while in the work of Robles et al.

[19], the measurements were averaged on 20 days of experiments for each tested condition, the simulated data were averaged on the entire season (3 months for spring and 3 months for autumn).

These

Best strategies to avoid GHG emissions: the key role of alkalinity

Simulations were run to explore optimal SRT and HRT for treating digestate. The (constant) liquid depth was set to 0.225 m, typical of raceways systems. The TA in the influent was set at its nominal value (10 mol m -3 ), typical of a diluted digestate.

Simulation outcomes are shown in Figure 4, while the 2D version of the graphs are reported in SI.5. As it was shown in Casagli et al. [START_REF] Casagli | Balancing Microalgae and Nitrifiers for Wastewater Treatment: Can Inorganic Carbon Limitation Cause an Environmental Threat[END_REF], despite pH regulation, a strong competition for inorganic carbon between algae and nitrifiers can take place since alkalinity is strongly reduced by nitrification. The consequence is a risk of N2O emission that can be very detrimental to the environment. In fact, it turns out that the N2O risk factor becomes nonzero for couples of HRT ≥ 4 d and SRT ≥ 2 d, reaching the highest values (35-45%) for HRT ≥ 8 d and SRT ≥ 4 d (Fig. 4C).

The maximum of algal productivity (20.18 gAlg m -2 d -1 ) was obtained for HRT=0.2 d and SRT=2 d (Fig. 4A). This pair of hydraulic and solid retention times is not risky as for N2O emission (risk factor is 0%). However, looking at the algae concentration and the efficiency in nitrogen recycling, both values are low (280 gCOD m -3 and 0.34% respectively, see Fig. 4B and4D). Indeed, for this low HRT, the nitrogen loading rate is high, leading to low efficiency of nitrogen recycling (EffN,Recyc). Very low SRT, typically lower than 1 d, whatever the HRT, lead to algal washout, with algal biomass concentration that goes to zero and so does the algal biomass productivity (Fig. 4A and4D). The maximum value for the efficiency in nitrogen recycling (10.4%) is obtained for a very different operational regime, i.e. HRT=10 d and SRT=3 d (Fig. 4B). However, these conditions are associated to an unsustainable N2O risk factor of 30% (Fig. 4C).

The higher algae and nitrifying bacteria concentrations, indeed lead to the exhaustion of the inorganic carbon pool. In addition, the algal productivity related to the highest efficiency in nitrogen recycling is reduced to 12.4 gAlg m -2 d -1 . This opposite trend reflects a classical contrasting tendency between productivity and efficiency in nutrient use, like for the nitrogen removal rate and removal efficiency in wastewater remediation.

This explains the conflicting trend for the Key Process Indicators (2.4.1 and 2.4.2).

Finally, the best operational conditions in terms of HRT and SRT for reaching a better trade-off between algae productivity and nitrogen recycling, while avoiding the risk of N2O emission, appear to be when HRT ranges from 2 to 5 d and SRT from 2 to 3 d.

However, under these conditions, the efficiency of nitrogen recycling is never higher than 6.6% (see Fig. 4A and4B).

In order to provide a better picture of the gain achievable by uncoupling HRT and SRT, Figure 5 reports a comparison between the algal productivity and the efficiency in nitrogen recycling when no solid/liquid separation is applied (i.e. HRT=SRT). In this case, HRT is adapted according to SRT, in order to maximize either algal productivity or nitrogen recycling efficiency. Figure 5A clearly shows that higher algae productivity are achieved for all the SRT values between 2 and 10 days, when operating the system with an HRT of 0.2 days. A doubling in algal productivity can be obtained by appropriately modulating HRT compared to the reference case of SRT=HRT. Figure 5B shows that the same is valid for the nitrogen recycling efficiency, that reaches higher values when these two operational parameters are decoupled (imposing HRT 10 days at varying SRT between 2 and 10 days). Nitrogen recycling efficiency could be tripled compared to the standard case when HRT=SRT.

Additional simulations were also run, by increasing the influent TA from 10 up to 15, 17.5, 20 and 22.5 mol m -3 (see Table SI TA supplementation appears to be particularly beneficial for algae-bacteria systems when the aim is to treat wastewater, since it allows to set higher HRT and SRT, favouring: i) slow growers, such as AOB and NOB; ii) the inorganic carbon availability to guarantee the full nitrification process, simultaneously avoiding the competition with algae; iii) reduced NH3 stripping conditions due to the improved TAN removal (see Fig. 7D1 and7D2) and, thus, lower TAN concentration in the reactor. The environmental benefit of this strategy is therefore evident.

It must be highlighted that alkalinity addition should be adjusted carefully to avoid a marked increase in CO2 emissions due to overloading the system with inorganic carbon, especially when the pH in the raceway is controlled with CO2 bubbling. Ideally, it should be regulated according to the CO2 saturation level in the pond, that is strictly related to water temperature. In fact, the driving-force in the gas-liquid exchange is given by the difference between the gas saturation level in the liquid phase and the actual concentration of the dissolved gas. With TA addition, the buffering capacity of the system increases, and so does the level of soluble inorganic carbon. The risk is that when an external source of CO2 is injected in the system for pH regulation, most of it is lost to the atmosphere.

When considering simultaneously the algal biomass productivity and the nitrogen recycling optimization, the situation does not change significantly. Regulating TA avoids N2O emissions, but it does not allow to simultaneously reach these two conflicting objectives.

Unravelling the influence of liquid depth

Two additional sets of simulations where run, considering the liquid depths of 0.12 and 0.06 m, first with the nominal alkalinity in the influent (10 mol m -3 ). Results are shown in Figure 6.

For δL=0.12 m, the maximum of algal productivity is 19.4 gAlg m -2 d -1 for HRT=0.2 d and SRT=2 d (Fig. 6A), which is slightly lower than for the nominal depth (δL=0.225m) for the same HRT and SRT. The corresponding concentration of nitrogen recycled in the algal biomass is doubled (0.6%, Fig. 6B). Indeed, at constant HRT, the incoming nitrogen load per surface unit is lower for lower depth, while the incoming nitrogen load per volume unit is constant. The lower areal productivity, even if the algal concentration in the reactor is higher (Fig. 7 A1), is due to the lower flowrate QOUT, necessarily associated to a lower δL for the same HRT, eventually resulting in a lower algal effluent flow rate. For this pair of HRT-SRT values, there is no risk of N2O emission, but the nitrogen recycling efficiency is still very low.

The maximum nitrogen recovered as algae biomass (11.8%) is obtained for HRT=10 d and SRT=3 d, but for this pair of values the N2O risk factor is 35% (Fig. 6B and6C). So far, the best working range of both HRT and SRT, for guarantying elevated algal biomass productivity and efficiency in nitrogen recycling is 2-5 d, similarly to the case δL=0.225m. In comparison, the nitrogen uptake from the algal biomass is higher and the best trade-off between algal biomass productivity (15.25 gAlg m -2 d -1 ) and nitrogen recycling (7.03%), avoiding N2O favourable conditions, is obtained for HRT=3 d and SRT=5 d.

Decreasing δL to 0.06 m (assumed to be the lowest operational depth for this kind of systems), the maximum of algal biomass productivity (15.2 gAlg m -2 d -1 ) is obtained for HRT=1 d and SRT=6 d (Fig. 6D), with a corresponding nitrogen recycling of 4.7% (Fig. 6E) and algae concentration in the reactor of 2603 gCOD m -3 (see Fig. 9A1). Compared to the the two-fold depth (δL=0.12 m), the maximum productivity is lower, even if the algal concentration is almost 5 times higher. This is due to a much lower nutrient influent load for 6 cm of liquid depth. Correspondingly, the nitrogen uptake in the algal biomass is higher. In addition, under this condition of HRT and SRT, the N2O risk factor is 0% (Fig. 6F). The best result for nitrogen recycling efficiency (12.4%) was found for HRT=10 d and SRT=3 d, but it corresponds to a N2O risk factor of 45% (Fig. 6F).

Finally, for this low liquid depth competition for inorganic carbon is stronger. It results that the N2O emission risk appears for a wider range of cases compared to δL=0.12 m.

In fact, while every value between 2 and 10 d can be applied for the SRT, the HRT should stay lower than 1 d, or being increased: i) up to 7, imposing SRT=2 d; ii) up to 3, imposing SRT=3 d; iii) up to 2, imposing SRT=4-5 d.

Additional simulations were run increasing the influent TA (see Figure 7 and8), in order to see the coupled effect of working at low liquid depth without alkalinity shortage in the reactor.

Indeed, the TA influent addition allows to explore a wider range of HRT and SRT (proportionally to the alkalinity added) without operating under inorganic carbon limitation. This allows finding a much better trade-off between algal productivity and nitrogen recycling efficiency. Total alkalinity must be at least 20 mol m -3 , to guaranty that the N2O risk factor is zero (Fig. 7C) for all HRT and SRT values.

Figure 8A and B show more explicitly which is the advantage for algal production and nitrogen recycling respectively, when decoupling HRT and SRT, compared to the scenario without membrane separation application. Adapting the HRT for each SRT value can lead to a threefold increase in algal productivity or in the nitrogen recycling rate.

Analysis of the optimal regimes

Figure 9 summarizes those conditions maximising algae productivity (left column) and nitrogen recycling (right column), for three liquid depths (δL: 0.22, 0.12 and 0.06 m).

Algae and heterotrophic bacteria are indeed favoured (see Fig. 9 A1) for the HRT-SRT pair that maximise the algae productivity (HRT 0.2 d, SRT 2d, δL 0.22 and 0.12 m), while nitrifying bacteria are washed-out from the system. Nonetheless, even if algae productivity reaches higher values, between 18 and 21 gAlg m -2 d -1 (Fig. 9 B1), the corresponding nitrogen recycling remains below 1% (Fig. 9 C1). The biomass composition drastically changes for HRT 1d, SRT 6 d and δL 0.06 m. Algae are more concentrated (with a concentration factor, compared to δL 0.12 and 0.22 m, ranging between 5 and 9.5, respectively), and nitrifying bacteria reach concentrations comparable to those of heterotrophic bacteria (Fig. 9 A1). In this case, the corresponding algae productivity is slightly lower (15 gAlg m -2 d -1 ), due to the decreased harvesting rate. However, the associated nitrogen recycling is higher (4.7%, see Fig. 9 C1), suggesting this operational configuration as the best trade-off between algae productivity and nitrogen recycling when the target objective is the algae productivity.

The beneficial effect of increasing the influent alkalinity is mostly evident for TAN removal rate (Fig. 9 D1), that is always higher than 10 gN m -2 d -1 (reaching 22.6 gN m -2 d -1 when operating at δL 0.22 m), compared to the cases where alkalinity was not added.

Also looking at the TAN removal rate in Fig. 9 D1, the best trade-off when maximizing algae productivity remains to operate the system at 0.06 m depth, HRT 1 d and SRT 3 d, since there is not a substantial difference comparing the TAN removal rate at 0.12 or 0.06 m.

The biomass composition reported in Fig. 9 A2, shows how the HRT-SRT pair maximizing the nitrogen recycling in the algal biomass (10 and 3 d, respectively) favour algae and nitrifying bacteria, while heterotrophic bacteria concentration remains low.

The highest algae productivity (14.3 gAlg m -2 d -1 , see Fig. 9 B2) is obtained with alkalinity addition for δL 0.22 m, also corresponding to the highest TAN removal rate (6.9 gN m -2 d -1 , Fig. 9 D2), while the best nitrogen recycling is reached for δL 0.06 m and with TA addition (Fig. 9C2). Results reported in Fig. 9 A2, B2, C2 and D2 suggest that the best trade-off between algae productivity and nitrogen recycling, when the objective target is the nitrogen recycling, by operating the system at δL 0.22 m and with TA influent addition. Therefore, by analysing in detail the cluster of conditions maximizing algae productivity and nitrogen recycling among the tested conditions, we can conclude that operating with HRT≤1 d and SRT≤6 d allows reaching the highest values of algae production, while operating at HRT 10 d and SRT 3 d allows maximising the nitrogen recovery in the algal biomass (up to 15%). Given the opposite trends of these two targets, the best trade-off is to run the system at lower liquid depth (0.06 m) when the target is maximising the algal productivity, while operating with higher liquid depth (0.12-22 m) when the objective is to maximize the nitrogen recycled in the algae.

Moreover, results highlighted that working at lower liquid depth allows: i) to maintain a higher algal biomass concentration in the reactor (requiring less energy for algae harvesting), with comparable biomass productivity; ii) to prevent massive growth of both heterotrophic and nitrifying bacteria, that can compete with algae for the main nutrients [START_REF] Chiaramonti | Review of energy balance in raceway ponds for microalgae cultivation: Re-thinking a traditional system is possible[END_REF].

Further step to target emerging contaminants

Microalgae can play a role in the treatment of emerging contaminants by bioadsorption, bioaccumulation and metabolic degradation [START_REF] Rempel | Current advances in microalgae-based bioremediation and other technologies for emerging contaminants treatment[END_REF]. There is currently no efficient mathematical model validated in outdoor conditions describing such processes. It is even likely that models should be molecule-specific to be accurate.

However, the ALBA model predicts the oxygen concentration in the medium, which gives an index of the intensity of the oxidative stress influencing the molecular degradation [START_REF] Matamoros | Capability of microalgaebased wastewater treatment systems to remove emerging organic contaminants: a pilot-scale study[END_REF]. Simulation results showed that, under the tested conditions, oxygen concentration could reach pics up to 22 mgO2.L -1 during the day. A modelling study can be further carried out, imposing as target objective to maximise the peaks of oxygen concentration in the reactor, also by exploring a wider range of HRT and SRT values. The model can therefore indirectly predict the benefit in terms of emerging contaminant removal due to oxidative stress. Further work must be carried out to quantify these aspects, and further integrate it with an estimation of the gain in emerging contaminants treatment with a separation system. This is an important point when computing the cost-benefit of the separation system.

As recently reviewed by You et al. [START_REF] You | Pollutants affect algae-bacteria interactions: A critical review[END_REF], pollutants (including pesticides, metals, engineered nanomaterials, pharmaceutical and personal care products, and aromatic pollutants) affect both bacteria and algae by interfering with their relationships. Cell-to-cell adhesion, substrate exchange and biodegradation of organic pollutants, enhancement of signal transduction, and horizontal transfer of tolerance genes are defence strategies in algal-bacterial systems, which should be further tamed and applied in wastewater treatment systems. Enhancing the treatment capacity by continuous addition of more efficient algae and/or bacteria species is then the next step. This complementary strategy of bioaugmentation could be combined to the separation system to further enhance the process performance. Numeric simulations can guide the process optimization and development for this new and promising strategy.

Long-term membrane application for algae bacteria systems

Membranes allow to decouple HRT and SRT and thus enhance the process performance by more than 20%, with the possibility of recovering the algal biomass for valuable applications [START_REF] Sutherland | Microalgal bioremediation of emerging contaminants -Opportunities and challenges[END_REF].

However, membrane biofouling due to heavy suspended solids load, the small size of microalgae and their poor sedimentation properties, combined with bacteria and exopolymeric substances can lead to a progressive decline in the permeate flux and an increase in the energy demand for operating the process, thus limiting long term operations [START_REF] Arbib | Optimization of pilot high rate algal ponds for simultaneous nutrient removal and lipids production[END_REF], [START_REF] Parsy | Impact of seasonal variations on Nannochloropsis oculata phototrophic productivity in an outdoor pilot scale raceway[END_REF]. Most of the membrane designs were developed for wastewater treatment [START_REF] Chiaramonti | Review of energy balance in raceway ponds for microalgae cultivation: Re-thinking a traditional system is possible[END_REF], [START_REF] You | Pollutants affect algae-bacteria interactions: A critical review[END_REF], with very different objectives, such as removal of pollutants and recovery of clean permeate vs recovery of valuable algal biomass [START_REF] Parsy | Impact of seasonal variations on Nannochloropsis oculata phototrophic productivity in an outdoor pilot scale raceway[END_REF], [START_REF] Novoa | Membrane Fouling in Algal Separation Processes: A Review of Influencing Factors and Mechanisms[END_REF]. The traditional fouling control techniques could be inefficient when targeting a by-product recovery [START_REF] Novoa | Membrane Fouling in Algal Separation Processes: A Review of Influencing Factors and Mechanisms[END_REF] and new directions have emerged in recent years, while algal production was coupled to wastewater treatment [START_REF] Liao | A review of membrane fouling and its control in algal-related membrane processes[END_REF]. Different fouling control techniques have been developed recently, supported by a better understanding of the interactions between the foulants and membrane surfaces.

The main strategies consist in changing the operating conditions, in pre-treating the feed (using coagulation, adsorption, or oxidation), in mechanically cleaning of the membranes, and in modifying the membrane surfaces by coatings or blending with polymers and nanoparticles [START_REF] Novoa | Membrane Fouling in Algal Separation Processes: A Review of Influencing Factors and Mechanisms[END_REF]. Choosing the adequate membrane cleaning mechanism is of major importance in terms of energy consumption, permeate production, and overall process performance [START_REF] Liao | A review of membrane fouling and its control in algal-related membrane processes[END_REF], [START_REF] Zhang | Membrane fouling in aerobic granular sludge (AGS)-membrane bioreactor (MBR): Effect of AGS size[END_REF]. The development of membranes for algae-bacteria systems with excellent bioactivity and associated to antifouling strategies is of crucial importance for long term applications.

Conclusion

This study showed the complexity and challenges for optimizing outdoor HRABP featured with solid/liquid separations systems, pointing out that the efficiency in nitrogen recycling and the algal productivity cannot be maximized simultaneously.

Thus, it is fundamental to first choose the target of the process (e.g. wastewater remediation, algal productivity, emerging contaminant treatment…), which will then be reached by decoupling HRT and SRT: a gain of 23% on algal productivity (up to 20.3 g.m -2 .day -1 ), and of 35% in nitrogen recycling could be obtained for the specific case study. When combined with depth optimization and alkalinity addition, nitrogen recycling rate can reach 15%, with an algal productivity staying above 16 g.m -2 .day -1 thus offering an interesting trade-off. Long-term application of efficient membrane modules, and the characterisation of their efficiency on emerging contaminant treatment, will then be crucial for future implementation of these strategies.

Figure 1 .Figure 2 .

 12 Figure 1. Total Suspended Solids productivity experimentally measured and computed with the ALBA model for Milan and Narbonne case studies according to the season (a: summer; b: autumn and winter).

4. 5 dBFigure 3 .

 53 Figure 3. Comparison between experimental [19] and simulation results (this study) for HRT uncoupled from SRT.

C: δL 22

 22 .5 cm; TA: 10 mol m -3 D: δL 22.5 cm; TA: 10 mol m -3

Figure 4 .Figure 5 .

 45 Figure 4. Simulation results at varying HRT (0.2-10 d) and SRT (0.2-10 d). Influent TA is 10 mol m -3 and depth, δL, is 0.225 m. A: algal biomass productivity [gAlg m -2 d -1 ]; B: Nitrogen recycling efficiency in the algal biomass [%]; C: N2O risk factor (percentage of time along the day for which N2O formation conditions occur, i.e. inorganic carbon < 0.2 molC m -3 ), [%]; D: Algal biomass concentration in the reactor [gCOD m -3 ]. The red line represents the scenario with no solid/liquid separation (i.e. HRT = SRT).

D: δL 6

 6 cm; TA: 10 mol m -3 E: δL 6 cm; TA: 10 mol m-3 F: δL 6 cm; TA: 10 mol m -3

Figure 6 . 3 Figure 7 .Figure 8 .Figure 9 .

 63789 Figure 6. Simulation results at varying HRT (0.2-10 d) and SRT (0.2-10 d); influent TA (10 mol m -3 ) and δL (0.12 m, A, B and C; 0.06 m, D, E and F). A, D: algal biomass productivity [gAlg m -2 d -1 ]; B, E: Nitrogen recycling efficiency in the algal biomass [%]; C, F: N2O risk factor (percentage of time along the day for which N2O formation conditions occur, i.e. inorganic carbon < 0.2 molC m -3 ), [%]. The red line represents the scenario with no solid/liquid separation (i.e. HRT = SRT).

Table 1 :

 1 Biomass productivity and nitrogen removarate reported in literature for various outdoor photobioreactors, different influent, HRT, SRT and weather conditions.

	95
	Wastewater characteristics (Av ± St.Dev)	Experimental set-up

**Total phosphorous value

Table 2 :

 2 Influent characteristics and associated model state variables.

Influent characteristics used for the simulations: liquid fraction of piggery digestate Description Symbol Value Unit

	Cations	S cat	10	mol m	-3
	Anions	S an	1e-005	mol m	-3
	pH	pH	8.5	-
	Total alkalinity	TA	500	mgCaCO3 L -1
	Total ammoniacal nitrogen	S NH	310.4	gN m	-3
	Nitrate	S NO3	11.7	gN m	-3
	Nitrite	S NO2	0	gN m	-3
	Inorganic soluble orthophosphates	S PO4	13.6	gP m	-3
	Inorganic carbon	S IC	266.1	gC m	-3
	Readily biodegradable organic matter	S S	133.3	gCOD m	-3
	Inert soluble organic matter	S I	247.6	gCOD m	-3
	Inert particulate organic matter	X I	118.9	gCOD m	-3
	Slowly biodegradable organic matter	X S	13.3	gCOD m	-3
	Dissolved oxygen	S O2	8	gO2 m	-3

  ], QIN is the inflow rate [m 3 d -1 ]; SNH,IN, SNO2,IN and SNO3,IN are the influent concentrations of ammoniacal nitrogen, nitrite and nitrate, respectively [gN m -3 ]; SS,IN, SI,IN, XS,IN and XI,IN are influent

	concentration of the state variables representing the organic matter concentration
	[gCOD m -3 ] as soluble degradable, soluble inert, particulate degradable, and particulate
	inert (I) components, respectively; iNSS, iNSI, iNXS and iNXI, [gN gCOD -1 ] quantify their
	nitrogen content.

  productivity is rarely experimentally measured in algae-bacteria consortia, since only TSS or VSS (Volatile Suspended Solids) are monitored. Both turned out to be accurately predicted, as shown in Figure1and 2, with p-values for TSS

	productivity and TAN removal rate predictions below 1E-6. The highest TSS productivity
	was reached in summer, for both Narbonne and Milan, with an average value of 16.0 ±
	5.6 [gTSS m -2 d -1 ] (Summer 2018, Narbonne), 18.2 ± 2.0 [gTSS m -2 d -1 ] (Summer 2019,
	Narbonne) and 12.2 ± 3.5 [gTSS m -2 d -1 ] (Summer 2016, Milan).
	During autumn and winter, the TSS productivity decreased in both experimentations,

and 3 for summer, autumn and winter, while spring results are shown in Supporting information, SI.4). It is worth mentioning that algal with an average value of 12.4 ± 3.1 and 12.7 ± 2.4 [gTSS m -2 d -1 ] in autumn and winter, respectively, in Narbonne and of 6.1 ± 3.1 [gTSS m -2 d -1 ] in autumn in Milan.

  .6.1 in SI.6). For an influent alkalinity increased up to 15 and 17.5 mol m -3 , inorganic carbon limiting conditions can still be reached when applying HRT ≥ 6 d and SRT ≥ 4 d, with a maximum value of N2O risk factor of 45% for HRT=10 d and SRT=5, 6 and 7 d. On the contrary, with 20 and 22.5 mol m -3 of influent TA, the N2O risk factor is always zero, indicating that this level of alkalinity (0.90 and 1.02 mol molNH4 +-1 respectively, corresponding to 2.51 and 2.82 gCaCO3 gNH4 +-1 ) is adequate to fully support both algae and nitrifiers carbon request.

A) Biomass productivity, in gVSS m -2 d -1 ; B) Nitrogen removal rate, in gN m -2 d -1 .

Acknowledgments

OB, FC and FB benefited from the support of the ADEME Biomsa project. EF wish to thank the Fondazione Cariplo (project: Polo delle Microalghe) for their financial support.

Author contribution statements

FC developed the modified model with feedback from OB. FC and OB implemented the simulation framework. FC wrote the manuscript with permanent feedback from OB.

EF, FB and OB guided FC in the definition of the case study framework. FB and OB obtained grants to support the study and supervised the work. In general, all authors provided critical feedback and helped shape the research, analysis and manuscript.